Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes :

  1. C MUCPR1 SOURCE PV 17/11/22 21:15:07 9630
  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. isyme=irigel(7,irige)
  226. if (ipt2.itypel.eq.22) then
  227. segact xmatri*mod
  228. else
  229. segact xmatri
  230. endif
  231. nod=noeled(/1)
  232. nop=noelep(/1)
  233. * if (nod.ne.nop) isyme=2
  234. do 55 iel=1,ipt2.num(/2)
  235. do 60 ip=1,nop
  236. if (icposp(ip).eq.0.or.icpr(ipt2.num(noelep(ip),iel)).eq.0)
  237. > then
  238. xbuffp(ip)=0.d0
  239. else
  240. * write (6,*) ' icposp icposd ',icposp(ip),icposd(ip)
  241. * write (6,*) ' ip xbuffp ',ip,xbuffp(ip)
  242. * write (6,*) ' icposp noelep ',icposp(ip),noelep(ip)
  243. * write (6,*) ' irige coerig ',irige,coerig(irige)
  244. * write (6,*) ' ipt2 ',ipt2.num(noelep(ip),iel)
  245. * write (6,*) ' icpr ',icpr(ipt2.num(noelep(ip),iel))
  246. * write (6,*) ' xvalp ',xvalp(icposp(ip),
  247. * > icpr(ipt2.num(noelep(ip),iel)))
  248. xbuffp(ip)=xvalp(icposp(ip),icpr(ipt2.num(noelep(ip),iel)))*
  249. > coerig(irige)
  250. endif
  251. 60 continue
  252.  
  253. * operation
  254. * on impose l'egalite des multiplicateurs de Lagrange dans mucpr2
  255. if (nop*nod.ne.0) then
  256. call mucpr2(nop,nod,re(1,1,iel),xbuffp,
  257. $ xbuffd,isyme)
  258. else
  259. do id=1,nod
  260. xbuffd(id)=0.D0
  261. enddo
  262. endif
  263. * transfert du buffer dual dans le vecteur dual etendu
  264. do 65 id=1,nod
  265. i1d=icposd(id)
  266. if (i1d.eq.0) then
  267. write (6,*) ' i1d nul '
  268. goto 65
  269. endif
  270. i2d=icpr(ipt2.num(noeled(id),iel))
  271. if (i2d.eq.0) goto 65
  272. xvald(i1d,i2d)=xvald(i1d,i2d)+xbuffd(id)
  273. 65 continue
  274. 55 continue
  275. * if (xmatri.ne.0) segdes xmatri
  276. segdes descr,ipt2,xmatri
  277. 50 continue
  278. segdes mrigid
  279. 51 continue
  280. *
  281. * transfert final dans le vecteur dual
  282. *
  283. do 80 isous=1,ipchp(/1)
  284. msoupo=ipchp(isous)
  285. mpoval=ipoval
  286. meleme=igeoc
  287. segact meleme
  288. do 82 i=1,nocomp(/2)
  289. do 83 j=1,nctod
  290. if (ifod(j).ne.noharm(i).and.ifodh(j).ne.noharm(i)) goto 83
  291. if (nocod(j).eq.nocomp(i)) goto 84
  292. 83 continue
  293. ifodp(i)=0
  294. goto 82
  295. call erreur(5)
  296. 84 continue
  297. if (ifodh(j).ne.0) noharm(i)=ifodh(j)-2**17
  298. * write (charm,fmt='(a4)') noharm(i)
  299. if (charm.eq.'NOHA') noharm(i)=nifour
  300. ifodp(i)=j
  301. 82 continue
  302. ieln=0
  303. ipt1=0
  304. do 85 iel=1,vpocha(/1)
  305. ieln=ieln+1
  306. if (ieln.ne.iel) num(1,ieln)=num(1,iel)
  307. ig=0
  308. do 86 i=1,vpocha(/2)
  309. if (ifodp(i).eq.0) then
  310. write (6,*) ' ifodp nul '
  311. goto 86
  312. endif
  313. xv=xvald(ifodp(i),icpr(num(1,iel)))
  314. vpocha(ieln,i)=xv
  315. * if (abs(xv).gt.1e-35) ig=1
  316. if (abs(xv).gt.xpetit) ig=1
  317. 86 continue
  318. if (ig.eq.0) then
  319. * il y a lieu de simplifier
  320. if (ipt1.eq.0) then
  321. segini,ipt1=meleme
  322. segdes meleme
  323. meleme=ipt1
  324. endif
  325. ieln=ieln-1
  326. endif
  327. 85 continue
  328. if (ieln.ne.vpocha(/1)) then
  329. nbnn=1
  330. nbelem=ieln
  331. nbref=0
  332. nbsous=0
  333. segadj meleme
  334. n=ieln
  335. nc=vpocha(/2)
  336. * write (6,*) ' chpoin simplifie ',vpocha(/1),ieln
  337. segadj mpoval
  338. endif
  339. if (meleme.ne.igeoc) call crech1(meleme,1)
  340. igeoc=meleme
  341. segdes mpoval,meleme
  342. segdes msoupo
  343. 80 continue
  344. segdes mchpoi,mrigid
  345. segsup itrav1,itrav2,itrav3,itrav4,itrav5,icpr
  346. end
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  

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