Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes :

  1. C MUCPR1 SOURCE PV 13/04/18 21:15:01 7771
  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. *
  46. segact mcoord
  47. nbpts=xcoor(/1)/(idim+1)
  48. *
  49. * creation du champ resultat
  50. *
  51. segini icpr
  52. ncpr=0
  53. segact mrigid
  54. mchpo2=imgeo2
  55. segact mchpo2
  56. nrigel=irigel(/2)
  57. segini,mchpoi=mchpo2
  58. segdes mchpo2
  59. nbz=ipchp(/1)
  60. do 1 isous=1,nbz
  61. msoup2=ipchp(isous)
  62. segini,msoupo=msoup2
  63. segdes msoup2
  64. ipchp(isous)=msoupo
  65. meleme=igeoc
  66. segact meleme
  67. n=num(/2)
  68. do 2 in=1,num(/1)
  69. do 2 im=1,num(/2)
  70. ip=num(in,im)
  71. if (icpr(ip).ne.0) goto 2
  72. ncpr=ncpr+1
  73. icpr(ip)=ncpr
  74. 2 continue
  75. nc=noharm(/1)
  76. segini mpoval
  77. ipoval=mpoval
  78. 1 continue
  79. *
  80. * constitution de la liste des composantes du chant primal
  81. *
  82. nctop=0
  83. segact mchpo1
  84. do 10 isoupo=1,mchpo1.ipchp(/1)
  85. msoup1=mchpo1.ipchp(isoupo)
  86. segact msoup1
  87. mpova1=msoup1.ipoval
  88. segact mpova1
  89. ipt1=msoup1.igeoc
  90. segact ipt1
  91. do 3 in=1,ipt1.num(/1)
  92. do 3 im=1,ipt1.num(/2)
  93. ip=ipt1.num(in,im)
  94. if (icpr(ip).ne.0) goto 3
  95. ncpr=ncpr+1
  96. icpr(ip)=ncpr
  97. 3 continue
  98.  
  99. nctop=nctop+msoup1.nocomp(/2)
  100. 10 continue
  101. segini itrav1
  102. nctop=0
  103. do 11 isous=1,mchpo1.ipchp(/1)
  104. msoup1=mchpo1.ipchp(isous)
  105. do 12 i=1,msoup1.nocomp(/2)
  106. do 13 j=1,nctop
  107. if (ifop(j).ne.msoup1.noharm(i)) goto 13
  108. if (nocop(j).eq.msoup1.nocomp(i)) goto 12
  109. 13 continue
  110. nctop=nctop+1
  111. nocop(nctop)=msoup1.nocomp(i)
  112. ifop(nctop)=msoup1.noharm(i)
  113. 12 continue
  114. 11 continue
  115. * write (6,*) ' mucpr1 nbpts,ncpr,nctop ',
  116. * > nbpts,ncpr,nctop
  117. *
  118. * constitution de la liste des composantes du chant dual
  119. *
  120. nctod=0
  121. do 15 isoupo=1,ipchp(/1)
  122. msoupo=ipchp(isoupo)
  123. mpoval=ipoval
  124. segact mpoval*mod
  125. nctod=nctod+nocomp(/2)
  126. 15 continue
  127. * write (6,*) ' dans mucpr1 nctod-1 ',nctod,ipchp(/1)
  128. segini itrav2
  129. nctod=0
  130. do 16 isous=1,ipchp(/1)
  131. msoupo=ipchp(isous)
  132. do 17 i=1,nocomp(/2)
  133. do 18 j=1,nctod
  134. if (ifod(j).ne.noharm(i)) goto 18
  135. if (nocod(j).eq.nocomp(i)) goto 17
  136. 18 continue
  137. nctod=nctod+1
  138. nocod(nctod)=nocomp(i)
  139. ifod(nctod)=noharm(i)
  140. 17 continue
  141. 16 continue
  142. * write (6,*) ' dans mucpr1 nctod-2 ',nctod
  143. *
  144. * expansion du champ primal
  145. *
  146. segini itrav3,itrav4
  147. do 20 isous=1,mchpo1.ipchp(/1)
  148. msoup1=mchpo1.ipchp(isous)
  149. mpova1=msoup1.ipoval
  150. ipt1=msoup1.igeoc
  151. do 22 i=1,msoup1.nocomp(/2)
  152. do 23 j=1,nctop
  153. if (ifop(j).ne.msoup1.noharm(i)) goto 23
  154. if (nocop(j).eq.msoup1.nocomp(i)) goto 24
  155. 23 continue
  156. call erreur(5)
  157. 24 continue
  158. do 25 iel=1,mpova1.vpocha(/1)
  159. xvalp(j,icpr(ipt1.num(1,iel)))=mpova1.vpocha(iel,i)
  160. 25 continue
  161. 22 continue
  162. 20 continue
  163. *
  164. * creation de itrav5
  165. *
  166. nbpma=0
  167. nbdma=0
  168. do 27 irige=1,nrigel
  169. descr=irigel(3,irige)
  170. segact descr
  171. nbpma=max(nbpma,lisinc(/2))
  172. nbdma=max(nbdma,lisdua(/2))
  173. 27 continue
  174. segini itrav5
  175. if (nbdma.eq.0) goto 51
  176.  
  177. *
  178. * travail effectif : boucle sur la matrice.
  179. *
  180. do 50 irige=1,nrigel
  181. * write (charm,fmt='(a4)') irigel(5,irige)
  182. descr=irigel(3,irige)
  183. segact descr
  184. * remplissage de icposp et icposd
  185. jsauv=0
  186. do 30 i=1,lisinc(/2)
  187. do 35 j=1,nctop
  188. if (ifop(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 35
  189. if (lisinc(i).eq.nocop(j)) goto 37
  190. 35 continue
  191. * inconnu du chp primal pas dans la raideur
  192. j=0
  193. 37 continue
  194. icposp(i)=j
  195. if (jsauv*j.ne.0.and.ifop(jsauv).ne.ifop(j)) call erreur(5)
  196. if (j.ne.0) jsauv=j
  197. 30 continue
  198. do 40 i=1,lisdua(/2)
  199. do 45 j=1,nctod
  200. if (ifod(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 45
  201. if (lisdua(i).eq.nocod(j)) goto 47
  202. 45 continue
  203. icposd(i)=0
  204. goto 40
  205. write (6,*) ' dans mucpr1 ',lisdua(/2),nctod
  206. write (6,*) irigel(5,irige),i,(lisdua(j),j=1,lisdua(/2))
  207. write (6,*) (nocod(j),ifod(j),j=1,nctod)
  208. call erreur(5)
  209. 47 continue
  210. if (charm.eq.'NOHA'.and.jsauv.ne.0) then
  211. ifodh(j)=ifop(jsauv)+2**17
  212. endif
  213. icposd(i)=j
  214. 40 continue
  215. * transfert des valeurs dans le buffer primal
  216. ipt2=irigel(1,irige)
  217. segact ipt2
  218. * write (6,*) ' ipt2 ',ipt2.num(/1),ipt2.num(/2)
  219. xmatri=irigel(4,irige)
  220. if (ipt2.itypel.eq.22) then
  221. segact xmatri*mod
  222. else
  223. segact xmatri
  224. endif
  225. * en cas de blocage, on a souvent la meme matrice
  226. * on optimise ce cas
  227. * xmatri=0
  228. do 55 iel=1,ipt2.num(/2)
  229. * if (xmatri.ne.imattt(iel)) then
  230. * if (xmatri.ne.0) segdes xmatri
  231. * xmatri=imattt(iel)
  232. * if (ipt2.itypel.eq.22) then
  233. * segact xmatri*mod
  234. * else
  235. * segact xmatri
  236. * endif
  237. * endif
  238. ixmat=xmatri
  239. do 60 ip=1,lisinc(/2)
  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. nod=noeled(/1)
  258. nop=noelep(/1)
  259. niu= min (re(/1),1)
  260. nio=min(re(/2),1)
  261. nii=re(/3)
  262. * write(6,*)'niu nio nii iel',niu,nio,nii,iel,nop,nod
  263. * operation
  264. * on impose l'egalité des multiplicateurs de Lagrange dans mucpr2
  265. if (nio*nio.ne.0) then
  266. call mucpr2(nop,nod,re(niu,nio,iel),xbuffp,
  267. $ xbuffd,ipt2.itypel)
  268. else
  269. do id=1,lisdua(/2)
  270. xbuffd(id)=0.D0
  271. enddo
  272. endif
  273. * transfert du buffer dual dans le vecteur dual etendu
  274. do 65 id=1,lisdua(/2)
  275. i1d=icposd(id)
  276. if (i1d.eq.0) then
  277. write (6,*) ' i1d nul '
  278. goto 65
  279. endif
  280. i2d=icpr(ipt2.num(noeled(id),iel))
  281. if (i2d.eq.0) goto 65
  282. xvald(i1d,i2d)=xvald(i1d,i2d)+xbuffd(id)
  283. 65 continue
  284. 55 continue
  285. * if (xmatri.ne.0) segdes xmatri
  286. segdes descr,ipt2,xmatri
  287. 50 continue
  288. segdes mrigid
  289. 51 continue
  290. *
  291. * transfert final dans le vecteur dual
  292. *
  293. do 80 isous=1,ipchp(/1)
  294. msoupo=ipchp(isous)
  295. mpoval=ipoval
  296. meleme=igeoc
  297. segact meleme
  298. do 82 i=1,nocomp(/2)
  299. do 83 j=1,nctod
  300. if (ifod(j).ne.noharm(i).and.ifodh(j).ne.noharm(i)) goto 83
  301. if (nocod(j).eq.nocomp(i)) goto 84
  302. 83 continue
  303. ifodp(i)=0
  304. goto 82
  305. call erreur(5)
  306. 84 continue
  307. if (ifodh(j).ne.0) noharm(i)=ifodh(j)-2**17
  308. * write (charm,fmt='(a4)') noharm(i)
  309. if (charm.eq.'NOHA') noharm(i)=nifour
  310. ifodp(i)=j
  311. 82 continue
  312. ieln=0
  313. ipt1=0
  314. do 85 iel=1,vpocha(/1)
  315. ieln=ieln+1
  316. if (ieln.ne.iel) num(1,ieln)=num(1,iel)
  317. ig=0
  318. do 86 i=1,vpocha(/2)
  319. if (ifodp(i).eq.0) then
  320. write (6,*) ' ifodp nul '
  321. goto 86
  322. endif
  323. xv=xvald(ifodp(i),icpr(num(1,iel)))
  324. vpocha(ieln,i)=xv
  325. * if (abs(xv).gt.1e-35) ig=1
  326. if (abs(xv).gt.xpetit) ig=1
  327. 86 continue
  328. if (ig.eq.0) then
  329. * il y a lieu de simplifier
  330. if (ipt1.eq.0) then
  331. segini,ipt1=meleme
  332. segdes meleme
  333. meleme=ipt1
  334. endif
  335. ieln=ieln-1
  336. endif
  337. 85 continue
  338. if (ieln.ne.vpocha(/1)) then
  339. nbnn=1
  340. nbelem=ieln
  341. nbref=0
  342. nbsous=0
  343. segadj meleme
  344. n=ieln
  345. nc=vpocha(/2)
  346. * write (6,*) ' chpoin simplifie ',vpocha(/1),ieln
  347. segadj mpoval
  348. endif
  349. if (meleme.ne.igeoc) call crech1(meleme,1)
  350. igeoc=meleme
  351. segdes mpoval,meleme
  352. segdes msoupo
  353. 80 continue
  354. segdes mchpoi,mrigid
  355. segsup itrav1,itrav2,itrav3,itrav4,itrav5,icpr
  356. end
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  

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