Télécharger exdiar.eso

Retour à la liste

Numérotation des lignes :

  1. C EXDIAR SOURCE BP208322 15/06/22 21:18:03 8543
  2. SUBROUTINE EXDIAR(MRIGID,MCHPOI)
  3. IMPLICIT REAL*8 (A-H,O-Z)
  4. IMPLICIT INTEGER (I-N)
  5. C*************************************************************************
  6. C Operateur EXDIAR
  7. C
  8. C OBJET : Extrait la diagonale d'une matrice.
  9. C
  10. C***********************************************************************
  11. C HISTORIQUE : 06/06/1 : Première version
  12. C
  13. C***********************************************************************
  14. -INC CCOPTIO
  15. -INC CCHAMP
  16. -INC SMCHPOI
  17. -INC SMRIGID
  18. -INC SMELEME
  19. -INC SMCOORD
  20. -INC TMTRAV
  21. *
  22. segment trav1
  23. character*4 compd(nbincd),compdd(nbincd)
  24. endsegment
  25. segment trav2
  26. integer cor1p(nligrp),cor1d(nligrd)
  27. endsegment
  28. *
  29. * Pour chaque segment DESCR, on a un segment CORES2 où
  30. * CORDIA donne les couples (NLIGRP,NLIGRD) qui sont des
  31. * termes diagonaux avec les cas particuliers suivant :
  32. * CORES2=0 : il n'y a pas de termes diagonaux dans ces matrices
  33. * élémentaires
  34. * CORES2=-1 : les termes diagonaux de la matrices élémentaire (carrée)
  35. * sont des termes diagonaux de la matrice globale
  36. * on a aussi a segment CORES3 où COMTRA indique
  37. * le nom d'inconnue primal+numéro d'harmonique suivant
  38. * le repérage induit par MTRA1 et MTRA2
  39. *
  40. SEGMENT CORES1
  41. INTEGER ITRAV2(NRIGEL)
  42. INTEGER IPCOR2(NRIGEL)
  43. INTEGER IPCOR3(NRIGEL)
  44. ENDSEGMENT
  45. SEGMENT CORES2
  46. INTEGER CORDIA(2,NLIGP)
  47. ENDSEGMENT
  48. SEGMENT CORES3
  49. INTEGER COMTRA(NLIGP)
  50. ENDSEGMENT
  51. *
  52. SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
  53. SEGMENT MTRA1
  54. * CHARACTER*4 ICOMP(0)
  55. INTEGER ICOMP(0)
  56. ENDSEGMENT
  57. SEGMENT MTRA2
  58. INTEGER MHAR(0)
  59. ENDSEGMENT
  60. logical descar
  61. *
  62. * Construisons d'abord l'ensemble des noms d'inconnues
  63. * primales, puis les duales associées (segment trav1)
  64. *
  65. C WRITE(IOIMP,*) 'Coucou exdiar'
  66. SEGACT MRIGID
  67. nbincd=0
  68. segini trav1
  69. NRIGEL=COERIG(/1)
  70. DO IRIG=1,NRIGEL
  71. DESCR=IRIGEL(3,IRIG)
  72. SEGACT DESCR
  73. do 45 iligrp=1,lisinc(/2)
  74. do 55 inc=1,nbincd
  75. if (compd(inc).eq.lisinc(iligrp)) goto 45
  76. 55 continue
  77. nbincd=nbincd+1
  78. segadj trav1
  79. compd(nbincd)=lisinc(iligrp)
  80. 45 continue
  81. SEGDES DESCR
  82. ENDDO
  83. do iincd=1,nbincd
  84. CALL PLACE(NOMDD,LNOMDD,idx,compd(iincd))
  85. IF (idx.EQ.0) THEN
  86. compdd(iincd)=compd(iincd)
  87. ELSE
  88. compdd(iincd)=NOMDU(idx)
  89. ENDIF
  90. ENDDO
  91. C WRITE(IOIMP,*) ' 1) compd'
  92. C WRITE (IOIMP,2019) (compd(I),I=1,compd(/2))
  93. C WRITE(IOIMP,*) ' 2) compdd'
  94. C WRITE (IOIMP,2019) (compdd(I),I=1,compdd(/2))
  95. *
  96. * Maintenant, on construit les CORES1 et CORES2 qui nous
  97. * permettent de savoir quels sont les termes diagonaux
  98. * dans les sous-matrices
  99. *
  100. SEGINI CORES1
  101. DO IRIG=1,NRIGEL
  102. DESCR=IRIGEL(3,IRIG)
  103. SEGACT DESCR
  104. * Tout d'abord le segment trav2 qui permet d'avoir
  105. * la correspondance sur les noms d'inconnues
  106. NLIGRP=LISINC(/2)
  107. NLIGRD=LISDUA(/2)
  108. segini trav2
  109. do 410 iligrp=1,nligrp
  110. do 420 incd=1,nbincd
  111. if (compd(incd).eq.lisinc(iligrp)) then
  112. cor1p(iligrp)=incd
  113. goto 410
  114. endif
  115. 420 continue
  116. 410 continue
  117. do 510 iligrd=1,nligrd
  118. do 520 incd=1,nbincd
  119. if (compdd(incd).eq.lisdua(iligrd)) then
  120. cor1d(iligrd)=incd
  121. goto 510
  122. endif
  123. 520 continue
  124. * Si on a pas trouvé dans la liste des duales, on regarde
  125. * dans la liste des primales.
  126. * Ceci est pratique en mécaflu où il arrive qu'on nomme les duales
  127. * pareil que les primales (ex. 'UX' 'UX')
  128. do 525 incd=1,nbincd
  129. if (compd(incd).eq.lisdua(iligrd)) then
  130. cor1d(iligrd)=incd
  131. goto 510
  132. endif
  133. 525 continue
  134. 510 continue
  135. * Voyons si le descripteur est carré, ce qui est un
  136. * cas courant qui simplifie notablement les choses
  137. descar=.false.
  138. if (nligrp.ne.nligrd) goto 600
  139. do iligr=1,nligrp
  140. if (cor1p(iligr).ne.cor1d(iligr)) goto 600
  141. if (noelep(iligr).ne.noeled(iligr)) goto 600
  142. enddo
  143. descar=.true.
  144. 600 continue
  145. * Si oui, tout est simple, sinon il faut rechercher
  146. * une correspondance primale-duale
  147. if (descar) then
  148. ipcor2(irig)=-1
  149. else
  150. nligp=nligrp
  151. segini cores2
  152. nligp=0
  153. do iligrp=1,nligrp
  154. do iligrd=1,nligrd
  155. if (cor1p(iligrp).eq.cor1d(iligrd).and.
  156. $ noelep(iligrp).eq.noeled(iligrd)) then
  157. nligp=nligp+1
  158. cordia(1,nligp)=iligrp
  159. cordia(2,nligp)=iligrd
  160. endif
  161. enddo
  162. enddo
  163. if (nligp.eq.0) then
  164. segsup cores2
  165. ipcor2(irig)=0
  166. else
  167. segadj cores2
  168. ipcor2(irig)=cores2
  169. endif
  170. endif
  171. itrav2(irig)=trav2
  172. C Write(ioimp,*) 'irig=',irig
  173. C Write(ioimp,*) 'descar=',descar
  174. C if (ipcor2.gt.0) then
  175. C write(ioimp,*) 'cordia(1, )'
  176. C WRITE (IOIMP,2020) (cordia(1,I),I=1,cordia(/2))
  177. C write(ioimp,*) 'cordia(2, )'
  178. C WRITE (IOIMP,2020) (cordia(2,I),I=1,cordia(/2))
  179. C else
  180. C WRITE(IOIMP,*) 'ipcor2=',ipcor2
  181. C endif
  182. SEGDES DESCR
  183. ENDDO
  184. *
  185. * Maintenant, on peut construire MTRA1, MTRA2
  186. * (repérage des couples noms d'inconnues primales,
  187. * harmoniques) ainsi que CORES3
  188. * et ICPR (repérage global->local des noeuds)
  189. * ce qui permettra dans une deuxième passe de remplir
  190. * le segment TMTRAV qui sera transformé en CHPOINT
  191. *
  192. SEGINI ICPR
  193. NNNOE=0
  194. SEGINI MTRA1,MTRA2
  195. NNIN=0
  196. DO 100 IRIG=1,NRIGEL
  197. DESCR=IRIGEL(3,IRIG)
  198. trav2=itrav2(irig)
  199. MELEME=IRIGEL(1,IRIG)
  200. cores2=ipcor2(irig)
  201. if (cores2.eq.0) goto 99
  202. SEGACT DESCR
  203. SEGACT MELEME
  204. if (cores2.eq.-1) then
  205. nligp=noelep(/1)
  206. else
  207. nligp=cordia(/2)
  208. endif
  209. segini cores3
  210. do 110 iligp=1,nligp
  211. if (cores2.eq.-1) then
  212. iligrp=iligp
  213. else
  214. iligrp=cordia(1,iligp)
  215. endif
  216. do ic=1,icomp(/1)
  217. if (icomp(ic).eq.cor1p(iligrp).and.
  218. $ mhar(ic).eq.irigel(5,irig)) then
  219. comtra(iligp)=ic
  220. goto 110
  221. endif
  222. enddo
  223. icomp(**)=cor1p(iligrp)
  224. mhar(**)=irigel(5,irig)
  225. nnin=nnin+1
  226. comtra(iligp)=nnin
  227. 110 continue
  228. C WRITE(IOIMP,*) ' 3) comtra'
  229. C WRITE (IOIMP,2020) (comtra(I),I=1,comtra(/1))
  230. ipcor3(irig)=cores3
  231. * Dans le cas cores2=-1
  232. * tous les noeuds du meleme sont référencés
  233. * sinon, on regarde ceux qui le sont
  234. if (cores2.eq.-1) then
  235. nno=num(/1)
  236. else
  237. nno=nligp
  238. endif
  239. do iel=1,num(/2)
  240. do ino=1,nno
  241. if (cores2.eq.-1) then
  242. ilno=ino
  243. else
  244. ilno=noelep(cordia(1,ino))
  245. endif
  246. ipt=num(ilno,iel)
  247. if (icpr(ipt).eq.0) then
  248. nnnoe=nnnoe+1
  249. icpr(ipt)=nnnoe
  250. endif
  251. enddo
  252. enddo
  253. 99 continue
  254. SEGDES MELEME
  255. SEGDES DESCR
  256. segsup trav2
  257. 100 CONTINUE
  258. C WRITE(IOIMP,*) 'nnin=',nnin
  259. C WRITE(IOIMP,*) 'nnin2=',icomp(/1)
  260. C WRITE(IOIMP,*) 'nnin3=',mhar(/1)
  261. C WRITE(IOIMP,*) ' 1) icomp'
  262. C WRITE (IOIMP,2020) (icomp(I),I=1,icomp(/1))
  263. C WRITE(IOIMP,*) ' 2) mhar'
  264. C WRITE (IOIMP,2020) (mhar(I),I=1,mhar(/1))
  265. C WRITE(IOIMP,*) 'nnnoe=',nnnoe
  266. C WRITE(IOIMP,*) ' 4) icpr'
  267. C WRITE (IOIMP,2020) (icpr(I),I=1,icpr(/1))
  268. *
  269. * Maintenant, on peut remplir MTRAV
  270. *
  271. SEGINI MTRAV
  272. DO ININ=1,NNIN
  273. INCO(ININ)=compd(ICOMP(ININ))
  274. NHAR(ININ)=MHAR(ININ)
  275. ENDDO
  276. SEGSUP MTRA1,MTRA2
  277. segsup trav1
  278. C xt=0.D0
  279. DO 1100 IRIG=1,NRIGEL
  280. MELEME=IRIGEL(1,IRIG)
  281. DESCR=IRIGEL(3,IRIG)
  282. XMATRI=IRIGEL(4,IRIG)
  283. cores2=ipcor2(irig)
  284. cores3=ipcor3(irig)
  285. if (cores2.eq.0) goto 1099
  286. SEGACT MELEME
  287. SEGACT DESCR
  288. SEGACT XMATRI
  289. nel=num(/2)
  290. if (cores2.eq.-1) then
  291. nligp=noelep(/1)
  292. else
  293. nligp=cordia(/2)
  294. endif
  295. do iel=1,nel
  296. do iligp=1,nligp
  297. if (cores2.eq.-1) then
  298. iligrp=iligp
  299. iligrd=iligrp
  300. else
  301. iligrp=cordia(1,iligp)
  302. iligrd=cordia(2,iligp)
  303. endif
  304. ININ=comtra(iligp)
  305. ino=noelep(iligrp)
  306. ipt=num(ino,iel)
  307. INNOE=icpr(ipt)
  308. C if (ipt.eq.2.and.lisinc(iligrp).eq.'UX ') then
  309. C write(ioimp,*) 'irig=',irig
  310. C write(ioimp,*) 'iel=',iel
  311. C write(ioimp,*) 'iligrp=',iligrp
  312. C write(ioimp,*) 'iligrd=',iligrd
  313. C xt=xt+RE(iligrd,iligrp,iel)
  314. C WRITE(IOIMP,*) 'RE=',RE(iligrd,iligrp,iel)
  315. C WRITE(IOIMP,*) 'xt=',xt
  316. C WRITE(IOIMP,*) 'inin=',inin
  317. C WRITE(IOIMP,*) 'innoe=',innoe
  318. C endif
  319. IGEO(INNOE)=ipt
  320. IBIN(ININ,INNOE)=1
  321. BB(ININ,INNOE)=BB(ININ,INNOE)+
  322. $ RE(iligrd,iligrp,iel)*COERIG(IRIG)
  323. enddo
  324. enddo
  325. if (cores2.ne.-1) segsup cores2
  326. segsup cores3
  327. SEGDES XMATRI
  328. SEGDES DESCR
  329. SEGDES MELEME
  330. 1099 continue
  331. 1100 CONTINUE
  332. SEGSUP ICPR
  333. SEGSUP CORES1
  334. SEGDES MRIGID
  335. *
  336. * Création du CHPOINT résultat
  337. *
  338. CALL CRECHP(MTRAV,MCHPOI)
  339. $
  340. SEGSUP MTRAV
  341. RETURN
  342. 2019 FORMAT (20(2X,A4) )
  343. 2020 FORMAT (20(2X,I4) )
  344. 2021 FORMAT (20(2X,L4) )
  345. 2022 FORMAT(10(1X,1PG12.5))
  346. *
  347. * End of EXDIAG
  348. *
  349. END
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales