C GRACO2 SOURCE PV 13/10/31 21:17:35 7850 SUBROUTINE GRACO2(MMATRX) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) * TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY * * COPIE DE CHOLE POUR FAIRE UN CHOLEVSKI INCOMPLET DE PLUS ON GARDE LA * MATRICE ASSEMBLEE ( SEMGENT LLIGN) * PARAMETER (LPCRAY=10000) INTEGER OOOVAL C C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX C SEULES LES VALEURS INITIALEMENT NON NULLES SONT CALCULEES C -INC CCOPTIO -INC SMMATRI -INC CCASSIS -INC CCHOLE SEGMENT KIVPO(IIMAX) SEGMENT KIVLO(IIMAX) SEGMENT ITMASQ(NBLIG) SEGMENT IMASQ(LMASQ) external graco5i nbthro=max(nbesc,1) nbthr=nbthro ** en fait la parallelisation n'apporte rien car la matrice est trop creuse ** nbthr=1 if (nbthr.gt.1) call threadii MACTIT=OOOVAL(1,1) prec=1D-30 nvstrm=0 MMATRI=MMATRX SEGACT,MMATRI*MOD MILIG1=IILIGN IASLIG= IILIGN SEGINI,MILIGN=MILIG1 IILIGN = MILIGN INO=ILIGN(/1) MDIA1=IDIAG IASDIA=IDIAG SEGINI,MDIAG=MDIA1 C------------------------ IDIAG = MDIAG C------------------------- NBLIG=INO SEGINI ITMASQ ** comme la matrice est beaucoup plus creuse qu'en direct, on prend un masdim plus petit precc=prec INC=DIAG(/1) INCC=INC MIMIK=IIMIK MINCPO=IINCPO SEGACT,MINCPO,MIMIK IPLUMI=IMIK(/2)*2 +4 IL2=0 LTRK=MAX(LPCRAY+1,OOOVAL(1,4)) IIMAX=((((IJMAX+IPLUMI) +LTRK)/LTRK)+1)*LTRK SEGINI KIVPO,KIVLO INEG=0 NBLAG=0 NENSLX=0 NVSTOC=0 NVSTOR=0 C C **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD, C **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS C C **** LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX C 1 CONTINUE IVALMA=IJMAX+IPLUMI IL1=IL2+1 IMIN=IL1 mactic=0 DO 2 I=IL1,INO LLIGN=ILIGN(I) SEGACT /ERR=3/LLIGN NA=IMMMM(/1) NBPAR=NA+1 NVALL=NJMAX SEGINI /ERR=3/ LIGN * recuperer la longueur du segment mactic=mactic+OOA(OOA(OOT+lign)+OOT-1) if ((mactic.gt.mactit/2.or.mactic.gt.mactit-nvstrm*12) > .AND.I.GT.IL1) THEN SEGSUP LIGN GOTO 3 endif NVSTOC=NVSTOC + NVALL IVALMA=IVALMA+NVALL C C **** DECOMPACTAGE C LMASQ=NVALL/MASDIM+2 SEGINI IMASQ IPA=1 valmi=1D50 valma=1D-50 DO 121 JPA=1,NA IVPO(JPA)=IPA KPA=IPPO(JPA+1)-IPPO(JPA) IPP=IPPO(JPA) IPPVV(JPA)=IPA-1 LPA=LDEB(JPA) LPA1=LPA-IPA DO 122 MPA=1,KPA LL=LINC(MPA+IPP) IPLA=LL-LPA1 VAL(IPLA)=XXVA(MPA+IPP) IMASQ(IPLA/MASDIM+1)=1 valma=max(valma,abs(val(ipla))) if (val(ipla).ne.0.D0) valmi=min(valmi,abs(val(ipla))) 122 CONTINUE * REDUCTION EVENTUELLE DE LA LIGNE A SA PARTIE SIGNIFICATIVE * > 1E-50*DIAGONALE **** *** VALREF=ABS(VAL(IPA+IMMMM(JPA)-LPA))*1D-50 ITROP=0 IPAHAU=IPA+IMMMM(JPA)-LPA ** DO 123 ILD=IPA,IPAHAU-1 ** IF (ABS(VAL(ILD)).GT.VALREF) GOTO 124 ** ITROP=ITROP+1 **23 CONTINUE *124 CONTINUE ** ITROP=0 ** IF (ITROP.NE.0) THEN ** LPA=LPA+ITROP ** JILD=IPAHAU-ITROP ** DO 125 ILD=IPA,JILD ** VAL(ILD)=VAL(ILD+ITROP) *125 CONTINUE ** DO 126 ILD=JILD+1,JILD+ITROP ** VAL(ILD)=0.D0 *126 CONTINUE ** NVALL=NVALL-ITROP ** NVSTOC=NVSTOC-ITROP ** IVALMA=IVALMA-ITROP ** ENDIF IPA=IPA+IMMMM(JPA)-LPA+1 IMMM(JPA)=IPNO(LPA) IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA) 121 CONTINUE * if (valmi.le.1d-20) then * write (6,*) ' valmi valma diag ',valmi,valma,val(val(/1)) * endif IF (NVALL.NE.VAL(/1))SEGADJ LIGN **** **** **** ** LMASQ=NVALL/MASDIM+1 ** SEGINI IMASQ ** DO ILD=1,VAL(/1) ** IF (ABS(VAL(ILD)).GT.VALREF) IMASQ(ILD/MASDIM+1)=1 ** ENDDO ITMASQ(I)=IMASQ IPREL=IMMMM(1) IDERL=IMMMM(NA) IPPVV(NA+1)=IPA-1 * SEGSUP LLIGN * SEGDES LLIGN ILIGN(I)=LIGN IF(IIMPI.EQ.1525) THEN WRITE( IOIMP,4987) I 4987 FORMAT (' NOEUD NUMERO ',I5) LL=VAL(/1) WRITE(IOIMP, 4926)(VAL(MPA),MPA=1,LL) 4926 FORMAT(' VAL ' , 10E11.4) ENDIF nbths=nbthr ** if (mactic/OOA(OOA(OOT+lign)+OOT-1).ge.(nbths*2)/3.and. if(i-il1+1.ge.nbths.and.mod(i+1-il1,nbths).eq.0 > .and.mactic.gt.NGMAXY*nbths) then il2=i GOTO 4 endif 2 CONTINUE IL2=INO GO TO 4 3 IL2=I-1 4 CONTINUE IF(IL2.GE.IL1) GO TO 40 C C **** APPEL AUX ERREURS MESSAGE PAS ASSEZ DE PLACE MEMOIRE C ITYP=48 CALL ERREUR(48) RETURN 40 CONTINUE lig1=ilign(imin) ipos=0 IM=IL2+1 DO 352 IH=IL2,IL1,-1 LIGN=ILIGN(IH) DO 354 JH=1,IMMM(/1) IM=MIN(IM,IMMM(JH)) 354 CONTINUE IMM=IM 352 CONTINUE 353 CONTINUE LIGN=ILIGN(IL1) IL11=IPREL C C **** BOUCLE *5* TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE C * if (ngmaxz.eq.ngmaxy) write (6,*) 'resolution in core ',il1,il2 * if (ngmaxz.ne.ngmaxy) write (6,*) 'resolution out of core ', * > il1,il2 ipos=0 iper=imin ider=imin-1 DO 5 I=IMIN,IL2 if (ierr.ne.0) return LIG1=ILIGN(I) IMASQ=ITMASQ(I) IF(I.LT.IL1) GO TO 7 C C ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A C ******* IPREL IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER C ******* LE TERME DIAGONAL C LIGN=LIG1 DO 156 KHG=1,IMMM(/1) II=IPREL-1+KHG IMMM(KHG)=0 NN=IPPVV(KHG+1) NNM1=IPPVV(KHG) N=NN-NNM1 DIAG(II)=VAL(NN) IF(N.EQ.1) GO TO 8 NMI=N-II IDEP=MAX(IL11,2-NMI) KIDEP=IDEP+NMI KI1=N-1 KQ=-NMI VAL(NN)=VAL(NN)+ # GRACO3(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI), # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG), # PREC,ittr(1)) if (imasq(nn/masdim+1).le.0) imasq(nn/masdim+1)=1 8 CONTINUE if (diag(ii).ne.0.d0) then if (val(nn).lt.abs(diag(ii))*0.5) val(nn)=diag(ii) endif IF(ITTR(II).EQ.0.AND. & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-10) GO TO 12 IF(ITTR(II).NE.0.AND. & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-6) GO TO 12 C C **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE C **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT C **** AU PREALABLE SUR CETTE INCONNUE. C VAL(NN)=DIAG(II)*1d-4 NENS=NENS+1 IMMM(KHG)=NENS if(ittr(ii).NE.0) NENSLX=NENSLX+1 12 CONTINUE IMASQ(NN/MASDIM+1)=1 DIAG(II)=VAL(NN) IF(DIAG(II).NE.0.d0) GO TO 41 KQ1=1+NNM1 KQN=N+NNM1 DO 16 LFG=KQ1,KQN if (imasq(lfg/masdim+1).le.0) goto 16 IF(VAL(LFG).NE.0.d0) GO TO 17 16 CONTINUE DIAG(II)=1.D0 if (ittr(ii).ne.0) diag(ii)=-1d0 VAL(NN)=DIAG(II) GO TO 41 17 CONTINUE C C **** ENVOI ERREUR MATRICE SINGUIERE C ITYP=49 INTERR(1)=I CALL ERREUR(49) RETURN C C **** ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS C ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE C 41 IF(DIAG(II).LT.0.D0) INEG=INEG+1 IF(ITTR(II).NE.0) NBLAG=NBLAG+1 156 CONTINUE NA=IMMM(/1) * * RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE) * izro=2 CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1), # NVALL,IPPVV(1),izro ,NA ,1D-30,imasq(1)) IF (NBPAR.NE.NA+1.OR.NVALL.NE.VAL(/1)) SEGADJ LIGN NVSTOR=NVSTOR+NVALL nvstrm=max(nvstrm,nvall) DO 143 LHG=1,NBPAR IVPO(2*LHG-1)=KIVPO(LHG) IVPO(2*LHG)=KIVLO(LHG) 143 CONTINUE segsup imasq imasq=itmasq(i) if (i.gt.1) then lig1=ilign(i-1) segdes lig1 endif ** WRITE (6,*) ' NOEU ',I,' DIAG INVERSE ',VAL(VAL(/1)) ** WRITE (6,*) ' LIGNE APRES COMPACTION VAL ' ** WRITE (6,*) (VAL(IK),IK=1,VAL(/1)) ** WRITE (6,*) 'KIVPO ' ** WRITE (6,*) (KIVPO(IK),IK=1,NBPAR) ** WRITE (6,*) 'KIVLO ' ** WRITE (6,*) (KIVLO(IK),IK=1,NBPAR) C **** ON TRIANGULARISE LES AUTRES LIGNES C il1=il1+1 if (il1.gt.il2) goto 5 LIG1=ILIGN(I) lign=ilign(il1) IL11=IPREL goto 7 GOTO 5 71 continue if (i.ne.il1-1) ngmaxz=2 000 000 000 * soit parce qu'on a fini, soit parce qu'on manque de memoire * il faut executer les lignes activees puis les desactiver * lancer les chole3 et attendre qu'ils soient finis if (ipos.ne.0) then * write (6,*) ' lancement thread ',iper,ider,il1,il2 if (iper.gt.ider) then write (6,*) ' erreur interne chole ' write (6,*) ' noeud ',i,' sur ',ino write (6,*) ' nombre de termes actuels ',nvstor call erreur(5) endif if (nbthr.gt.1) then do ith=2,nbthr call threadid(ith,graco5i) enddo call graco5i(1) do ith=2,nbthr call threadif(ith) enddo * en multithread il peut y avoir n'importe quoi dans oov(1) du * aux acces simultanés et ca crache gemat. donc : oov(1)=0 else do ith=1,nbthr call graco5i(ith) enddo endif endif ipos=0 iderf=ider-1 if (ider.ne.il1-1) iderf=ider do il=iderf,iper,-1 lign=ilign(il) segdes lign * write (6,*) ' desactivation ligne ',il enddo iper=ider+1 * write (6,*) ' iper ider il1 ',iper,ider,il1 if (iper.ne.il1) goto 7 goto 5 7 CONTINUE ipos=ipos+1 SEGACT/err=71/LIG1 ider=i if (i.eq.il1-1) goto 71 5 CONTINUE * DO 11 I=IL1,IL2 * LIGN=ILIGN(I) SEGDES,LIGN * 11 CONTINUE IF(IL2.LT.INO) GO TO 1 C ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF C ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE INEG=INEG-NBLAG * on ne compte pas 2 fois les multiplicateurs qui vont etre * eliminés lors de la resolution car mode d'ensemble INEG=INEG-(NBLAG-NENSLX) if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX 4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS', 1' DE LAGRANGE DETECTES') IF(IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC 4821 FORMAT( ' NOMBRE DE VALEURS DANS LE PROFIL',I12) IF(IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR 4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12) SEGDES,MINCPO SEGDES,MIMIK SEGDES,MMATRI SEGDES,MILIGN SEGDES,MDIAG MMATRX=MMATRI SEGSUP KIVPO,KIVLO,ITMASQ if (nbthr.gt.1) call threadis RETURN END