1 : C CHOLE SOURCE PV 13/10/31 21:17:33 7850 2 : SUBROUTINE CHOLE(MMATRX,PREC) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : * TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY 6 : PARAMETER (LPCRAY=10000) 7 : INTEGER OOOVAL 8 : C 9 : C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX 10 : C 11 : -INC CCOPTIO 12 : -INC SMMATRI 13 : -INC CCASSIS 14 : -INC CCHOLE 15 : SEGMENT KIVPO(IIMAX) 16 : SEGMENT KIVLO(IIMAX) 17 : SEGMENT ITMASQ(NBLIG) 18 : SEGMENT IMASQ(LMASQ) 19 : 20 : external chole3i 21 : SAVE IPASV 22 : DATA IPASV/0/ 23 : * character*8 zen 24 : * equivalence (zen,izen) 25 : iposm=0 26 : * zen='CPU'//char(0) 27 : * le=4 28 : * nbthro=0 29 : nvaor=0 30 : nbthro=max(nbesc,1) 30 : if (nbthro.gt.1) call threadii 31 : nbthr=nbthro 32 : * nouvelle methode de gestion de l'espace memoire necessitee par la parallelisation 33 : * memoire vive totale 34 : MACTIT=OOOVAL(1,1) 35 : * un bloc de memoire fera au plus macti/2 36 : nvstrm=0 37 : MMATRI=MMATRX 38 : SEGACT,MMATRI*MOD 39 : PRCHLV=PREC 40 : MILIGN=IILIGN 41 : SEGACT,MILIGN*MOD 42 : INO=ILIGN(/1) 43 : MDIAG=IDIAG 44 : SEGACT,MDIAG*MOD 45 : NBLIG=INO 46 : SEGINI ITMASQ 47 : precc=prec 48 : INC=DIAG(/1) 49 : nvstrm=inc 50 : INCC=INC 51 : MIMIK=IIMIK 52 : MINCPO=IINCPO 53 : SEGACT,MINCPO,MIMIK 54 : IPLUMI=IMIK(/2)*2 +4 55 : IL2=0 56 : IIMAX=IJMAX+IPLUMI 57 : SEGINI KIVPO,KIVLO 58 : INEG=0 59 : NBLAG=0 60 : NENSLX=0 61 : NVSTOC=0 62 : NVSTOR=0 63 : NGMAXZ=NGMAXY 64 : C 65 : C **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD, 66 : C **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS 67 : C 68 : C **** LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX 69 : C 70 : 1 CONTINUE 71 : IVALMA=IJMAX+IPLUMI 72 : IL1=IL2+1 73 : * write (6,*) ' il1 il2 apres 1 ',il1,il2 74 : IMIN=IL1 75 : mactic=0 76 : DO 2 I=IL1,INO 77 : LLIGN=ILIGN(I) 78 : SEGACT /ERR=3/LLIGN 79 : NA=IMMMM(/1) 80 : NBPAR=NA+1 81 : NVALL=NJMAX 82 : SEGINI /ERR=3/ LIGN 83 : * recuperer la longueur du segment 84 : mactic=mactic+OOA(OOA(OOT+lign)+OOT-1) 85 : if ((mactic.gt.mactit/2.or.mactic.gt.mactit-nvstrm*12) 86 : > .AND.I.GT.IL1) THEN 87 : SEGSUP LIGN 88 : GOTO 3 89 : endif 90 : NVSTOC=NVSTOC + NVALL 91 : IVALMA=IVALMA+NVALL 92 : nvaor = nvaor + xxva(/1) 93 : C 94 : C **** DECOMPACTAGE 95 : C 96 : LMASQ=NVALL/MASDIM+2 97 : SEGINI IMASQ 98 : IPA=1 99 : DO 121 JPA=1,NA 100 : IVPO(JPA)=IPA 101 : KPA=IPPO(JPA+1)-IPPO(JPA) 102 : IPP=IPPO(JPA) 103 : IPPVV(JPA)=IPA-1 104 : LPA=LDEB(JPA) 105 : LPA1=LPA-IPA 106 : DO 122 MPA=1,KPA 107 : LL=LINC(MPA+IPP) 108 : IPLA=LL-LPA1 109 : VAL(IPLA)=XXVA(MPA+IPP) 110 : IMASQ(IPLA/MASDIM+1)=1 111 : 122 CONTINUE 112 : IPA=IPA+IMMMM(JPA)-LPA+1 113 : IMMM(JPA)=IPNO(LPA) 114 : IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA) 115 : 121 CONTINUE 116 : **** **** **** 117 : ITMASQ(I)=IMASQ 118 : IPREL=IMMMM(1) 119 : IDERL=IMMMM(NA) 120 : IPPVV(NA+1)=IPA-1 121 : SEGSUP LLIGN 122 : ILIGN(I)=LIGN 123 : nbths=nbthr 124 : * nb de ligne multiple du nb de threads 125 : if (mod(i+1-il1,nbths).eq.0 126 : > .and.mactic.gt.NGMAXZ*nbths) then 127 : * write (6,*) ' chole il1 il2 ',il1,i 128 : il2=i 129 : GOTO 4 130 : endif 131 : 2 CONTINUE 132 : IL2=INO 133 : GO TO 4 134 : 3 IL2=I-1 135 : 4 CONTINUE 136 : IF(IL2.GE.IL1) GO TO 40 137 : C 138 : C **** APPEL AUX ERREURS MESSAGE PAS ASSEZ DE PLACE MEMOIRE 139 : C 140 : ITYP=48 140 : CALL ERREUR(48) 140 : if (nbthr.gt.1) call threadis 141 : RETURN 142 : 40 CONTINUE 143 : IM=IL2+1 144 : DO 352 IH=IL2,IL1,-1 145 : LIGN=ILIGN(IH) 146 : DO 354 JH=1,IMMM(/1) 147 : IM=MIN(IM,IMMM(JH)) 148 : 354 CONTINUE 149 : IMM=IM 150 : 352 CONTINUE 151 : 353 CONTINUE 152 : LIGN=ILIGN(IL1) 153 : IL11=IPREL 154 : C 155 : C **** BOUCLE *5* TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE 156 : C 157 : * if (ngmaxz.eq.ngmaxy) write (6,*) 'resolution in core ',il1,il2 158 : * if (ngmaxz.ne.ngmaxy) write (6,*) 'resolution out of core ', 159 : * > il1,il2 160 : lig1=ilign(imin) 161 : ipos=0 162 : iper=imin 163 : ider=imin-1 164 : ngmaxz=ngmaxy 165 : DO 5 I=IMIN,IL2 166 : IMASQ=ITMASQ(I) 167 : LIG1=ILIGN(I) 168 : IF(I.LT.IL1) GO TO 7 169 : C 170 : C ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A 171 : C ******* IPREL IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER 172 : C ******* LE TERME DIAGONAL 173 : C 174 : LIGN=LIG1 175 : DO 156 KHG=1,IMMM(/1) 176 : II=IPREL-1+KHG 177 : IMMM(KHG)=0 178 : NN=IPPVV(KHG+1) 179 : NNM1=IPPVV(KHG) 180 : N=NN-NNM1 181 : DIAG(II)=VAL(NN) 182 : IF(N.EQ.1) GO TO 8 183 : NMI=N-II 184 : IDEP=MAX(IL11,2-NMI) 185 : KIDEP=IDEP+NMI 186 : KI1=N-1 187 : KQ=-NMI 188 : VAL(NN)=VAL(NN)+ 189 : # CHOLE1(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI), 190 : # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG), 191 : # PREC) 192 : if (imasq(nn/masdim+1).le.0) imasq(nn/masdim+1)=1 193 : 8 CONTINUE 194 : IF(ITTR(II).EQ.0.AND. 195 : & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.D-10) GO TO 12 196 : IF(ITTR(II).NE.0.AND. 197 : & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.D-6 ) GO TO 12 198 : * il faut mettre une valeur plus grande sur les LX car on a un probleme de conditionnement 199 : * sur le calcul des reactions en cas de 2 relations presque identique 200 : C 201 : C **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE 202 : C **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT 203 : C **** AU PREALABLE SUR CETTE INCONNUE. 204 : C 205 : * write (6,*) ' chole mode d ensemble ittr ligne ', 206 : * > ittr(ii),ii,diag(ii),val(nn) 207 : * on garde le signe car il fau un moins sur les ML 208 : if (ittr(ii).eq.0) then 209 : VAL(NN)=DIAG(II)*1d0 210 : else 211 : VAL(NN)=val(nn)-4.D0*DIAG(II) 212 : NENSLX=NENSLX+1 213 : endif 214 : NENS=NENS+1 215 : IMMM(KHG)=NENS 216 : 12 CONTINUE 217 : DIAG(II)=VAL(NN) 218 : IF(DIAG(II).NE.0.d0) GO TO 41 219 : KQ1=1+NNM1 220 : KQN=N+NNM1 221 : DO 16 LFG=KQ1,KQN 222 : if (imasq(lfg/masdim+1).le.0) goto 16 223 : IF(VAL(LFG).NE.0.d0) GO TO 17 224 : 16 CONTINUE 225 : DIAG(II)=1.d0 226 : if (ittr(ii).ne.0) diag(ii)=-1.d0 227 : VAL(NN)=DIAG(II) 228 : GO TO 41 229 : 17 CONTINUE 230 : * write (6,*) ' chole apres 17 ',val(lfg) 231 : diag(ii)=val(lfg) 232 : if (ittr(ii).ne.0) diag(ii)=-abs(diag(ii)) 233 : val(nn)=diag(ii) 234 : goto 41 235 : C 236 : C **** ENVOI ERREUR MATRICE SINGUIERE 237 : C 238 : ITYP=49 239 : INTERR(1)=I 239 : CALL ERREUR(49) 239 : if (nbthr.gt.1) call threadis 240 : RETURN 241 : C 242 : C **** ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS 243 : C ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE 244 : C 245 : 41 IF(DIAG(II).LT.0.D0) INEG=INEG+1 246 : IF(ITTR(II).NE.0) NBLAG=NBLAG+1 247 : 156 CONTINUE 248 : NA=IMMM(/1) 249 : * 250 : * RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE) 251 : * 251 : CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1), 252 : # NVALL,IPPVV(1),IZROSF,NA,PREC,imasq(1)) 253 : * 254 : IF (NBPAR.NE.NA+1.OR.NVALL.NE.VAL(/1)) SEGADJ LIGN 255 : NVSTOR=NVSTOR+NVALL 256 : nvstrm=max(nvstrm,nvall) 257 : DO 143 LHG=1,NBPAR 258 : IVPO(2*LHG-1)=KIVPO(LHG) 259 : IVPO(2*LHG)=KIVLO(LHG) 260 : 143 CONTINUE 261 : segsup imasq 262 : imasq=itmasq(i) 263 : if (i.gt.1) then 264 : lig1=ilign(i-1) 265 : segdes lig1 266 : endif 267 : c 268 : C **** ON TRIANGULARISE LES AUTRE LIGNES 269 : C 270 : il1=il1+1 271 : if (il1.gt.il2) goto 5 272 : LIG1=ILIGN(I) 273 : lign=ilign(il1) 274 : IL11=IPREL 275 : goto 7 276 : 71 continue 277 : if (i.ne.il1-1) then 278 : write (6,*) ' passage au grand ngmax ',i,il1,il2 279 : ngmaxz=2 000 000 000 280 : endif 281 : * soit parce qu'on a fini, soit parce qu'on manque de memoire 282 : * il faut executer les lignes activees puis les desactiver 283 : * lancer les chole3 et attendre qu'ils soient finis 284 : if (ipos.ne.0) then 285 : * write (6,*) ' lancement thread ',iper,ider,il1,il2 286 : if (iper.gt.ider) then 287 : write (6,*) ' erreur interne chole ' 287 : call erreur(5) 288 : endif 289 : if (nbthr.ge.1) then 290 : do ith=2,nbthr 290 : call threadid(ith,chole3i) 291 : enddo 291 : call chole3i(1) 292 : do ith=2,nbthr 292 : call threadif(ith) 293 : enddo 294 : * en multithread il peut y avoir n'importe quoi dans oov(1) du 295 : * aux acces simultan�s et ca crache gemat. donc : 296 : oov(1)=0 297 : else 298 : do ith=1,nbthr 298 : call chole3i(ith) 299 : enddo 300 : endif 301 : endif 302 : iposm=max(iposm,ipos) 303 : ipos=0 304 : iderf=ider-1 305 : if (ider.ne.il1-1) iderf=ider 306 : do il=iderf,iper,-1 307 : lign=ilign(il) 308 : segdes lign 309 : * write (6,*) ' desactivation ligne ',il 310 : enddo 311 : iper=ider+1 312 : * write (6,*) ' iper ider il1 ',iper,ider,il1 313 : if (iper.ne.il1) goto 7 314 : goto 5 315 : 7 CONTINUE 316 : SEGACT/err=71/LIG1 317 : ipos=ipos+1 318 : ider=i 319 : if (i.eq.il1-1) goto 71 320 : 5 CONTINUE 321 : * write (6,*) ' il1 il2 apres 5 ',il1,il2 322 : DO 11 I=IL1,IL2 323 : LIGN=ILIGN(I) 324 : SEGDES,LIGN 325 : 11 CONTINUE 326 : IF(IL2.LT.INO) GO TO 1 327 : C ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF 328 : C ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE 329 : c INEG=INEG-NBLAG 330 : * on ne compte pas 2 fois les multiplicateurs qui vont etre 331 : * eliminés lors de la resolution car mode d'ensemble 332 : INEG=INEG-(NBLAG-NENSLX) 333 : if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX 334 : 4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS', 335 : 1' DE LAGRANGE DETECTES') 336 : 337 : if (iimpi.eq.1) WRITE(IOIMP,4821) NVSTOC 338 : 4821 FORMAT( ' NOMBRE DE VALEURS DANS LE PROFIL',I12) 339 : if (iimpi.eq.1) WRITE(IOIMP,4822) NVSTOR 340 : 4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12) 341 : if (iimpi.eq.1) WRITE(IOIMP,4823) NVaor 342 : 4823 FORMAT( ' NOMBRE DE VALEURS initiales',I9) 343 : INTERR(1)=NVSTOR 344 : reaerr(1)=nvstor/inc**(4./3) 344 : if (ipasv.eq.0) CALL ERREUR(-278) 345 : IPASV=1 346 : SEGDES,MINCPO 347 : SEGDES,MIMIK 348 : SEGDES,MMATRI 349 : SEGDES,MILIGN 350 : SEGDES,MDIAG 351 : MMATRX=MMATRI 352 : SEGSUP KIVPO,KIVLO,ITMASQ 353 : * write (6,*) ' chole ipos max ',iposm 353 : if (nbthr.gt.1) call threadis 354 : RETURN 355 : END 356 : 357 : 358 : 359 : 360 : 361 : 362 : 363 : 364 : 365 : 366 : 367 : 368 : 369 : 370 : 371 : 372 : 373 : 374 : 375 : 376 : 377 : 378 : 379 : 380 : 381 : 382 : 383 : 384 : 385 : 386 : 387 : 388 : 389 : 390 : 391 : 392 : 393 : 394 : 395 : 396 : 397 :
© Cast3M 2003 - All rights reserved.
Disclaimer