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 155 : call erreur(5) 156 : 24 continue 157 : do 25 iel=1,mpova1.vpocha(/1) 158 : xvalp(j,icpr(ipt1.num(1,iel)))=mpova1.vpocha(iel,i) 159 : 25 continue 160 : 22 continue 161 : 20 continue 162 : * 163 : * creation de itrav5 164 : * 165 : nbpma=0 166 : nbdma=0 167 : do 27 irige=1,nrigel 168 : descr=irigel(3,irige) 169 : segact descr 170 : nbpma=max(nbpma,lisinc(/2)) 171 : nbdma=max(nbdma,lisdua(/2)) 172 : 27 continue 173 : segini itrav5 174 : if (nbdma.eq.0) goto 51 175 : 176 : * 177 : * travail effectif : boucle sur la matrice. 178 : * 179 : do 50 irige=1,nrigel 180 : * write (charm,fmt='(a4)') irigel(5,irige) 181 : descr=irigel(3,irige) 182 : segact descr 183 : * remplissage de icposp et icposd 184 : jsauv=0 185 : do 30 i=1,lisinc(/2) 186 : do 35 j=1,nctop 187 : if (ifop(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 35 188 : if (lisinc(i).eq.nocop(j)) goto 37 189 : 35 continue 190 : * inconnu du chp primal pas dans la raideur 191 : j=0 192 : 37 continue 193 : icposp(i)=j 193 : if (jsauv*j.ne.0.and.ifop(jsauv).ne.ifop(j)) call erreur(5) 194 : if (j.ne.0) jsauv=j 195 : 30 continue 196 : do 40 i=1,lisdua(/2) 197 : do 45 j=1,nctod 198 : if (ifod(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 45 199 : if (lisdua(i).eq.nocod(j)) goto 47 200 : 45 continue 201 : icposd(i)=0 202 : goto 40 203 : write (6,*) ' dans mucpr1 ',lisdua(/2),nctod 204 : write (6,*) irigel(5,irige),i,(lisdua(j),j=1,lisdua(/2)) 205 : write (6,*) (nocod(j),ifod(j),j=1,nctod) 205 : call erreur(5) 206 : 47 continue 207 : if (charm.eq.'NOHA'.and.jsauv.ne.0) then 208 : ifodh(j)=ifop(jsauv)+2**17 209 : endif 210 : icposd(i)=j 211 : 40 continue 212 : * transfert des valeurs dans le buffer primal 213 : ipt2=irigel(1,irige) 214 : segact ipt2 215 : * write (6,*) ' ipt2 ',ipt2.num(/1),ipt2.num(/2) 216 : xmatri=irigel(4,irige) 217 : if (ipt2.itypel.eq.22) then 218 : segact xmatri*mod 219 : else 220 : segact xmatri 221 : endif 222 : * en cas de blocage, on a souvent la meme matrice 223 : * on optimise ce cas 224 : * xmatri=0 225 : do 55 iel=1,ipt2.num(/2) 226 : * if (xmatri.ne.imattt(iel)) then 227 : * if (xmatri.ne.0) segdes xmatri 228 : * xmatri=imattt(iel) 229 : * if (ipt2.itypel.eq.22) then 230 : * segact xmatri*mod 231 : * else 232 : * segact xmatri 233 : * endif 234 : * endif 235 : ixmat=xmatri 236 : do 60 ip=1,lisinc(/2) 237 : if (icposp(ip).eq.0.or.icpr(ipt2.num(noelep(ip),iel)).eq.0) 238 : > then 239 : xbuffp(ip)=0.d0 240 : else 241 : * write (6,*) ' icposp icposd ',icposp(ip),icposd(ip) 242 : * write (6,*) ' ip xbuffp ',ip,xbuffp(ip) 243 : * write (6,*) ' icposp noelep ',icposp(ip),noelep(ip) 244 : * write (6,*) ' irige coerig ',irige,coerig(irige) 245 : * write (6,*) ' ipt2 ',ipt2.num(noelep(ip),iel) 246 : * write (6,*) ' icpr ',icpr(ipt2.num(noelep(ip),iel)) 247 : * write (6,*) ' xvalp ',xvalp(icposp(ip), 248 : * > icpr(ipt2.num(noelep(ip),iel))) 249 : xbuffp(ip)=xvalp(icposp(ip),icpr(ipt2.num(noelep(ip),iel)))* 250 : > coerig(irige) 251 : endif 252 : 60 continue 253 : 254 : nod=noeled(/1) 255 : nop=noelep(/1) 256 : niu= min (re(/1),1) 257 : nio=min(re(/2),1) 258 : nii=re(/3) 259 : * write(6,*)'niu nio nii iel',niu,nio,nii,iel,nop,nod 260 : * operation 261 : * on impose l'egalité des multiplicateurs de Lagrange dans mucpr2 262 : if (nio*nio.ne.0) then 262 : call mucpr2(nop,nod,re(niu,nio,iel),xbuffp, 263 : $ xbuffd,ipt2.itypel) 264 : else 265 : do id=1,lisdua(/2) 266 : xbuffd(id)=0.D0 267 : enddo 268 : endif 269 : * transfert du buffer dual dans le vecteur dual etendu 270 : do 65 id=1,lisdua(/2) 271 : i1d=icposd(id) 272 : if (i1d.eq.0) then 273 : write (6,*) ' i1d nul ' 274 : goto 65 275 : endif 276 : i2d=icpr(ipt2.num(noeled(id),iel)) 277 : if (i2d.eq.0) goto 65 278 : xvald(i1d,i2d)=xvald(i1d,i2d)+xbuffd(id) 279 : 65 continue 280 : 55 continue 281 : * if (xmatri.ne.0) segdes xmatri 282 : segdes descr,ipt2,xmatri 283 : 50 continue 284 : segdes mrigid 285 : 51 continue 286 : * 287 : * transfert final dans le vecteur dual 288 : * 289 : do 80 isous=1,ipchp(/1) 290 : msoupo=ipchp(isous) 291 : mpoval=ipoval 292 : meleme=igeoc 293 : segact meleme 294 : do 82 i=1,nocomp(/2) 295 : do 83 j=1,nctod 296 : if (ifod(j).ne.noharm(i).and.ifodh(j).ne.noharm(i)) goto 83 297 : if (nocod(j).eq.nocomp(i)) goto 84 298 : 83 continue 299 : ifodp(i)=0 300 : goto 82 300 : call erreur(5) 301 : 84 continue 302 : if (ifodh(j).ne.0) noharm(i)=ifodh(j)-2**17 303 : * write (charm,fmt='(a4)') noharm(i) 304 : if (charm.eq.'NOHA') noharm(i)=nifour 305 : ifodp(i)=j 306 : 82 continue 307 : ieln=0 308 : ipt1=0 309 : do 85 iel=1,vpocha(/1) 310 : ieln=ieln+1 311 : if (ieln.ne.iel) num(1,ieln)=num(1,iel) 312 : ig=0 313 : do 86 i=1,vpocha(/2) 314 : if (ifodp(i).eq.0) then 315 : write (6,*) ' ifodp nul ' 316 : goto 86 317 : endif 318 : xv=xvald(ifodp(i),icpr(num(1,iel))) 319 : vpocha(ieln,i)=xv 320 : * if (abs(xv).gt.1e-35) ig=1 321 : if (abs(xv).gt.xpetit) ig=1 322 : 86 continue 323 : if (ig.eq.0) then 324 : * il y a lieu de simplifier 325 : if (ipt1.eq.0) then 326 : segini,ipt1=meleme 327 : segdes meleme 328 : meleme=ipt1 329 : endif 330 : ieln=ieln-1 331 : endif 332 : 85 continue 333 : if (ieln.ne.vpocha(/1)) then 334 : nbnn=1 335 : nbelem=ieln 336 : nbref=0 337 : nbsous=0 338 : segadj meleme 339 : n=ieln 340 : nc=vpocha(/2) 341 : * write (6,*) ' chpoin simplifie ',vpocha(/1),ieln 342 : segadj mpoval 343 : endif 343 : if (meleme.ne.igeoc) call crech1(meleme,1) 344 : igeoc=meleme 345 : segdes mpoval,meleme 346 : segdes msoupo 347 : 80 continue 348 : segdes mchpoi,mrigid 349 : segsup itrav1,itrav2,itrav3,itrav4,itrav5,icpr 350 : end 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 :
© Cast3M 2003 - All rights reserved.
Disclaimer