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