1 : C CHOMOD SOURCE PV 13/10/31 21:17:34 7850 2 : C 3 : C DECOMPOSITION DE CROUT MODIFIEE POUR LE CALCUL DES MATRICES 4 : C CONDENSEES D'UN SUPER-ELEMENT 5 : C RESULTAT DANS XMATRI 6 : C 7 : SUBROUTINE CHOMOD(MMATRX,NBNNMA,SNTT,SNTO,XMATRI,NLIGRA) 8 : IMPLICIT INTEGER(I-N) 9 : IMPLICIT REAL*8(A-H,O-Z) 10 : * TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY 11 : PARAMETER (LPCRAY=10000) 12 : INTEGER OOOVAL,ILIMITE 13 : SEGMENT ITMASQ(NBLIG) 14 : SEGMENT IMASQ(LMASQ) 15 : external chole3i 16 : C 17 : LOGICAL TEST1,LIMIT 18 : C 19 : C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX 20 : C 21 : -INC CCOPTIO 22 : -INC SMMATRI 23 : -INC SMRIGID 24 : -INC CCASSIS 25 : -INC CCHOLE 26 : C 27 : SEGMENT SNTO 28 : INTEGER NTOTMA(NN) 29 : ENDSEGMENT 30 : C 31 : SEGMENT SNTT 32 : INTEGER NTTMAI(NN) 33 : ENDSEGMENT 34 : C 35 : SEGMENT KIVPO(IIMAX) 36 : SEGMENT KIVLO(IIMAX) 37 : SAVE IPASV 38 : DATA IPASV/0/ 39 : DATA PREC /1D-18/ 40 : * character*8 zen 41 : * equivalence (zen,izen) 42 : * multithread 43 : iposm=0 44 : * zen='CPU'//char(0) 45 : * le=4 46 : * nbthro=0 47 : nvaor=0 47 : * call ooozen(izen,le) 48 : * if (le.eq.1) read(zen,fmt='(i1)') nbthro 49 : * if (le.eq.2) read(zen,fmt='(i2)') nbthro 50 : * if (nbthro.eq.1) nbthro=0 51 : * nbthro=max(nbesc,2*nbthro,1) 52 : * nbthro=min(nbthro,64) 53 : nbthro=max( nbesc,1) 53 : if (nbthro.gt.1) call threadii 54 : nbthr=nbthro 55 : * nouvelle methode de gestion de l'espace memoire necessitee par la para 56 : * memoire vive totale 57 : MACTIT=OOOVAL(1,1) 58 : * un bloc de memoire fera au plus macti/2 59 : nvstrm=0 60 : MMATRI=MMATRX 61 : SEGACT,MMATRI*MOD 62 : PRCHLV=PREC 63 : MILIGN=IILIGN 64 : SEGACT,MILIGN*MOD 65 : INO=ILIGN(/1) 66 : MDIAG=IDIAG 67 : SEGACT,MDIAG*MOD 68 : NBLIG=INO 69 : SEGINI ITMASQ 70 : precc=prec 71 : INC=DIAG(/1) 72 : nvstrm=inc 73 : INCC=INC 74 : MIMIK=IIMIK 75 : MINCPO=IINCPO 76 : SEGACT,MINCPO,MIMIK 77 : IPLUMI=IMIK(/2)*2 +4 78 : IL2=0 79 : IIMAX=IJMAX+IPLUMI 80 : SEGINI KIVPO,KIVLO 81 : NLIGRP=NLIGRA 82 : NLIGRD=NLIGRA 83 : NGMAXZ=NGMAXY 84 : nelrig=1 85 : SEGINI,XMATRI 86 : C 87 : *** ILIMITE = nb de noeuds qui ne sont pas maitres 88 : C 89 : SEGACT,SNTO 90 : SEGACT,SNTT 91 : ILIMITE = INO - (NTOTMA(/1) + NTTMAI(/1)) 92 : TEST1=.FALSE. 93 : IF(ILIMITE.EQ.0) THEN 94 : IL2=0 95 : TEST1=.TRUE. 96 : ENDIF 97 : INEG=0 98 : NBLAG=0 99 : NENSLX=0 100 : NVSTOC=0 101 : NVSTOR=0 102 : 103 : LIMIT=.FALSE. 104 : C 105 : C **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD, 106 : C **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUD 107 : C 108 : C **** LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX 109 : C 110 : 1 CONTINUE 111 : IF (TEST1) LIMIT=.TRUE. 112 : C 113 : IVALMA=IJMAX+IPLUMI 114 : IL1=IL2+1 115 : IMIN=IL1 116 : 117 : mactic=0 118 : DO 2 I=IL1,INO 119 : VALREG=1D50 120 : LLIGN=ILIGN(I) 121 : SEGACT /ERR=3/LLIGN 122 : NA=IMMMM(/1) 123 : NBPAR=NA+1 124 : NVALL=NJMAX 125 : * IF(IVALMA+NVALL.GT.NGMAXZ.AND.I.GT.IL1 ) GO TO 3 126 : SEGINI /ERR=3/ LIGN 127 : * recuperer la longueur du segment 128 : mactic=mactic+OOA(OOA(OOT+lign)+OOT-1) 129 : if ((mactic.gt.mactit/2.or.mactic.gt.mactit-nvstrm*12) 130 : > .AND.I.GT.IL1) THEN 131 : SEGSUP LIGN 132 : GOTO 3 133 : endif 134 : NVSTOC=NVSTOC + NVALL 135 : IVALMA=IVALMA+NVALL 136 : nvaor=nvaor+xxva(/1) 137 : C 138 : C **** DECOMPACTAGE 139 : C 140 : IPA=1 141 : DO 121 JPA=1,NA 142 : IVPO(JPA)=IPA 143 : KPA=IPPO(JPA+1)-IPPO(JPA) 144 : IPP=IPPO(JPA) 145 : IPPVV(JPA)=IPA-1 146 : LPA=LDEB(JPA) 147 : LPA1=LPA-IPA 148 : DO 122 MPA=1,KPA 149 : LL=LINC(MPA+IPP) 150 : IPLA=LL-LPA1 151 : VAL(IPLA)=XXVA(MPA+IPP) 152 : 122 CONTINUE 153 : * 154 : * REDUCTION EVENTUELLE DE LA LIGNE A SA PARTIE SIGNIFICATIVE 155 : * > 1E-15*DIAGONALE 156 : * 157 : VALREF=ABS(VAL(IPA+IMMMM(JPA)-LPA))*prec 158 : VALREG=MIN(VALREG,VALREF) 159 : ITROP=0 160 : IPAHAU=IPA+IMMMM(JPA)-LPA 161 : DO 123 ILD=IPA,IPAHAU-1 162 : IF (ABS(VAL(ILD)).GT.VALREF) GOTO 124 163 : ITROP=ITROP+1 164 : 123 CONTINUE 165 : 124 CONTINUE 166 : IF (ITROP.NE.0) THEN 167 : LPA=LPA+ITROP 168 : JILD=IPAHAU-ITROP 169 : DO 125 ILD=IPA,JILD 170 : VAL(ILD)=VAL(ILD+ITROP) 171 : 125 CONTINUE 172 : DO 126 ILD=JILD+1,JILD+ITROP 173 : VAL(ILD)=0.D0 174 : 126 CONTINUE 175 : NVALL=NVALL-ITROP 176 : NVSTOC=NVSTOC-ITROP 177 : IVALMA=IVALMA-ITROP 178 : ENDIF 179 : IPA=IPA+IMMMM(JPA)-LPA+1 180 : IMMM(JPA)=IPNO(LPA) 181 : IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA) 182 : 121 CONTINUE 183 : IF (NVALL.NE.VAL(/1)) SEGADJ LIGN 184 : LMASQ=NVALL/MASDIM+1 185 : SEGINI IMASQ 186 : DO ILD=1,VAL(/1) 187 : IF (ABS(VAL(ILD)).GT.VALREG) IMASQ(ILD/MASDIM+1)=1 188 : ENDDO 189 : ITMASQ(I)=IMASQ 190 : **** **** **** 191 : IPREL=IMMMM(1) 192 : IDERL=IMMMM(NA) 193 : IPPVV(NA+1)=IPA-1 194 : SEGSUP LLIGN 195 : ILIGN(I)=LIGN 196 : IF(IIMPI.EQ.1525) THEN 197 : WRITE( IOIMP,4987) I 198 : 4987 FORMAT (' NOEUD NUMERO ',I5) 199 : LL=VAL(/1) 200 : WRITE(IOIMP, 4926)(VAL(MPA),MPA=1,LL) 201 : 4926 FORMAT(' VAL ' , 10E11.4) 202 : ENDIF 203 : nbths=nbthr 204 : if (mactic/OOA(OOA(OOT+lign)+OOT-1).ge.(nbths*2)/3.and. 205 : > i-il1+1.ge.2*nbths.and.mod(i-il1+1,nbths).eq.0 206 : > .and.mactic.gt.NGMAXZ*nbths) then 207 : il2=i 208 : GOTO 4 209 : endif 210 : C 211 : IF ((I.EQ.(ILIMITE)).AND.(TEST1.EQV..FALSE.)) THEN 212 : TEST1=.TRUE. 213 : IL2 = I 214 : GOTO 4 215 : ENDIF 216 : 2 CONTINUE 217 : IL2=INO 218 : GO TO 4 219 : 3 IL2=I-1 220 : 4 CONTINUE 221 : C 222 : IF(IL2.GE.IL1) GO TO 40 223 : C 224 : C **** APPEL AUX ERREURS MESSAGE PAS ASSEZ DE PLACE MEMOIRE 225 : C 226 : ITYP=48 226 : CALL ERREUR(48) 226 : if (nbthr.gt.1) call threadis 227 : RETURN 228 : 40 CONTINUE 229 : IM=IL2+1 230 : DO 352 IH=IL2,IL1,-1 231 : LIGN=ILIGN(IH) 232 : DO 354 JH=1,IMMM(/1) 233 : IM=MIN(IM,IMMM(JH)) 234 : 354 CONTINUE 235 : IMM=IM 236 : 352 CONTINUE 237 : 353 CONTINUE 238 : LIGN=ILIGN(IL1) 239 : IL11=IPREL 240 : C 241 : C **** BOUCLE *5* TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE 242 : C 243 : * if (ngmaxz.eq.ngmaxy) write (6,*) 'resolution in core ',il1,il2 244 : * if (ngmaxz.ne.ngmaxy) write (6,*) 'resolution out of core ', 245 : * > il1,il2 246 : lig1=ilign(imin) 247 : ipos=0 248 : iper=imin 249 : ider=imin-1 250 : ngmaxz=ngmaxy 251 : DO 5 I=IMIN,IL2 252 : IMASQ=ITMASQ(I) 253 : LIG1=ILIGN(I) 254 : IF(I.LT.IL1) GO TO 7 255 : C____________ 256 : C 257 : C ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A 258 : C ******* IPREL IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER 259 : C ******* LE TERME DIAGONAL 260 : C 261 : LIGN=LIG1 262 : DO 156 KHG=1,IMMM(/1) 263 : II=IPREL-1+KHG 264 : IMMM(KHG)=0 265 : NN=IPPVV(KHG+1) 266 : NNM1=IPPVV(KHG) 267 : N=NN-NNM1 268 : DIAG(II)=VAL(NN) 269 : IF(N.EQ.1) THEN 270 : IF(LIMIT) THEN 271 : RE(II-NBNNMA,II-NBNNMA,1)=VAL(NN) 272 : GO TO 41 273 : ELSE 274 : GO TO 8 275 : ENDIF 276 : ENDIF 277 : NMI=N-II 278 : KI1=N-1 279 : KQ=-NMI 280 : IF(LIMIT) THEN 281 : C 282 : C ****** CAS NOEUD MAITRE 283 : C 284 : IDEP=MAX(IL11,1-NMI) 285 : KIDEP=IDEP+NMI 286 : VAL(NN)=VAL(NN)+ 287 : # SUPER1(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI), 288 : # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,NBNNMA,XMATRI, 289 : # IMASQ(1),1+IPPVV(KHG)) 290 : RE(II-NBNNMA,II-NBNNMA,1)=VAL(NN) 291 : GOTO 41 292 : ENDIF 293 : C 294 : C ****** CAS NOEUD NON MAITRE 295 : C 296 : IDEP=MAX(IL11,2-NMI) 297 : KIDEP=IDEP+NMI 298 : VAL(NN)=VAL(NN)+ 299 : # CHOLE1(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI), 300 : # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG), 301 : # prec) 302 : C 303 : 8 CONTINUE 304 : IF(ITTR(II).EQ.0.AND. 305 : & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-10) GO TO 12 306 : IF(ITTR(II).NE.0.AND. 307 : & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-6 ) GO TO 12 308 : C 309 : C **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE 310 : C **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT 311 : C **** AU PREALABLE SUR CETTE INCONNUE. 312 : C 313 : * write (6,*) ' chomod nens ittr val diag ',ittr(ii),nens+1, 314 : * > val(nn),diag(ii),ii,ipno(ii),ilimite 315 : if (ittr(ii).eq.0) then 316 : VAL(NN)=DIAG(II)*1d0 317 : else 318 : VAL(NN)=val(nn)-4.D0*DIAG(II) 319 : NENSLX=NENSLX+1 320 : endif 321 : NENS=NENS+1 322 : IMMM(KHG)=NENS 323 : 12 CONTINUE 324 : IMASQ(NN/MASDIM+1)=1 325 : DIAG(II)=VAL(NN) 326 : IF(DIAG(II).NE.0.d0) GO TO 41 327 : KQ1=1+NNM1 328 : KQN=N+NNM1 329 : DO 16 LFG=KQ1,KQN 330 : IF(VAL(LFG).NE.0.d0) GO TO 17 331 : 16 CONTINUE 332 : DIAG(II)=1.d0 333 : if (ittr(ii).ne.0) diag(ii)=-1.d0 334 : VAL(NN)=diag(ii) 335 : GO TO 41 336 : 17 CONTINUE 337 : diag(ii)=val(lfg) 338 : if (ittr(ii).ne.0) diag(ii)=-abs(diag(ii)) 339 : val(nn)=diag(ii) 340 : goto 41 341 : C 342 : C **** ENVOI ERREUR MATRICE SINGULIERE 343 : C 344 : ITYP=49 345 : INTERR(1)=I 345 : CALL ERREUR(49) 345 : if (nbthr.gt.1) call threadis 346 : RETURN 347 : C 348 : C **** ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS 349 : C ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE 350 : C 351 : 41 IF(DIAG(II).LT.0.D0) INEG=INEG+1 352 : IF(ITTR(II).NE.0) NBLAG=NBLAG+1 353 : 156 CONTINUE 354 : NA=IMMM(/1) 355 : * 356 : * RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE) 357 : * 357 : CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1), 358 : # NVALL,IPPVV(1),IZROSF,NA,PREC,imasq(1)) 359 : * 360 : IF (NBPAR.NE.NA+1.OR.NVALL.NE.VAL(/1)) SEGADJ LIGN 361 : NVSTOR=NVSTOR+NVALL 362 : nvstrm=max(nvstrm,nvall) 363 : DO 143 LHG=1,NBPAR 364 : IVPO(2*LHG-1)=KIVPO(LHG) 365 : IVPO(2*LHG)=KIVLO(LHG) 366 : 143 CONTINUE 367 : segsup imasq 368 : imasq=itmasq(i) 369 : if (i.gt.1) then 370 : lig1=ilign(i-1) 371 : segdes lig1 372 : endif 373 : C 374 : C **** ON TRIANGULARISE LES AUTRES LIGNES 375 : C 376 : il1=il1+1 377 : if (il1.gt.il2) goto 5 378 : LIG1=ILIGN(I) 379 : lign=ilign(il1) 380 : IL11=IPREL 381 : goto 7 382 : 72 continue 383 : * write (6,*) ' chomod segact ',i,ILIMITE,il1,il2 384 : 71 continue 385 : if (i.ne.il1-1) ngmaxz=2 000 000 000 386 : * soit parce qu'on a fini, soit parce qu'on manque de memoire 387 : * il faut executer les lignes activees puis les desactiver 388 : * lancer les chole3 et attendre qu'ils soient finis 389 : if (ipos.ne.0) then 390 : * write (6,*) ' lancement thread ',iper,ider,il1,il2 391 : if (iper.gt.ider) then 392 : write (6,*) ' erreur interne chole ' 392 : call erreur(5) 393 : endif 394 : IF(I.LE.ILIMITE) THEN 395 : if (nbthr.gt.1) then 396 : do ith=1,min(il2-il1+1,nbthr-1) 396 : call threadid(ith,chole3i) 397 : enddo 397 : call chole3i(nbthr) 398 : do ith=1,min(il2-il1+1,nbthr-1) 398 : call threadif(ith) 399 : enddo 400 : * en multithread il peut y avoir n'importe quoi dans oov(1) du 401 : * aux acces simultan�s et ca crache gemat. donc : 402 : oov(1)=0 403 : else 404 : do ith=1,nbthr 404 : call chole3i(ith) 405 : enddo 406 : endif 407 : endif 408 : iposm=max(iposm,ipos) 409 : ipos=0 410 : iderf=ider-1 411 : if (ider.ne.il1-1) iderf=ider 412 : do il=iderf,iper,-1 413 : lign=ilign(il) 414 : segdes lign 415 : * write (6,*) ' desactivation ligne ',il 416 : enddo 417 : iper=ider+1 418 : * write (6,*) ' iper ider il1 ',iper,ider,il1 419 : if (i.eq.ilimite) goto 5 420 : if (iper.ne.il1) goto 7 421 : ENDIF 422 : goto 5 423 : 7 CONTINUE 424 : SEGACT/err=72/LIG1 425 : ipos=ipos+1 426 : ider=i 427 : if (i.eq.il1-1.and.i.le.ilimite) goto 71 428 : IF(I.EQ.ILIMITE) GOTO 71 429 : IF(I.GT.ILIMITE) THEN 430 : IPPR=LIG1.IPREL 431 : IDDR=LIG1.IDERL 432 : C 433 : C 434 : DO 10 JBI=IL1,IL2 435 : LIGN=ILIGN(JBI) 436 : IMASQ=ITMASQ(JBI) 437 : C 438 : C ****** CAS NOEUD MAITRE 439 : C 439 : CALL SUPER3(IPREL,IDERL,IPPVV(1),IPPR,IDDR, 440 : # IVPO(1),LIG1.IPPVV(1),VAL(1), 441 : # LIG1.VAL(1),LIG1.IVPO(1),NBNNMA,XMATRI, 442 : # imasq(1)) 443 : IF(IMM.GT.I) GO TO 106 444 : 10 CONTINUE 445 : 106 CONTINUE 446 : IF(I.LT.IL1) then 447 : SEGDES,LIG1 448 : endif 449 : ENDIF 450 : 5 CONTINUE 451 : DO 11 I=IL1,IL2 452 : LIGN=ILIGN(I) 453 : SEGDES,LIGN 454 : 11 CONTINUE 455 : IF(IL2.LT.INO) GO TO 1 456 : 457 : C ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF 458 : C ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE 459 : c INEG=INEG-NBLAG 460 : * on ne compte pas 2 fois les multiplicateurs qui vont etre 461 : * eliminés lors de la resolution car mode d'ensemble 462 : INEG=INEG-(NBLAG-NENSLX) 463 : if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX 464 : 4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS', 465 : 1' DE LAGRANGE DETECTES') 466 : 467 : IF (IIMPI.EQ.1)WRITE(IOIMP,4821) NVSTOC 468 : 4821 FORMAT( ' CHOMOD NOMBRE DE VALEURS DANS LE PROFIL',I9) 469 : IF (IIMPI.EQ.1)WRITE(IOIMP,4822) NVSTOR 470 : 4822 FORMAT( ' CHOMOD NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I9) 471 : IF (IIMPI.EQ.1)WRITE(IOIMP,4823) NVaor 472 : 4823 FORMAT( ' CHOMOD NOMBRE DE VALEURS initiales',I9) 473 : INTERR(1)=NVSTOR 474 : reaerr(1)=nvstor/inc**(4./3) 474 : if (ipasv.eq.0) CALL ERREUR(-278) 475 : IPASV=1 476 : SEGDES,MINCPO 477 : SEGDES,MIMIK 478 : SEGDES,MMATRI 479 : SEGDES,MILIGN 480 : SEGDES,MDIAG 481 : SEGDES,SNTO 482 : SEGDES,SNTT 483 : SEGDES,XMATRI 484 : MMATRX=MMATRI 485 : SEGSUP KIVPO,KIVLO,ITMASQ 485 : if (nbthr.gt.1) call threadis 486 : RETURN 487 : END 488 : 489 : 490 : 491 : 492 : 493 : 494 : 495 : 496 : 497 :
© Cast3M 2003 - All rights reserved.
Disclaimer