1 : C GRACO2 SOURCE PV 13/10/31 21:17:35 7850 2 : SUBROUTINE GRACO2(MMATRX) 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 : * 7 : * COPIE DE CHOLE POUR FAIRE UN CHOLEVSKI INCOMPLET DE PLUS ON GARDE LA 8 : * MATRICE ASSEMBLEE ( SEMGENT LLIGN) 9 : * 10 : PARAMETER (LPCRAY=10000) 11 : INTEGER OOOVAL 12 : C 13 : C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX 14 : C SEULES LES VALEURS INITIALEMENT NON NULLES SONT CALCULEES 15 : C 16 : -INC CCOPTIO 17 : -INC SMMATRI 18 : -INC CCASSIS 19 : -INC CCHOLE 20 : SEGMENT KIVPO(IIMAX) 21 : SEGMENT KIVLO(IIMAX) 22 : SEGMENT ITMASQ(NBLIG) 23 : SEGMENT IMASQ(LMASQ) 24 : external graco5i 25 : nbthro=max(nbesc,1) 26 : nbthr=nbthro 27 : ** en fait la parallelisation n'apporte rien car la matrice est trop creuse 28 : ** nbthr=1 28 : if (nbthr.gt.1) call threadii 29 : MACTIT=OOOVAL(1,1) 30 : prec=1D-30 31 : nvstrm=0 32 : MMATRI=MMATRX 33 : SEGACT,MMATRI*MOD 34 : MILIG1=IILIGN 35 : IASLIG= IILIGN 36 : SEGINI,MILIGN=MILIG1 37 : IILIGN = MILIGN 38 : INO=ILIGN(/1) 39 : MDIA1=IDIAG 40 : IASDIA=IDIAG 41 : SEGINI,MDIAG=MDIA1 42 : C------------------------ 43 : IDIAG = MDIAG 44 : C------------------------- 45 : NBLIG=INO 46 : SEGINI ITMASQ 47 : ** comme la matrice est beaucoup plus creuse qu'en direct, on prend un masdim plus petit 48 : precc=prec 49 : INC=DIAG(/1) 50 : INCC=INC 51 : MIMIK=IIMIK 52 : MINCPO=IINCPO 53 : SEGACT,MINCPO,MIMIK 54 : IPLUMI=IMIK(/2)*2 +4 55 : IL2=0 56 : LTRK=MAX(LPCRAY+1,OOOVAL(1,4)) 57 : IIMAX=((((IJMAX+IPLUMI) +LTRK)/LTRK)+1)*LTRK 58 : SEGINI KIVPO,KIVLO 59 : INEG=0 60 : NBLAG=0 61 : NENSLX=0 62 : NVSTOC=0 63 : NVSTOR=0 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 : IMIN=IL1 74 : mactic=0 75 : DO 2 I=IL1,INO 76 : LLIGN=ILIGN(I) 77 : SEGACT /ERR=3/LLIGN 78 : NA=IMMMM(/1) 79 : NBPAR=NA+1 80 : NVALL=NJMAX 81 : SEGINI /ERR=3/ LIGN 82 : * recuperer la longueur du segment 83 : mactic=mactic+OOA(OOA(OOT+lign)+OOT-1) 84 : if ((mactic.gt.mactit/2.or.mactic.gt.mactit-nvstrm*12) 85 : > .AND.I.GT.IL1) THEN 86 : SEGSUP LIGN 87 : GOTO 3 88 : endif 89 : NVSTOC=NVSTOC + NVALL 90 : IVALMA=IVALMA+NVALL 91 : C 92 : C **** DECOMPACTAGE 93 : C 94 : LMASQ=NVALL/MASDIM+2 95 : SEGINI IMASQ 96 : IPA=1 97 : valmi=1D50 98 : valma=1D-50 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 : valma=max(valma,abs(val(ipla))) 112 : if (val(ipla).ne.0.D0) valmi=min(valmi,abs(val(ipla))) 113 : 122 CONTINUE 114 : * REDUCTION EVENTUELLE DE LA LIGNE A SA PARTIE SIGNIFICATIVE 115 : * > 1E-50*DIAGONALE 116 : **** *** 117 : VALREF=ABS(VAL(IPA+IMMMM(JPA)-LPA))*1D-50 118 : ITROP=0 119 : IPAHAU=IPA+IMMMM(JPA)-LPA 120 : ** DO 123 ILD=IPA,IPAHAU-1 121 : ** IF (ABS(VAL(ILD)).GT.VALREF) GOTO 124 122 : ** ITROP=ITROP+1 123 : **23 CONTINUE 124 : *124 CONTINUE 125 : ** ITROP=0 126 : ** IF (ITROP.NE.0) THEN 127 : ** LPA=LPA+ITROP 128 : ** JILD=IPAHAU-ITROP 129 : ** DO 125 ILD=IPA,JILD 130 : ** VAL(ILD)=VAL(ILD+ITROP) 131 : *125 CONTINUE 132 : ** DO 126 ILD=JILD+1,JILD+ITROP 133 : ** VAL(ILD)=0.D0 134 : *126 CONTINUE 135 : ** NVALL=NVALL-ITROP 136 : ** NVSTOC=NVSTOC-ITROP 137 : ** IVALMA=IVALMA-ITROP 138 : ** ENDIF 139 : IPA=IPA+IMMMM(JPA)-LPA+1 140 : IMMM(JPA)=IPNO(LPA) 141 : IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA) 142 : 121 CONTINUE 143 : * if (valmi.le.1d-20) then 144 : * write (6,*) ' valmi valma diag ',valmi,valma,val(val(/1)) 145 : * endif 146 : IF (NVALL.NE.VAL(/1))SEGADJ LIGN 147 : **** **** **** 148 : ** LMASQ=NVALL/MASDIM+1 149 : ** SEGINI IMASQ 150 : ** DO ILD=1,VAL(/1) 151 : ** IF (ABS(VAL(ILD)).GT.VALREF) IMASQ(ILD/MASDIM+1)=1 152 : ** ENDDO 153 : ITMASQ(I)=IMASQ 154 : IPREL=IMMMM(1) 155 : IDERL=IMMMM(NA) 156 : IPPVV(NA+1)=IPA-1 157 : * SEGSUP LLIGN 158 : * SEGDES LLIGN 159 : ILIGN(I)=LIGN 160 : IF(IIMPI.EQ.1525) THEN 161 : WRITE( IOIMP,4987) I 162 : 4987 FORMAT (' NOEUD NUMERO ',I5) 163 : LL=VAL(/1) 164 : WRITE(IOIMP, 4926)(VAL(MPA),MPA=1,LL) 165 : 4926 FORMAT(' VAL ' , 10E11.4) 166 : ENDIF 167 : nbths=nbthr 168 : ** if (mactic/OOA(OOA(OOT+lign)+OOT-1).ge.(nbths*2)/3.and. 169 : if(i-il1+1.ge.nbths.and.mod(i+1-il1,nbths).eq.0 170 : > .and.mactic.gt.NGMAXY*nbths) then 171 : il2=i 172 : GOTO 4 173 : endif 174 : 2 CONTINUE 175 : IL2=INO 176 : GO TO 4 177 : 3 IL2=I-1 178 : 4 CONTINUE 179 : IF(IL2.GE.IL1) GO TO 40 180 : C 181 : C **** APPEL AUX ERREURS MESSAGE PAS ASSEZ DE PLACE MEMOIRE 182 : C 183 : ITYP=48 183 : CALL ERREUR(48) 184 : RETURN 185 : 40 CONTINUE 186 : lig1=ilign(imin) 187 : ipos=0 188 : IM=IL2+1 189 : DO 352 IH=IL2,IL1,-1 190 : LIGN=ILIGN(IH) 191 : DO 354 JH=1,IMMM(/1) 192 : IM=MIN(IM,IMMM(JH)) 193 : 354 CONTINUE 194 : IMM=IM 195 : 352 CONTINUE 196 : 353 CONTINUE 197 : LIGN=ILIGN(IL1) 198 : IL11=IPREL 199 : C 200 : C **** BOUCLE *5* TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE 201 : C 202 : * if (ngmaxz.eq.ngmaxy) write (6,*) 'resolution in core ',il1,il2 203 : * if (ngmaxz.ne.ngmaxy) write (6,*) 'resolution out of core ', 204 : * > il1,il2 205 : ipos=0 206 : iper=imin 207 : ider=imin-1 208 : DO 5 I=IMIN,IL2 209 : if (ierr.ne.0) return 210 : LIG1=ILIGN(I) 211 : IMASQ=ITMASQ(I) 212 : IF(I.LT.IL1) GO TO 7 213 : C 214 : C ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A 215 : C ******* IPREL IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER 216 : C ******* LE TERME DIAGONAL 217 : C 218 : LIGN=LIG1 219 : DO 156 KHG=1,IMMM(/1) 220 : II=IPREL-1+KHG 221 : IMMM(KHG)=0 222 : NN=IPPVV(KHG+1) 223 : NNM1=IPPVV(KHG) 224 : N=NN-NNM1 225 : DIAG(II)=VAL(NN) 226 : IF(N.EQ.1) GO TO 8 227 : NMI=N-II 228 : IDEP=MAX(IL11,2-NMI) 229 : KIDEP=IDEP+NMI 230 : KI1=N-1 231 : KQ=-NMI 232 : VAL(NN)=VAL(NN)+ 233 : # GRACO3(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI), 234 : # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG), 235 : # PREC,ittr(1)) 236 : if (imasq(nn/masdim+1).le.0) imasq(nn/masdim+1)=1 237 : 8 CONTINUE 238 : if (diag(ii).ne.0.d0) then 239 : if (val(nn).lt.abs(diag(ii))*0.5) val(nn)=diag(ii) 240 : endif 241 : 242 : IF(ITTR(II).EQ.0.AND. 243 : & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-10) GO TO 12 244 : IF(ITTR(II).NE.0.AND. 245 : & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-6) GO TO 12 246 : C 247 : C **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE 248 : C **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT 249 : C **** AU PREALABLE SUR CETTE INCONNUE. 250 : C 251 : VAL(NN)=DIAG(II)*1d-4 252 : NENS=NENS+1 253 : IMMM(KHG)=NENS 254 : if(ittr(ii).NE.0) NENSLX=NENSLX+1 255 : 12 CONTINUE 256 : IMASQ(NN/MASDIM+1)=1 257 : DIAG(II)=VAL(NN) 258 : IF(DIAG(II).NE.0.d0) GO TO 41 259 : KQ1=1+NNM1 260 : KQN=N+NNM1 261 : DO 16 LFG=KQ1,KQN 262 : if (imasq(lfg/masdim+1).le.0) goto 16 263 : IF(VAL(LFG).NE.0.d0) GO TO 17 264 : 16 CONTINUE 265 : DIAG(II)=1.D0 266 : if (ittr(ii).ne.0) diag(ii)=-1d0 267 : VAL(NN)=DIAG(II) 268 : GO TO 41 269 : 17 CONTINUE 270 : C 271 : C **** ENVOI ERREUR MATRICE SINGUIERE 272 : C 273 : ITYP=49 274 : INTERR(1)=I 274 : CALL ERREUR(49) 275 : RETURN 276 : C 277 : C **** ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS 278 : C ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE 279 : C 280 : 41 IF(DIAG(II).LT.0.D0) INEG=INEG+1 281 : IF(ITTR(II).NE.0) NBLAG=NBLAG+1 282 : 156 CONTINUE 283 : NA=IMMM(/1) 284 : * 285 : * RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE) 286 : * 287 : izro=2 287 : CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1), 288 : # NVALL,IPPVV(1),izro ,NA ,1D-30,imasq(1)) 289 : 290 : IF (NBPAR.NE.NA+1.OR.NVALL.NE.VAL(/1)) SEGADJ LIGN 291 : NVSTOR=NVSTOR+NVALL 292 : nvstrm=max(nvstrm,nvall) 293 : DO 143 LHG=1,NBPAR 294 : IVPO(2*LHG-1)=KIVPO(LHG) 295 : IVPO(2*LHG)=KIVLO(LHG) 296 : 143 CONTINUE 297 : segsup imasq 298 : imasq=itmasq(i) 299 : if (i.gt.1) then 300 : lig1=ilign(i-1) 301 : segdes lig1 302 : endif 303 : ** WRITE (6,*) ' NOEU ',I,' DIAG INVERSE ',VAL(VAL(/1)) 304 : ** WRITE (6,*) ' LIGNE APRES COMPACTION VAL ' 305 : ** WRITE (6,*) (VAL(IK),IK=1,VAL(/1)) 306 : ** WRITE (6,*) 'KIVPO ' 307 : ** WRITE (6,*) (KIVPO(IK),IK=1,NBPAR) 308 : ** WRITE (6,*) 'KIVLO ' 309 : ** WRITE (6,*) (KIVLO(IK),IK=1,NBPAR) 310 : 311 : C **** ON TRIANGULARISE LES AUTRES LIGNES 312 : C 313 : il1=il1+1 314 : if (il1.gt.il2) goto 5 315 : LIG1=ILIGN(I) 316 : lign=ilign(il1) 317 : IL11=IPREL 318 : goto 7 319 : GOTO 5 320 : 71 continue 321 : if (i.ne.il1-1) ngmaxz=2 000 000 000 322 : * soit parce qu'on a fini, soit parce qu'on manque de memoire 323 : * il faut executer les lignes activees puis les desactiver 324 : * lancer les chole3 et attendre qu'ils soient finis 325 : if (ipos.ne.0) then 326 : * write (6,*) ' lancement thread ',iper,ider,il1,il2 327 : if (iper.gt.ider) then 328 : write (6,*) ' erreur interne chole ' 329 : write (6,*) ' noeud ',i,' sur ',ino 330 : write (6,*) ' nombre de termes actuels ',nvstor 330 : call erreur(5) 331 : endif 332 : if (nbthr.gt.1) then 333 : do ith=2,nbthr 333 : call threadid(ith,graco5i) 334 : enddo 334 : call graco5i(1) 335 : do ith=2,nbthr 335 : call threadif(ith) 336 : enddo 337 : * en multithread il peut y avoir n'importe quoi dans oov(1) du 338 : * aux acces simultanés et ca crache gemat. donc : 339 : oov(1)=0 340 : else 341 : do ith=1,nbthr 341 : call graco5i(ith) 342 : enddo 343 : endif 344 : endif 345 : ipos=0 346 : iderf=ider-1 347 : if (ider.ne.il1-1) iderf=ider 348 : do il=iderf,iper,-1 349 : lign=ilign(il) 350 : segdes lign 351 : * write (6,*) ' desactivation ligne ',il 352 : enddo 353 : iper=ider+1 354 : * write (6,*) ' iper ider il1 ',iper,ider,il1 355 : if (iper.ne.il1) goto 7 356 : goto 5 357 : 7 CONTINUE 358 : ipos=ipos+1 359 : SEGACT/err=71/LIG1 360 : ider=i 361 : if (i.eq.il1-1) goto 71 362 : 5 CONTINUE 363 : * DO 11 I=IL1,IL2 364 : * LIGN=ILIGN(I) 365 : SEGDES,LIGN 366 : * 11 CONTINUE 367 : IF(IL2.LT.INO) GO TO 1 368 : C ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF 369 : C ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE 370 : INEG=INEG-NBLAG 371 : * on ne compte pas 2 fois les multiplicateurs qui vont etre 372 : * eliminés lors de la resolution car mode d'ensemble 373 : INEG=INEG-(NBLAG-NENSLX) 374 : if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX 375 : 4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS', 376 : 1' DE LAGRANGE DETECTES') 377 : 378 : IF(IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC 379 : 4821 FORMAT( ' NOMBRE DE VALEURS DANS LE PROFIL',I12) 380 : IF(IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR 381 : 4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12) 382 : SEGDES,MINCPO 383 : SEGDES,MIMIK 384 : SEGDES,MMATRI 385 : SEGDES,MILIGN 386 : SEGDES,MDIAG 387 : MMATRX=MMATRI 388 : SEGSUP KIVPO,KIVLO,ITMASQ 388 : if (nbthr.gt.1) call threadis 389 : RETURN 390 : END 391 : 392 : 393 : 394 : 395 : 396 : 397 : 398 : 399 : 400 : 401 : 402 : 403 : 404 : 405 : 406 : 407 : 408 : 409 :
© Cast3M 2003 - All rights reserved.
Disclaimer