Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes :

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

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