Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes :

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

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