1 : C GRACO8 SOURCE PV 11/04/11 21:15:06 6900 2 : C 3 : SUBROUTINE GRACO8(MMATRX,MVECTX,NOEN) 4 : C 5 : C **** EXECUTE LA SOLUTION X DE (Lt D L) X=F 6 : C 7 : CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB 8 : CMB 9 : CMB Plutot la solution de L.D.Lt ou L.D.Mt (cas non symétrique) 10 : CMB Elle devrait dons s'appeller DESMON et non MONDES. 11 : CMB 12 : CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB 13 : C 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : PARAMETER (LPCRAY=10000) 17 : INTEGER OOOVAL 18 : 19 : C ... Variable décrivant l'EXIStence du Triangle Supérieur différent 20 : C de l'inférieur (cas non symétrique) ... 21 : LOGICAL EXISTS 22 : 23 : -INC SMMATRI 24 : -INC SMELEME 25 : -INC SMVECTD 26 : -INC CCOPTIO 27 : 28 : SEGMENT ITTRV(MULRE) 29 : SEGMENT,ITRAA(NENS) 30 : SEGMENT/BID/(BIDON(IIMAX+10)*D) 31 : 32 : if (noen.eq.0) then 33 : nbelem=0 34 : nbnn=1 35 : nbsous=0 36 : nbref=0 37 : segini ipt8 38 : ipt8.itypel=1 39 : endif 39 : CALL OOOMRU(1) 40 : IVALMA=0 41 : IF(IIMPI.EQ.3) WRITE(IOIMP,1000)MMATRX,MVECTX 42 : 1000 FORMAT(' SUBROUTINE MONDES : POINTEUR DE LA MATRICE=',I5, 43 : 1 ' POINTEUR DU VECTEUR=',I5) 44 : C 45 : C **** ACTIVATION DES SEGMENTS 46 : C 47 : MMATRI=MMATRX 48 : SEGACT,MMATRI*MOD 49 : * BEC00SE OPTIMISEUR 50 : ITRAA=MMATRI 51 : IF(NENS.NE.0) SEGINI,ITRAA 52 : NE1=NENS 53 : IAA=0 54 : 55 : EXISTS=.FALSE. 56 : IF(IILIGS.NE.0) THEN 57 : MILIGN=IILIGS 58 : SEGACT MILIGN 59 : IF(ILIGN(/1).GT.0) EXISTS=.TRUE. 60 : ENDIF 61 : 62 : MILIGN=IILIGN 63 : SEGACT,MILIGN 64 : INC=IPNO(/1) 65 : 66 : MVECTD=MVECTX 67 : SEGACT MVECTD*MOD 68 : 69 : MDNOR=IDNORM 70 : SEGACT MDNOR 71 : IF(IDNORD.GT.0) THEN 72 : MDNO1=IDNORD 73 : SEGACT MDNO1 74 : ELSE 75 : MDNO1=MDNOR 76 : ENDIF 77 : 78 : MDIAG=IDIAG 79 : SEGACT,MDIAG 80 : 81 : C ... MULRE = nombre de seconds membres ... 82 : MULRE = VECTBB(/1)/INC 83 : C ... MUINC ne servira que comme une borne sur les indices de boucles ... 84 : MUINC = ( MULRE-1)*INC 85 : SEGINI ITTRV 86 : 87 : C ... Multiplication du second membre par les DNOR(.) (à cause du 88 : C conditionnement des matrices de blocages) ... 89 : DNORMA=0.D0 90 : DO 450 K=0,MUINC,INC 91 : DO 45 I=1,INC 92 : *** VECTBB(I+K)=VECTBB(I+K)*MDNO1.DNOR(I) 93 : DNORMA=MAX(DNORMA,ABS(VECTBB(I+K))) 94 : 45 CONTINUE 95 : 450 CONTINUE 96 : DNORMA=DNORMA*1D-28 97 : 98 : C ... Détection du premier terme important du second membre ... 99 : II=0 100 : DO 451 K=0,MUINC,INC 101 : II=II+1 102 : DO 452 I=1,INC 103 : IF(ABS(VECTBB(I+K)).LT.DNORMA) GO TO 452 104 : C ... Le numéro du noeud concerné par le premier terme important va à ITTRV ... 105 : ITTRV(II)=IPNO(I) 106 : GO TO 451 107 : 452 CONTINUE 108 : 451 CONTINUE 109 : 110 : C ... NNOE = nombre de noeuds concernés ... 111 : NNOE=ILIGN(/1) 112 : IDEP=NNOE 113 : C ... On cherche le minimum (non nul) des ITTRV que l'on met dans IDEP ... 114 : DO 453 I=1,MULRE 115 : IF(ITTRV(I).LT.IDEP) THEN 116 : IF(ITTRV(I).NE.0) THEN 117 : IDEP=ITTRV(I) 118 : ELSE 119 : IDEP=1 120 : GO TO 4530 121 : ENDIF 122 : ENDIF 123 : 453 CONTINUE 124 : 4530 CONTINUE 125 : SEGDES MDNOR 126 : C 127 : C **** DESCENTE: ON RESOU L*C=B. EN FAIT ON STOCKE C DANS B 128 : C 129 : LTRK=MAX(1+LPCRAY,OOOVAL(1,4)) 130 : IIMAX=(((IJMAX +LTRK)/LTRK)+1)*LTRK 131 : * test si la place disponible permet de tout mettre en memoire 132 : xplds=oooval(1,1)-oooval(2,3) 133 : if (real(nnoe)*ijmax.lt.xplds) iimax=0 134 : *** SEGINI BID 135 : C ... IDEP : on commence par ce noeud car le second membre est nul pour 136 : C tous les précédents ... 137 : DO 610 IVS=IDEP,NNOE 138 : LIGN=ILIGN(IVS) 139 : SEGACT /ERR=50/ LIGN 140 : IVALMA=IVALMA+VAL(/1) 141 : *pvpv IF (IVALMA.GT.NGMAXY) GOTO 50 142 : NA=IMMM(/1) 143 : IPRELL=IPREL 144 : DO 611 IA=1,NA 145 : C ... Si l'inconnue présente un mode d'ensemble ... 146 : IF(IMMM(IA).NE.0) THEN 147 : C ... On incrémente le compteur et 148 : IAA=IAA+1 149 : C ... on met le N° de la ligne dans ITRAA ... 150 : ITRAA(IAA)=IPRELL+IA-1 151 : ENDIF 152 : 611 CONTINUE 153 : C ... IFIB sera vu dans MONDE2 par le COMMON/CMOND2/ ... 154 : IFIB=IVPO(/1) 154 : CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1), 155 : > NA,IPREL,MULRE,INC,IVS,LLOM,IFIB) 156 : 610 CONTINUE 157 : C ... Si on n'a pas eu de pb de mémoire on passe par là ... 158 : *** SEGSUP BID 159 : C ... ILMAX vaut le dernier noeud qui a pu être traité ... 160 : ILMAX=NNOE 161 : C ... On va à la division par la diagonale ... 162 : GOTO 55 163 : 164 : 50 CONTINUE 165 : C ... Si on est là, c'est à cause des pb de mémoire ... 166 : *** SEGSUP BID 167 : C ... IVS est le N° du LIGN qui n'a pas pu être traité ... 168 : C ... ILMAX = N° du dernier traité ... 169 : ILMAX=IVS-1 170 : C ... IILMAX = N° du premier à traiter ... 171 : IILMAX=IVS 172 : DO 54 IVS=IILMAX,NNOE 173 : 58 CONTINUE 174 : LIGN=ILIGN(IVS) 175 : C ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ... 176 : SEGACT /ERR=56/ LIGN 177 : GOTO 57 178 : 56 CONTINUE 179 : * EN CAS DE PROBLEME ON FAIT UN PEU DE PLACE 180 : C ... Si on a toujours des pb de mémoire, on désactive le LIGN N° ILMAX, 181 : LIGN=ILIGN(ILMAX) 182 : SEGDES LIGN*(NOMOD,MRU) 183 : C ... puis ce dernier est décrémenté ... 184 : ILMAX=ILMAX-1 184 : IF (ILMAX.EQ.1) CALL ERREUR(5) 185 : C ... et on recommence ... 186 : GOTO 58 187 : C ... On est là si tout s'est bien passé avec SEGACT ... 188 : 57 CONTINUE 189 : NA=IMMM(/1) 190 : IPRELL=IPREL 191 : DO 612 IA=1,NA 192 : IF(IMMM(IA).EQ.0) GOTO 612 193 : IAA=IAA+1 194 : ITRAA(IAA)=IPRELL+IA-1 195 : 612 CONTINUE 196 : IFIB = IVPO(/1) 196 : CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1), 197 : > NA,IPREL,MULRE,INC,IVS,LLOM,IFIB) 198 : SEGDES,LIGN*(NOMOD,MRU) 199 : 54 CONTINUE 200 : 55 CONTINUE 201 : C 202 : C ... À cet endroit la descente est finie. L'état des LIGN est le suivant : 203 : C ... Ceux de IDEP à ILMAX sont actifs, les autres (si ILMAX < NNOE) sont désactivés ... 204 : C 205 : C **** DIVISION PAR LE TERME DIAGONAL **** 206 : C 207 : ** vmin=1D50 208 : ** vmax=0D0 209 : DO 120 K=0,MUINC,INC 210 : DO 12 I=1,INC 211 : J=I+K 212 : VECTBB(J)=VECTBB(J)/DIAG(I) 213 : ** vmin=min(vmin,abs(vectbb(j))) 214 : ** vmax=max(vmax,abs(vectbb(j))) 215 : 12 CONTINUE 216 : ** write (6,*) ' graco8 vmin vmax ',vmin,vmax 217 : 120 CONTINUE 218 : SEGDES,MDIAG 219 : 220 : C 221 : C **** MONTEE **** 222 : C 223 : C ... Si la matrice est asymétrique ... 224 : IF(EXISTS) THEN 225 : 226 : C ... On commence par désactiver les LIGN qui sont restés actifs ... 227 : DO 2000 I=IDEP,ILMAX 228 : LIGN=ILIGN(I) 229 : SEGDES,LIGN*(NOMOD,MRU) 230 : 2000 CONTINUE 231 : 232 : C ... Puis on désactive MILIGN ... 233 : SEGDES,MILIGN 234 : 235 : C ... On change de MILIGN ... 236 : MILIGN=IILIGS 237 : SEGACT,MILIGN 238 : 239 : C ... Et on active des LIGN pointés par ce dernier ... 240 : ILMAX=IDEP-1 241 : DO 2100 I=IDEP,NNOE 242 : LIGN=ILIGN(I) 243 : SEGACT /ERR=2110/ LIGN 244 : ILMAX = I 245 : 2100 CONTINUE 246 : C ... Et on passe à la montée ... 247 : GOTO 3000 248 : 249 : C ... Si on n'a même pas activé le premier, alors ADIOS ... 249 : 2110 IF(ILMAX.EQ.IDEP-1) CALL ERREUR(5) 250 : 251 : 3000 CONTINUE 252 : ENDIF 253 : 254 : 255 : J=NNOE 256 : C ... Début de la pseudo boucle sur les LIGN qui sont restés désactivés ... 257 : 70 CONTINUE 258 : IF (J.EQ.ILMAX) GOTO 72 259 : LIGN=ILIGN(J) 260 : SEGACT /ERR=68/ LIGN 261 : GO TO 67 262 : C ... Si on a des pb avec activation, on désactive et on diminue ILMAX ... 263 : 68 CONTINUE 264 : LIGN = ILIGN( ILMAX) 265 : SEGDES LIGN*(NOMOD,MRU) 266 : ILMAX=ILMAX-1 266 : IF(ILMAX.EQ.1) CALL ERREUR (5) 267 : C ... puis on recommence la tentative ... 268 : GO TO 70 269 : C ... On a réussi à activer ... 270 : 67 CONTINUE 271 : NA=IMMM(/1) 272 : IPRELL=IPREL 273 : IFIB=IVPO(/1) 273 : CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1), 274 : > NA,IPREL,MULRE,INC,IFIB) 275 : SEGDES,LIGN*(NOMOD,MRU) 276 : J = J-1 277 : GO TO 70 278 : 72 CONTINUE 279 : CC FIN DE PSEUDO BOUCLE J = INC ,ILMAX+1,-1 - Vieux commentaire et Faux !!! 280 : CC FIN DE PSEUDO BOUCLE J = NNOE ,ILMAX+1,-1 281 : 282 : C ... Dans cette boucle on parcourt les LIGN qui sont restés actifs ... 283 : DO 13 J=ILMAX,IDEP,-1 284 : LIGN=ILIGN(J) 285 : NA=IMMM(/1) 286 : IPRELL=IPREL 287 : IFIB=IVPO(/1) 287 : CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1), 288 : > NA,IPREL,MULRE,INC,IFIB) 289 : *pvpv SEGDES,LIGN*(NOMOD,MRU) 290 : 13 CONTINUE 291 : 292 : C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ... 293 : DO 113 J=IDEP-1,1,-1 294 : LIGN=ILIGN(J) 295 : SEGACT LIGN 296 : NA=IMMM(/1) 297 : IPRELL=IPREL 298 : DO 1131 ILM=1,NA 299 : IF(IMMM(ILM).EQ.0) GO TO 1131 300 : IAA=IAA+1 301 : ITRAA(IAA)=IPRELL+ILM-1 302 : 1131 CONTINUE 303 : IFIB=IVPO(/1) 303 : CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1), 304 : > NA,IPREL,MULRE,INC,IFIB) 305 : SEGDES,LIGN*(NOMOD,MRU) 306 : 113 CONTINUE 307 : 308 : C 309 : C **** ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE 310 : C 311 : C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ... 312 : NENS=IAA 313 : NBEN=0 314 : IF(NENS.EQ.0) GO TO 26 315 : 316 : MINCPO=IINCPO 317 : MIMIK=IIMIK 318 : SEGACT,MINCPO,MIMIK 319 : MELEME=IGEOMA 320 : SEGACT,MELEME 321 : NNOE=INCPO(/2) 322 : IINC1=INCPO(/1) 323 : 324 : C ... Boucle sur les seconds membres ... 325 : DO 300 KI=0,MUINC,INC 326 : C ... XMA sera le maximum des val. abs. des termes (multiplicateurs 327 : C compris) du résultat partiel (avant la multiplication par MDNOR) 328 : C N° KI+1 multiplié par 1.e-10 ... 329 : XMA=0. 330 : DO 30 KK=1,INC 331 : I=KK+KI 332 : XMA=MAX(XMA,ABS(VECTBB(I))) 333 : 30 CONTINUE 334 : XMA=XMA*1.D-28 335 : 336 : C ... Boucle sur les modes d'ensemble ... 337 : ** write (6,*) ' nens dans graco8 ',nens 338 : DO 21 IA=1,NENS 339 : I1=ITRAA(IA) 340 : J=IPNO(I1) 341 : DO 22 K=1,IINC1 342 : IF(INCPO(K,J).EQ.I1) GO TO 23 343 : 22 CONTINUE 343 : CALL ERREUR(5) 344 : RETURN 345 : 346 : 23 CONTINUE 347 : C ... Si ce n'est pas un multiplicateur, le déplacement doit être 348 : C < à XMA, sinon ERREUR ... 349 : IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA) then 350 : vectbb(i1+ki)=0.d0 351 : GO TO 20 352 : endif 353 : C ... Si c'est un multiplicateur, le multiplicateur doit être 354 : C < à 1.e+7*XMA, sinon ERREUR ... 355 : if (ittr(i1).ne.0) then 356 : i2=ittr(i1) 357 : if (abs(vectbb(i1+ki)+vectbb(i2+ki)).le.xma*1d5 ) then 358 : vectbb(i1+ki)=0.d0 359 : vectbb(i1+ki)=0.d0 360 : goto 20 361 : endif 362 : endif 363 : 364 : * IF(ITTR(I1).NE.0.AND.ABS(VECTBB(I1+KI)).LE.XMA*1E4) then 365 : * vectbb(i1+ki)=0.d0 366 : * GO TO 20 367 : * endif 368 : C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme 369 : C ... Si option NOEN on accumule les pts dans un meleme 370 : IF(NOEN.EQ.0) THEN 371 : nbelem=ipt8.num(/2)+1 372 : segadj ipt8 373 : ipt8.num(1,nbelem)=num(1,j) 374 : ENDIF 375 : NBEN=NBEN+1 376 : MOTERR(1:4)=IMIK(K) 377 : INTERR(1)=NUM(1,J) 378 : IF(NOEN.EQ.0) THEN 379 : GO TO 21 380 : ENDIF 381 : ITY=149 381 : CALL ERREUR(ITY) 382 : RETURN 383 : 384 : 20 CONTINUE 385 : C ... Messages d'information ... 386 : JJK=NUM(1,J) 387 : KIK=KI/INC +1 388 : 389 : IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1) then 390 : * if (imik(k).ne.'LX ') 391 : WRITE(IOIMP,41) JJK,IMIK(K) 392 : endif 393 : 394 : IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) 395 : & WRITE(IOIMP,42) JJK,IMIK(K) 396 : 397 : IF(ITTR(I1).NE.1 .AND. MULRE.NE.1) 398 : & WRITE(IOIMP,410)KIK,JJK,IMIK(K) 399 : 400 : IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1) 401 : & WRITE(IOIMP,420)KIK,JJK,IMIK(K) 402 : 403 : 21 CONTINUE 404 : 405 : 300 CONTINUE 406 : 407 : 41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ', 408 : * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ', 409 : * 'LA SUSDITE') 410 : 42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ', 411 : * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ', 412 : * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE ') 413 : 410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD ' 414 : *,I6,' INCONNUE ', 415 : * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ', 416 : * 'LA SUSDITE') 417 : 420 FORMAT(' VECTEUR NUMERO ',I3,/, 418 : * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ', 419 : * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ', 420 : * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE ') 421 : 422 : SEGDES,MELEME 423 : SEGDES,MINCPO 424 : SEGDES,MIMIK 425 : 26 CONTINUE 426 : * indiquer si besoin le nombre de modes d'ensemble excites 427 : * et le maillage des noeuds freinés 428 : IF (NOEN.EQ.0) then 429 : * write (6,*) ' mondes nben ',nben 429 : CALL ECRENT(NBEN) 429 : call ecrobj('MAILLAGE',ipt8) 430 : segdes ipt8 431 : endif 432 : MDNOR=IDNORM 433 : SEGACT,MDNOR 434 : DO 350 K=0,MUINC,INC 435 : DO 35 I = 1, INC 436 : *** VECTBB(I+K)=VECTBB(I+K)*DNOR(I) 437 : 35 CONTINUE 438 : * on force l'egalite des multiplicateurs de lagrange 439 : * en tenant compte de dnor 440 : DO 36 I = 1, INC 441 : if (ITTR(I).ne.0) then 442 : * write (6,*) ' dans mondes ',i,ittr(i) 443 : if (vectbb(i+k).eq.0.d0.or.vectbb(ittr(i)+k).eq.0.d0) then 444 : * write (6,*) ' mondes vectbbs ',vectbb(i+k),vectbb(ittr(i)+k) 445 : vectbb(i+k)=0.d0 446 : vectbb(ittr(i)+k)=0.d0 447 : goto 36 448 : endif 449 : dn1=dnor(i) 450 : dn2=dnor(ittr(i)) 451 : vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k)*dn2/dn1)/2.D0 452 : vectbb(ittr(i)+k)=vectbb(i+k)*dn1/dn2 453 : endif 454 : 36 CONTINUE 455 : 350 CONTINUE 456 : ** SEGDES,MDNOR 457 : C ... On ne désactive MDNO1 que dans le cas où il est 458 : C différent de MDNOR ... 459 : IF(IDNORD.GT.0) THEN 460 : ** SEGDES,MDNO1 461 : ENDIF 462 : ** SEGDES MMATRI 463 : ** SEGDES,MILIGN 464 : ** SEGDES,MVECTD 465 : SEGSUP ITTRV 466 : IF(NE1.NE.0) SEGSUP,ITRAA 467 : 468 : IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD 469 : 1001 FORMAT(' SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5) 470 : 470 : CALL OOOMRU(0) 471 : RETURN 472 : END 473 : C 474 : 475 : 476 : 477 : 478 : 479 : 480 : 481 : 482 : 483 : 484 : 485 : 486 : 487 : 488 : 489 : 490 : 491 : 492 :
© Cast3M 2003 - All rights reserved.
Disclaimer