Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes :

  1. C MUCPR1 SOURCE PV 20/05/23 21:15:00 10616
  2. subroutine mucpr1(mchpo1,mrigid,mchpoi)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. * multiplication rapide d'une matrice par un champoint dans le cas
  6. * ou on connait deja la structure du chanpoint resultant, qui a par
  7. * exemple ete obtenu par un appel precedent a mucpri
  8. *
  9. * entree : mrigid, mchpo1
  10. * sortie : mchpoi
  11. *
  12. -INC CCREEL
  13.  
  14. -INC PPARAM
  15. -INC CCOPTIO
  16. -INC SMCHPOI
  17. -INC SMELEME
  18. -INC SMCOORD
  19. -INC SMRIGID
  20. segment itrav1
  21. * liste des inconnues primales dans le champ etendu
  22. character*4 nocop(nctop)
  23. integer ifop(nctop)
  24. endsegment
  25. segment itrav2
  26. * liste des inconnues duales dans le champ etendu
  27. character*4 nocod(nctod)
  28. integer ifod(nctod),ifodh(nctod),ifodp(nctod)
  29. endsegment
  30. segment itrav3
  31. * champ primal etendu
  32. real*8 xvalp(nctop,ncpr )
  33. endsegment
  34. segment itrav4
  35. * champ dual etendu
  36. real*8 xvald(nctod,ncpr )
  37. endsegment
  38. segment itrav5
  39. * position inconnues primales et duales de la raideur dans les champs etendus
  40. integer icposp(nbpma),icposd(nbdma)
  41. * buffers primal et dual
  42. real*8 xbuffp(nbpma),xbuffd(nbdma)
  43. character*4 charm
  44. endsegment
  45. segment icpr(nbpts)
  46.  
  47. character*4 cnoha
  48. integer*4 inoha
  49. equivalence(cnoha,inoha)
  50. data cnoha/'NOHA'/
  51.  
  52. *
  53. *
  54. * creation du champ resultat
  55. *
  56. segini icpr
  57. ncpr=0
  58. segact mrigid
  59. mchpo2=imgeo2
  60. segact mchpo2
  61. nrigel=irigel(/2)
  62. segini,mchpoi=mchpo2
  63. nbz=ipchp(/1)
  64. do 1 isous=1,nbz
  65. msoup2=ipchp(isous)
  66. segini,msoupo=msoup2
  67. ipchp(isous)=msoupo
  68. meleme=igeoc
  69. segact meleme
  70. n=num(/2)
  71. do im=1,num(/2)
  72. do in=1,num(/1)
  73. ip=num(in,im)
  74. if (icpr(ip).eq.0) then
  75. ncpr=ncpr+1
  76. icpr(ip)=ncpr
  77. endif
  78. enddo
  79. enddo
  80. nc=noharm(/1)
  81. segini mpoval
  82. ipoval=mpoval
  83. 1 continue
  84. *
  85. * constitution de la liste des composantes du chant primal
  86. *
  87. nctop=0
  88. segact mchpo1
  89. do 10 isoupo=1,mchpo1.ipchp(/1)
  90. msoup1=mchpo1.ipchp(isoupo)
  91. segact msoup1
  92. mpova1=msoup1.ipoval
  93. segact mpova1
  94. ipt1=msoup1.igeoc
  95. segact ipt1
  96. do im=1,ipt1.num(/2)
  97. do in=1,ipt1.num(/1)
  98. ip=ipt1.num(in,im)
  99. if (icpr(ip).eq.0) then
  100. ncpr=ncpr+1
  101. icpr(ip)=ncpr
  102. endif
  103. enddo
  104. enddo
  105.  
  106. nctop=nctop+msoup1.nocomp(/2)
  107. 10 continue
  108. segini itrav1
  109. nctop=0
  110. do 11 isous=1,mchpo1.ipchp(/1)
  111. msoup1=mchpo1.ipchp(isous)
  112. do 12 i=1,msoup1.nocomp(/2)
  113. do 13 j=1,nctop
  114. if (ifop(j).ne.msoup1.noharm(i)) goto 13
  115. if (nocop(j).eq.msoup1.nocomp(i)) goto 12
  116. 13 continue
  117. nctop=nctop+1
  118. nocop(nctop)=msoup1.nocomp(i)
  119. ifop(nctop)=msoup1.noharm(i)
  120. 12 continue
  121. 11 continue
  122. * write (6,*) ' mucpr1 nbpts,ncpr,nctop ',
  123. * > nbpts,ncpr,nctop
  124. *
  125. * constitution de la liste des composantes du chant dual
  126. *
  127. nctod=0
  128. do 15 isoupo=1,ipchp(/1)
  129. msoupo=ipchp(isoupo)
  130. mpoval=ipoval
  131. segact mpoval*mod
  132. nctod=nctod+nocomp(/2)
  133. 15 continue
  134. * write (6,*) ' dans mucpr1 nctod-1 ',nctod,ipchp(/1)
  135. segini itrav2
  136. nctod=0
  137. do 16 isous=1,ipchp(/1)
  138. msoupo=ipchp(isous)
  139. do 17 i=1,nocomp(/2)
  140. do 18 j=1,nctod
  141. if (ifod(j).ne.noharm(i)) goto 18
  142. if (nocod(j).eq.nocomp(i)) goto 17
  143. 18 continue
  144. nctod=nctod+1
  145. nocod(nctod)=nocomp(i)
  146. ifod(nctod)=noharm(i)
  147. 17 continue
  148. 16 continue
  149. * write (6,*) ' dans mucpr1 nctod-2 ',nctod
  150. *
  151. * expansion du champ primal
  152. *
  153. segini itrav3,itrav4
  154. do 20 isous=1,mchpo1.ipchp(/1)
  155. msoup1=mchpo1.ipchp(isous)
  156. mpova1=msoup1.ipoval
  157. ipt1=msoup1.igeoc
  158. do 22 i=1,msoup1.nocomp(/2)
  159. do 23 j=1,nctop
  160. if (ifop(j).ne.msoup1.noharm(i)) goto 23
  161. if (nocop(j).eq.msoup1.nocomp(i)) goto 24
  162. 23 continue
  163. call erreur(5)
  164. 24 continue
  165. do 25 iel=1,mpova1.vpocha(/1)
  166. xvalp(j,icpr(ipt1.num(1,iel)))=mpova1.vpocha(iel,i)
  167. 25 continue
  168. 22 continue
  169. 20 continue
  170. *
  171. * creation de itrav5
  172. *
  173. nbpma=0
  174. nbdma=0
  175. do 27 irige=1,nrigel
  176. descr=irigel(3,irige)
  177. segact descr
  178. nbpma=max(nbpma,lisinc(/2))
  179. nbdma=max(nbdma,lisdua(/2))
  180. 27 continue
  181. segini itrav5
  182. if (nbdma.eq.0) goto 51
  183.  
  184. *
  185. * travail effectif : boucle sur la matrice.
  186. *
  187. do 50 irige=1,nrigel
  188. * write (charm,fmt='(a4)') irigel(5,irige)
  189. descr=irigel(3,irige)
  190. segact descr
  191. * remplissage de icposp et icposd
  192. jsauv=0
  193. do 30 i=1,lisinc(/2)
  194. do 35 j=1,nctop
  195. if (ifop(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 35
  196. if (lisinc(i).eq.nocop(j)) goto 37
  197. 35 continue
  198. * inconnu du chp primal pas dans la raideur
  199. j=0
  200. 37 continue
  201. icposp(i)=j
  202. if (jsauv*j.ne.0) then
  203. if (ifop(jsauv).ne.ifop(j)) call erreur(5)
  204. endif
  205. if (j.ne.0) jsauv=j
  206. 30 continue
  207. do 40 i=1,lisdua(/2)
  208. do 45 j=1,nctod
  209. if (ifod(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 45
  210. if (lisdua(i).eq.nocod(j)) goto 47
  211. 45 continue
  212. icposd(i)=0
  213. goto 40
  214. write (6,*) ' dans mucpr1 ',lisdua(/2),nctod
  215. write (6,*) irigel(5,irige),i,(lisdua(j),j=1,lisdua(/2))
  216. write (6,*) (nocod(j),ifod(j),j=1,nctod)
  217. call erreur(5)
  218. 47 continue
  219. if (charm.eq.'NOHA'.and.jsauv.ne.0) then
  220. ifodh(j)=ifop(jsauv)+2**17
  221. endif
  222. icposd(i)=j
  223. 40 continue
  224. * transfert des valeurs dans le buffer primal
  225. ipt2=irigel(1,irige)
  226. segact ipt2
  227. * write (6,*) ' ipt2 ',ipt2.num(/1),ipt2.num(/2)
  228. xmatri=irigel(4,irige)
  229. isyme=irigel(7,irige)
  230. if (ipt2.itypel.eq.22) then
  231. segact xmatri*mod
  232. else
  233. segact xmatri
  234. endif
  235. nod=noeled(/1)
  236. nop=noelep(/1)
  237. * if (nod.ne.nop) isyme=2
  238. do 55 iel=1,ipt2.num(/2)
  239. do 60 ip=1,nop
  240. if (icposp(ip).eq.0.or.icpr(ipt2.num(noelep(ip),iel)).eq.0)
  241. > then
  242. xbuffp(ip)=0.d0
  243. else
  244. * write (6,*) ' icposp icposd ',icposp(ip),icposd(ip)
  245. * write (6,*) ' ip xbuffp ',ip,xbuffp(ip)
  246. * write (6,*) ' icposp noelep ',icposp(ip),noelep(ip)
  247. * write (6,*) ' irige coerig ',irige,coerig(irige)
  248. * write (6,*) ' ipt2 ',ipt2.num(noelep(ip),iel)
  249. * write (6,*) ' icpr ',icpr(ipt2.num(noelep(ip),iel))
  250. * write (6,*) ' xvalp ',xvalp(icposp(ip),
  251. * > icpr(ipt2.num(noelep(ip),iel)))
  252. xbuffp(ip)=xvalp(icposp(ip),icpr(ipt2.num(noelep(ip),iel)))*
  253. > coerig(irige)
  254. endif
  255. 60 continue
  256.  
  257. * operation
  258. * on impose l'egalite des multiplicateurs de Lagrange dans mucpr2
  259. if (nop*nod.ne.0) then
  260. call mucpr2(nop,nod,re(1,1,iel),xbuffp,
  261. $ xbuffd,isyme)
  262. else
  263. do id=1,nod
  264. xbuffd(id)=0.D0
  265. enddo
  266. endif
  267. * transfert du buffer dual dans le vecteur dual etendu
  268. do 65 id=1,nod
  269. i1d=icposd(id)
  270. if (i1d.eq.0) then
  271. write (6,*) ' i1d nul '
  272. goto 65
  273. endif
  274. i2d=icpr(ipt2.num(noeled(id),iel))
  275. if (i2d.eq.0) goto 65
  276. xvald(i1d,i2d)=xvald(i1d,i2d)+xbuffd(id)
  277. 65 continue
  278. 55 continue
  279. * if (xmatri.ne.0) segdes xmatri
  280. segdes descr,xmatri
  281. 50 continue
  282. segdes mrigid
  283. 51 continue
  284. *
  285. * transfert final dans le vecteur dual
  286. *
  287. do 80 isous=1,ipchp(/1)
  288. msoupo=ipchp(isous)
  289. mpoval=ipoval
  290. meleme=igeoc
  291. segact meleme
  292. do 82 i=1,nocomp(/2)
  293. do 83 j=1,nctod
  294. if (ifod(j).ne.noharm(i).and.ifodh(j).ne.noharm(i)) goto 83
  295. if (nocod(j).eq.nocomp(i)) goto 84
  296. 83 continue
  297. ifodp(i)=0
  298. goto 82
  299. call erreur(5)
  300. 84 continue
  301. if (ifodh(j).ne.0) noharm(i)=ifodh(j)-2**17
  302. * write (charm,fmt='(a4)') noharm(i)
  303. if (charm.eq.'NOHA') noharm(i)=nifour
  304. ifodp(i)=j
  305. 82 continue
  306. ieln=0
  307. ipt1=0
  308. do 85 iel=1,vpocha(/1)
  309. ieln=ieln+1
  310. if (ieln.ne.iel) num(1,ieln)=num(1,iel)
  311. ig=0
  312. do 86 i=1,vpocha(/2)
  313. if (ifodp(i).eq.0) then
  314. write (6,*) ' ifodp nul '
  315. goto 86
  316. endif
  317. xv=xvald(ifodp(i),icpr(num(1,iel)))
  318. vpocha(ieln,i)=xv
  319. * if (abs(xv).gt.1e-35) ig=1
  320. if (abs(xv).gt.xpetit) ig=1
  321. 86 continue
  322. if (ig.eq.0) then
  323. * il y a lieu de simplifier
  324. if (ipt1.eq.0) then
  325. segini,ipt1=meleme
  326. meleme=ipt1
  327. endif
  328. ieln=ieln-1
  329. endif
  330. 85 continue
  331. if (ieln.ne.vpocha(/1)) then
  332. nbnn =1
  333. nbelem=ieln
  334. nbref =0
  335. nbsous=0
  336. segadj,meleme
  337. n=ieln
  338. nc=vpocha(/2)
  339. * write (6,*) ' chpoin simplifie ',vpocha(/1),ieln
  340. segadj mpoval
  341. endif
  342. if (meleme.ne.igeoc) call crech1(meleme,1)
  343. igeoc=meleme
  344. segact,mpoval*nomod,meleme*nomod,msoupo*nomod
  345. 80 continue
  346. segact,mchpoi*nomod
  347. segdes,mrigid
  348. segsup itrav1,itrav2,itrav3,itrav4,itrav5,icpr
  349. end
  350.  
  351.  
  352.  
  353.  
  354.  

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