1 : C SEPARM SOURCE FANDEUR 13/02/07 21:15:14 7694 2 : subroutine separm(mrigid,ri1,ri2,nounil,lagdua) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : -INC SMRIGID 6 : -INC SMCOORD 7 : -INC SMELEME 8 : -INC CCOPTIO 9 : * extrait de mrigid dans ri1 la partie des relations pouvant etre 10 : * traitees comme des dependances (elimination des inconnues) 11 : * il faut que la premiere inconnue d'une dependance n'apparaisse pas 12 : * dans les autres dependances 13 : * juillet 2003 passage aux simples multiplicateurs de Lagrange, mais 14 : * separm dualise ceux qu'il garde 15 : *-INC CCHAMP 16 : segment trav1 17 : character*4 compp(lcomp1) 18 : C* character*4 compp(nbincp) 19 : endsegment 20 : segment trav2 21 : integer ielim(nbno,nbincp),iautr(nbno,nbincp),icomb(nbno,nbincp) 22 : endsegment 23 : segment trav3 24 : integer posinc(nligrp) 25 : endsegment 26 : character*4 co1 27 : integer ri1s,ri2s,ri1f 28 : 29 : * write (6,*) ' dans separm nouinl ',nounil 30 : nbdep=0 31 : nbpiv=0 32 : nbinc=0 33 : nbpvt=0 34 : nbno=xcoor(/1)/(idim+1) 35 : segact mrigid 36 : * remplissage du segment trav1 (liste des noms de composantes) 37 : lcomp1 = 50 38 : segini trav1 39 : nbincp=0 40 : do 10 irig=1,irigel(/2) 41 : if (irigel(6,irig).ne.0.and.irigel(6,irig).ne.2 42 : & .and.nounil.eq.0) goto 10 43 : meleme=irigel(1,irig) 44 : segact meleme 45 : if (itypel.ne.22) goto 20 46 : descr=irigel(3,irig) 47 : segact descr 48 : do 40 nligrp=2,lisinc(/2) 49 : C print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp) 50 : do 50 inc=1,nbincp 51 : if (compp(inc).eq.lisinc(nligrp)) goto 40 52 : 50 continue 53 : nbincp=nbincp+1 54 : if (nbincp.GT.lcomp1) then 55 : lcomp1 = lcomp1 + 50 56 : segadj trav1 57 : endif 58 : compp(nbincp)=lisinc(nligrp) 59 : 40 continue 60 : segdes,descr 61 : 20 continue 62 : segdes,meleme 63 : 10 continue 64 : c* lcomp1 = nbincp 65 : c* segadj trav1 66 : * on construit ri1 et ri2 en notant dans le segment trav2 ce qu'on fait 67 : segini trav2 68 : * on interdit d'eliminer les inconnues apparaissant dans le frottement 69 : do 60 irig=1,irigel(/2) 70 : if (irigel(6,irig).eq.0.or.nounil.eq.1) goto 60 71 : meleme=irigel(1,irig) 72 : segact meleme 73 : if (itypel.ne.22) goto 60 74 : descr=irigel(3,irig) 75 : segact descr 76 : do 70 nligrp=2,lisinc(/2) 77 : do 80 inc=1,nbincp 78 : C print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp) 79 : C & ,' compp(',inc,')=',compp(inc) 80 : if (compp(inc).eq.lisinc(nligrp)) goto 90 81 : 80 continue 82 : c* write(ioimp,*) '- ',lisinc(nligrp),' non trouve dans trav1' 83 : goto 70 84 : 90 continue 85 : ip=noelep(nligrp) 86 : do 100 j=1,num(/2) 87 : C IF(num(ip,j).GT.ielim(/1).OR.inc.GT.ielim(/2).OR. 88 : C * num(ip,j).GT.iautr(/1).OR.inc.GT.iautr(/2))THEN 89 : C print *,'BUG DANS SEPARM : ', 90 : C & num(ip,j),'>',ielim(/1),iautr(/1), 91 : C & inc,'>',ielim(/2),iautr(/2) 92 : C ENDIF 93 : ielim(num(ip,j),inc)=1 94 : iautr(num(ip,j),inc)=1 95 : 100 continue 96 : 70 continue 97 : segdes descr,meleme 98 : 60 continue 99 : do 700 irig=1,irigel(/2) 100 : if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 690 101 : meleme=irigel(1,irig) 102 : descr=irigel(3,irig) 103 : xmatri=irigel(4,irig) 104 : segact meleme 105 : if (itypel.ne.22) goto 690 106 : segact descr 107 : * remplir posinc pour accelerer les tests sur les composantes 108 : nligrp=lisinc(/2) 109 : segini trav3 110 : do 730 inc=2,nligrp 111 : do 720 incp=1,nbincp 112 : if (lisinc(inc).eq.compp(incp)) then 113 : posinc(inc)=incp 114 : goto 730 115 : endif 116 : 720 continue 116 : call erreur(5) 117 : 730 continue 118 : do 750 ince=2,lisinc(/2) 119 : incp=posinc(ince) 120 : do 740 j=1,num(/2) 121 : ip=num(noelep(ince),j) 122 : icomb(ip,incp)=icomb(ip,incp)+1 123 : 740 continue 124 : 750 continue 125 : segsup trav3 126 : 690 continue 127 : 700 continue 128 : ipass=1 129 : 2000 continue 130 : segini,ri1=mrigid 131 : segini,ri2=mrigid 132 : do 200 irig=1,irigel(/2) 133 : if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 190 134 : meleme=irigel(1,irig) 135 : descr=irigel(3,irig) 136 : Xmatri=irigel(4,irig) 137 : segact meleme 138 : if (itypel.ne.22) goto 190 139 : segact descr,Xmatri 140 : * remplir posinc pour accelerer les tests sur les composantes 141 : nligrp=lisinc(/2) 142 : segini trav3 143 : do 230 inc=2,nligrp 144 : do 220 incp=1,nbincp 145 : if (lisinc(inc).eq.compp(incp)) then 146 : posinc(inc)=incp 147 : goto 230 148 : endif 149 : 220 continue 149 : call erreur(5) 150 : 230 continue 151 : nbdepe=0 152 : * pour arriver à eliminer, il est interessant de faire pivoter les matrices 153 : * on va le faire en se basant sur le milieu du paquet 154 : iel=num(/2)/2+1 155 : incf=0 156 : do 1000 ince=2,lisinc(/2) 157 : ipp=num(noelep(ince),iel) 158 : incpp=posinc(ince) 159 : * cette inconnue apparait elle dejà dans une autre relation 160 : if (iautr(ipp,incpp).eq.1) goto 1000 161 : * cette matrice touche t elle une autre relation 162 : do 1240 inc=2,lisinc(/2) 163 : if (inc.eq.ince) goto 1240 164 : incp=posinc(inc) 165 : ip=num(noelep(inc), iel) 166 : if (ielim(ip,incp).eq.1) goto 1000 167 : 1240 continue 168 : * le pivot est il correct 169 : * xmatri=imattt(iel) 170 : * segact xmatri 171 : remax=abs(re(ince,1,iel)) 172 : do 1270 ilig=2,re(/1) 173 : remax=max(remax,abs(re(ilig,1,iel))) 174 : 1270 continue 175 : remax=remax*1e-2 176 : if (abs(re(ince,1,iel)).lt.remax) goto 1000 177 : * segdes xmatri 178 : * ok on pivote si il y a lieu 179 : * write (6,*) ' ipp incpp icomb ',ipp,incpp,icomb(ipp,incpp) 180 : if (icomb(ipp,incpp).eq.1) goto 1010 181 : if (incf.eq.0) incf=ince 182 : incf=ince 183 : goto 1000 184 : 1000 continue 185 : ince=2 186 : if (incf.ne.0) ince=incf 187 : 1010 continue 188 : if (ince.ne.2) then 189 : * on pivote 2 et ince 190 : * write (6,*) ' on pivote avec la ligne ',ince 191 : nbpvt=nbpvt+1 192 : incp=posinc(ince) 193 : posinc(ince)=posinc(2) 194 : posinc(2)=incp 195 : segini,des1=descr 196 : segdes descr 197 : descr=des1 198 : ri1.irigel(3,irig)=descr 199 : ri2.irigel(3,irig)=descr 200 : * write (6,*) noelep(2),noelep(ince),lisinc(2),lisinc(ince) 201 : co1=lisinc(2) 202 : lisinc(2)=lisinc(ince) 203 : lisinc(ince)=co1 204 : co1=lisdua(2) 205 : lisdua(2)=lisdua(ince) 206 : lisdua(ince)=co1 207 : noe=noelep(2) 208 : noelep(2)=noelep(ince) 209 : noelep(ince)=noe 210 : noe=noeled(2) 211 : noeled(2)=noeled(ince) 212 : noeled(ince)=noe 213 : segini,Xmatr1=Xmatri 214 : segdes Xmatri 215 : Xmatri=Xmatr1 216 : ri1.irigel(4,irig)=Xmatri 217 : ri2.irigel(4,irig)=Xmatri 218 : do 1100 im=1,re(/3) 219 : * xmatri=imattt(im) 220 : * segini,xmatr1=xmatri 221 : * segdes xmatri 222 : * xmatri=xmatr1 223 : * imattt(im)=xmatri 224 : do 1130 il=1,re(/1) 225 : ret=re(il,2,im) 226 : re(il,2,im)=re(il,ince,im) 227 : re(il,ince,im)=ret 228 : 1130 continue 229 : do 1160 il=1,re(/2) 230 : ret=re(2,il,im) 231 : re(2,il,im)=re(ince,il,im) 232 : re(ince,il,im)=ret 233 : 1160 continue 234 : * segdes xmatri 235 : 1100 continue 236 : endif 237 : * peut etre une dependance 238 : segini,ipt1=meleme 239 : segini,xmatr1=xmatri 240 : do 210 j=1,num(/2) 241 : ip=ipt1.num(noelep(2),j) 242 : * cette matrice contient t elle des inconnues 243 : if (noelep(/1).le.1) then 244 : ipt1.num(1,j)=0 245 : goto 210 246 : endif 247 : incp=posinc(2) 248 : if (ipass.eq.1.and.icomb(ipp,incpp).ne.1) then 249 : ipt1.num(1,j)=0 250 : nbinc=nbinc+1 251 : goto 210 252 : endif 253 : if (iautr(ip,incp).eq.1) then 254 : ipt1.num(1,j)=0 255 : nbinc=nbinc+1 256 : goto 210 257 : endif 258 : * cette matrice touche t elle une autre relation 259 : do 240 inc=2,lisinc(/2) 260 : incp=posinc(inc) 261 : ip=ipt1.num(noelep(inc),j) 262 : if (ielim(ip,incp).eq.1) then 263 : ipt1.num(1,j)=0 264 : nbinc=nbinc+1 265 : goto 210 266 : endif 267 : 240 continue 268 : * le pivot est il correct 269 : * xmatri=imattt(j) 270 : * segact xmatri 271 : remax=abs(re(2,1,j)) 272 : do 270 ilig=2,re(/1) 273 : remax=max(remax,abs(re(ilig,1,j))) 274 : 270 continue 275 : remax=remax*1e-2 276 : if (abs(re(2,1,j)).lt.remax) then 277 : ipt1.num(1,j)=0 278 : nbpiv=nbpiv+1 279 : goto 210 280 : endif 281 : * segdes xmatri 282 : * ok on la prend en dependance 283 : nbdepe=nbdepe+1 284 : nbdep=nbdep+1 285 : * write (6,*) ' inconnu eliminee ',ipt1.num(noelep(2),j), 286 : * > lisinc(2) 287 : ielim(ipt1.num(noelep(2),j),posinc(2))=1 288 : do 280 inc=2,lisinc(/2) 289 : iautr(ipt1.num(noelep(inc),j),posinc(inc))=1 290 : 280 continue 291 : 210 continue 292 : segsup trav3 293 : * il faut maintenant tasser ipt1 et creer ipt2 294 : nbnn=num(/1) 295 : nbelem=num(/2)-nbdepe 296 : * write(6,*) ' nombre d elements crees pour ipt2 ',nbelem 297 : nbsous=0 298 : nbref=0 299 : segini ipt2 300 : ipt2.itypel=itypel 301 : nelrig=nbelem 302 : nligrd=re(/1) 303 : nligrp=re(/2) 304 : segini xmatr2 305 : idec=0 306 : do 300 j=1,num(/2) 307 : if (ipt1.num(1,j).eq.0) then 308 : idec=idec+1 309 : do 310 i=1,num(/1) 310 : ipt2.num(i,idec)=num(i,j) 311 : 310 continue 312 : do io=1,nligrp 313 : do iu=1,nligrd 314 : xmatr2.re(iu,io,idec)=re(iu,io,j) 315 : enddo 316 : enddo 317 : * imatr2.imattt(idec)=imattt(j) 318 : else 319 : do 320 i=1,num(/1) 320 : ipt1.num(i,j-idec)=num(i,j) 321 : 320 continue 322 : do io=1,nligrp 323 : do iu=1,nligrd 324 : xmatr1.re(iu,io,j-idec)=re(iu,io,j) 325 : enddo 326 : enddo 327 : * imatr1.imattt(j-idec)=imattt(j) 328 : endif 329 : 300 continue 330 : nbelem=nbdepe 331 : segadj ipt1 332 : nelrig=nbelem 333 : segadj xmatr1 334 : ri1.irigel(1,irig)=ipt1 335 : ri2.irigel(1,irig)=ipt2 336 : ri1.irigel(4,irig)=xmatr1 337 : ri2.irigel(4,irig)=xmatr2 338 : segdes ipt1,ipt2,meleme,xmatr1,xmatr2 339 : goto 200 340 : 190 continue 341 : * pas une dependance 342 : ri1.irigel(1,irig)=0 343 : 200 continue 344 : * compression de ri1 345 : idec=0 346 : do 500 irig=1,irigel(/2) 347 : meleme=ri1.irigel(1,irig) 348 : if (meleme.eq.0) then 349 : idec=idec+1 350 : else 351 : segact meleme 352 : if (num(/2).eq.0) then 353 : idec=idec+1 354 : else 355 : do 510 ir=1,irigel(/1) 356 : ri1.irigel(ir,irig-idec)=ri1.irigel(ir,irig) 357 : 510 continue 358 : ri1.coerig(irig-idec)=ri1.coerig(irig) 359 : endif 360 : segdes meleme 361 : endif 362 : 500 continue 363 : nrige=max(irigel(/1),8) 364 : nrigel=irigel(/2)-idec 365 : * write (6,*) ' dimension de ri1 ',nrige,nrigel 366 : segadj ri1 367 : do irig=1,nrigel 368 : ri1.irigel(8,irig)=1 369 : enddo 370 : * compression de ri2 371 : idec=0 372 : do 600 irig=1,irigel(/2) 373 : meleme=ri2.irigel(1,irig) 374 : if (meleme.eq.0) then 375 : idec=idec+1 376 : else 377 : segact meleme 378 : if (num(/2).eq.0) then 379 : idec=idec+1 380 : else 381 : do 610 ir=1,irigel(/1) 382 : ri2.irigel(ir,irig-idec)=ri2.irigel(ir,irig) 383 : 610 continue 384 : ri2.coerig(irig-idec)=ri2.coerig(irig) 385 : endif 386 : segdes meleme 387 : endif 388 : 600 continue 389 : nrige=max(irigel(/1),8) 390 : nrigel=irigel(/2)-idec 391 : * write (6,*) ' dimension de ri2 ',nrige,nrigel 392 : segadj ri2 393 : do irig=1,nrigel 394 : ri2.irigel(8,irig)=0 395 : enddo 396 : * On va voir si on ne peut pas faire pivoter ri2 397 : if (ipass.eq.2) goto 2001 398 : ri1s=ri1 399 : ri2s=ri2 400 : mrigis=mrigid 401 : mrigid=ri2 402 : nbpiv=0 403 : nbinc=0 404 : ipass=2 405 : goto 2000 406 : 2001 continue 407 : mrigid=mrigis 407 : call fusrig(ri1s,ri1,ri1f) 408 : ri1=ri1f 409 : segsup trav1,trav2 410 : if (iimpi.ne.0) then 411 : write (6,*) ' nombre de relations éliminées ',nbdep 412 : write (6,*) ' nombre de relations gardées à cause du pivot ', 413 : > nbpiv 414 : write (6,*) ' nombre de relations gardées car non indépendantes ', 415 : > nbinc 416 : write (6,*) ' nombre de paquets pivotés ', 417 : > nbpvt 418 : endif 419 : * dualisation des multiplicateurs de Lagrange dans ri2 420 : * si on a des conditions unilatérales, on ne dualise pas, ce sera fait 421 : * dans le resou de unilater 422 : iunil=0 423 : segact ri2 424 : do ir=1,ri2.irigel(/2) 425 : if (ri2.irigel(6,ir).ne.0) iunil=1 426 : enddo 427 : if ((iunil.eq.0.or.nounil.eq.1).and.lagdua.ge.0) 427 : $ call dbblx(ri2,lagdua) 428 : * write (6,*) ' matrice ri1 ' 428 : * call prrigi(ri1) 429 : * write (6,*) ' matrice ri2 ' 429 : * call prrigi(ri2) 430 : 431 : return 432 : end 433 : 434 :
© Cast3M 2003 - All rights reserved.
Disclaimer