Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes :

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

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