C LDMT3 SOURCE PV 13/11/12 21:15:00 7864 SUBROUTINE LDMT3(MMATRX,PREC) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) * TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY PARAMETER (LPCRAY=10000) INTEGER OOOVAL POINTEUR LLL.LIGN, MMM.LIGN POINTEUR LIG2.LIGN POINTEUR L.LLIGN, M.LLIGN POINTEUR LL.MILIGN, MM.MILIGN SEGMENT ITEMP REAL*8 P(INC) ENDSEGMENT POINTEUR R.ITEMP,W.ITEMP C C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX C -INC CCOPTIO -INC SMMATRI -INC CCHOLE SEGMENT KIVPO(IIMAX) SEGMENT KIVLO(IIMAX) SEGMENT ITMASQ(NBLIG) SEGMENT ITMASI(NBLIGI) SEGMENT IMASQ(LMASQ) SEGMENT IMASQI(LMASQI) SAVE IPASV DATA IPASV/0/ * prec=1.d-20 * nouvelle methode de gestion de l'espace memoire necessitee par la parallelisation * memoire vive totale MACTIT=OOOVAL(1,1) * un bloc de memoire fera au plus macti/3 nvstrm=0 MMATRI=MMATRX SEGACT,MMATRI*MOD PRCHLV=PREC LL=IILIGN MM=IILIGS SEGACT,MM*MOD,LL*MOD INO=MM.ILIGN(/1) MDIAG=IDIAG SEGACT,MDIAG*MOD NBLIG=INO NBLIGI=INO SEGINI ITMASQ,ITMASI INC=DIAG(/1) INCC=INC MIMIK=IIMIK MINCPO=IINCPO SEGACT,MINCPO,MIMIK IPLUMI=IMIK(/2)*2 +4 IL2=0 IIMAX=IJMAX+IPLUMI SEGINI KIVPO,KIVLO INEG=0 NBLAG=0 NENSLX=0 NVSTOC=0 NVSTOR=0 NVSTIC=0 NVSTIR=0 NGMAXZ=NGMAXY 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 IVALMI=IJMAX+IPLUMI IL1=IL2+1 IMIN=IL1 IMINII=IL1 mactic=0 DO 2 I=IL1,INO VALREG=1D50 L=LL.ILIGN(I) SEGACT /ERR=3/L M=MM.ILIGN(I) SEGACT /ERR=3/M NA=M.IMMMM(/1) NBPAR=NA+1 NVALL=M.NJMAX NVALLL=NVALL IF(IVALMA+NVALL.GT.NGMAXZ.AND.I.GT.IL1 ) GO TO 3 SEGINI /ERR=3/MMM NVALL= L.NJMAX NVILL=NVALL NA=L.IMMMM(/1) NBPAR=NA+1 IF(IVALMI+NVALL.GT.NGMAXZ.AND.I.GT.IL1 ) GO TO 3 SEGINI /ERR=3/LLL * recuperer la longueur du segment mactic=mactic+OOA(OOA(OOT+MMM)+OOT-1)+OOA(OOA(OOT+LLL)+OOT-1) if ((mactic.gt.mactit/3.or.mactic.gt.mactit-nvstrm*4) > .AND.I.GT.IL1) THEN SEGSUP LLL,MMM GOTO 3 endif NVSTOC=NVSTOC + NVALLL NVSTIC=NVSTIC+NVALL IVALMA=IVALMA+NVALLL IVALMI=IVALMI+NVALL NVALL=NVALLL C C **** DECOMPACTAGE C IPA=1 IPI=1 DO 121 JPA=1,NA MMM.IVPO(JPA)=IPA LLL.IVPO(JPA)=IPI KPA=M.IPPO(JPA+1)-M.IPPO(JPA) KPI=L.IPPO(JPA+1)-L.IPPO(JPA) IPP=M.IPPO(JPA) IPPI= L.IPPO(JPA) MMM.IPPVV(JPA)=IPA-1 LLL.IPPVV(JPA)=IPI-1 LPA=M.LDEB(JPA) LPAI=L.LDEB(JPA) LPA1=LPA-IPA LPA1I=LPAI-IPI DO 122 MPA=1,KPA LLO=M.LINC(MPA+IPP) IPLA=LLO-LPA1 MMM.VAL(IPLA)=M.XXVA(MPA+IPP) 122 CONTINUE DO 1220 MPA=1,KPI LLO=L.LINC(MPA+IPPI) IPLA=LLO-LPA1I LLL.VAL(IPLA)=L.XXVA(MPA+IPPI) 1220 CONTINUE * REDUCTION EVENTUELLE DE LA LIGNE A SA PARTIE SIGNIFICATIVE * > 1E-15*DIAGONALE **** *** VALREF=ABS(MMM.VAL(IPA+M.IMMMM(JPA)-LPA))*PREC VALREG=MIN(VALREG,VALREF) ITROP=0 ITRIP=0 IPAHAU=IPA+M.IMMMM(JPA)-LPA IPIHAU=IPI+L.IMMMM(JPA)-LPAI ** DO 123 ILD=IPA,IPAHAU-1 ** IF (ABS(MMM.VAL(ILD)).GT.VALREF) GOTO 124 ** ITROP=ITROP+1 *123 CONTINUE *124 CONTINUE ** DO 1230 ILD=IPI,IPIHAU-1 ** IF (ABS(LLL.VAL(ILD)).GT.VALREF) GOTO 1240 ** ITRIP=ITRIP+1 *230 CONTINUE *240 CONTINUE ** IF (ITROP.NE.0) THEN ** LPA=LPA+ITROP ** JILD=IPAHAU-ITROP ** DO 125 ILD=IPA,JILD ** MMM.VAL(ILD)=MMM.VAL(ILD+ITROP) *125 CONTINUE ** DO 126 ILD=JILD+1,JILD+ITROP ** MMM.VAL(ILD)=0.D0 *126 CONTINUE ** NVALL=NVALL-ITROP ** NVSTOC=NVSTOC-ITROP ** IVALMA=IVALMA-ITROP ** ENDIF ** IF (ITRIP.NE.0) THEN ** LPAI=LPAI+ITRIP ** JILD=IPIHAU-ITRIP ** DO 1250 ILD=IPI,JILD ** LLL.VAL(ILD)=LLL.VAL(ILD+ITRIP) 1250 CONTINUE ** DO 1260 ILD=JILD+1,JILD+ITRIP ** LLL.VAL(ILD)=0.D0 1260 CONTINUE ** NVILL=NVILL-ITRIP ** NVSTIC=NVSTIC-ITRIP ** IVALMI=IVALMI-ITRIP ** ENDIF IPA=IPA+M.IMMMM(JPA)-LPA+1 MMM.IMMM(JPA)=MM.IPNO(LPA) IF(IMIN.GT.MM.IPNO(LPA)) IMIN=MM.IPNO(LPA) IPI=IPI+L.IMMMM(JPA)-LPAI+1 LLL.IMMM(JPA)=LL.IPNO(LPAI) IF(IMINII.GT.LL.IPNO(LPAI)) IMINII=LL.IPNO(LPAI) 121 CONTINUE IF (NVALL.NE.MMM.VAL(/1))SEGADJ MMM NVALLL=NVALL NVALL=NVILL IF (NVILL.NE.LLL.VAL(/1))SEGADJ LLL NVALL=NVALLL **** **** **** LMASQ=NVALL/MASDIM+2 LMASQI=NVILL/MASDIM+2 SEGINI IMASQ,IMASQI DO ILD=1,MMM.VAL(/1) IF (ABS(MMM.VAL(ILD)).GT.VALREG) IMASQ(ILD/MASDIM+1)=1 ENDDO DO ILD=1,LLL.VAL(/1) IF (ABS(LLL.VAL(ILD)).GT.VALREG) IMASQI(ILD/MASDIM+1)=1 ENDDO ITMASQ(I)=IMASQ ITMASI(I)=IMASQI MMM.IPREL=M.IMMMM(1) MMM.IDERL=M.IMMMM(NA) LLL.IPREL=L.IMMMM(1) LLL.IDERL=L.IMMMM(NA) MMM.IPPVV(NA+1)=IPA-1 LLL.IPPVV(NA+1)=IPI-1 SEGSUP L,M MM.ILIGN(I)=MMM LL.ILIGN(I)=LLL * IF(I.GE.1000000) THEN * WRITE( IOIMP,4987) I * 4987 FORMAT (' NOEUD NUMERO ',I5 , ' colonne puis ligne') * LLO=MMM.VAL(/1) * write(6,*) 'M IMMM ', (MMM.IMMM(IU),IU=1,MMM.IMMM(/1) ) * WRITE(IOIMP, 4926)(MMM.VAL(MPA),MPA=1,LLO) *4926 FORMAT(' M VAL ' , 6E11.4) * write(6,7632) ( MMM.IPPVV(IU),iu=1,MMM.IPPVV(/1)) *7632 FORMAT(' M IPPVV',10i6) * LLO=LLL.VAL(/1) * write(6,*) 'L IMMM ', (LLL.IMMM(IU),IU=1,LLL.IMMM(/1) ) * WRITE(IOIMP, 4927)(LLL.VAL(MPA),MPA=1,LLO) *4927 FORMAT(' L VAL ' , 6E11.4) * write(6,7633) ( LLL.IPPVV(IU),iu=1,LLL.IPPVV(/1)) *7633 FORMAT(' L IPPVV',10i6) * 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 IB=IL2+1 MMM=MM.ILIGN(IL2) MMM.IMM=IB IL2M=IL2-1 IF(IL2M.LT.IL1) GO TO 353 DO 352 IH=IL2M,IL1,-1 MMM=MM.ILIGN(IH+1) IM=MMM.IMMM(1) DO 354 JH=1,MMM.IMMM(/1) IM=MIN(IM,MMM.IMMM(JH)) 354 CONTINUE LIG1=MM.ILIGN(IH) LIG1.IMM=MIN(MMM.IMM,IM) 352 CONTINUE 353 CONTINUE C MMM=ILIGN(IL1) IB=IL2+1 LLL=LL.ILIGN(IL2) LLL.IMM=IB IL2M=IL2-1 IF(IL2M.LT.IL1) GO TO 3530 DO 3520 IH=IL2M,IL1,-1 LLL=LL.ILIGN(IH+1) IM=LLL.IMMM(1) DO 3540 JH=1,LLL.IMMM(/1) IM=MIN(IM,LLL.IMMM(JH)) 3540 CONTINUE LIG2=LL.ILIGN(IH) LIG2.IMM=MIN(LLL.IMM,IM) 3520 CONTINUE 3530 CONTINUE MMM=MM.ILIGN(IL1) IL11=MMM.IPREL LLL=LL.ILIGN(IL1) IL22=LLL.IPREL C C **** BOUCLE *5* TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE C lig1=MM.ilign(imin) call oooeta(lig1,ieta) lig2=LL.ilign(iminII) call oooeta(lig2,ietai) if (ieta.eq.3.OR.IETAI.EQ.3) then * on cherche a reprendre des lignes migres sur disque * on passe dans le mode ou ngmaxz est grand ! ngmaxz=2000000000 endif IMINA=MIN(IMIN,IMINII) DO 5 I=IMINA,IL2 IMASQ=ITMASQ(I) IMASQI=ITMASI(I) LIG1=MM.ILIGN(I) LIG2=LL.ILIGN(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 C on s'occupe d'abord de la partie superieur LIGN=LIG1 NAA= IMMM(/1) DO 156 KHG=1,NAA LIGN=LIG1 II=IPREL-1+KHG IMMM(KHG)=0 NN=IPPVV(KHG+1) NNS=NN NNM1=IPPVV(KHG) NNM1S=NNM1 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)+ #CHOLI1(LL.ILIGN,LIG2,VAL(1+IPPVV(KHG)),DIAG(1-NMI),LL.IPNO(1-NMI), # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG), # PREC,1) 8 CONTINUE LIGN=LIG2 IMMM(KHG)=0 NN=IPPVV(KHG+1) NNM1=IPPVV(KHG) N=NN-NNM1 IF(N.EQ.1) GO TO 88 NMI=N-II IDEP=MAX(IL22,2-NMI) KIDEP=IDEP+NMI KI1=N-1 KQ=-NMI VAL(NN)=VAL(NN)+ #CHOLI1(MM.ILIGN,LIG1,VAL(1+IPPVV(KHG)),DIAG(1-NMI),MM.IPNO(1-NMI), # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQI(1),1+IPPVV(KHG), # PREC,2) 88 CONTINUE LIGN=LIG1 IF(LL.ITTR(II).EQ.0.AND. & ABS(LIG2.VAL(NN)).GT.ABS(DIAG(II))*1.E-10) GO TO 12 IF(LL.ITTR(II).NE.0.AND. & ABS(LIG2.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 LIG2.VAL(NN)=DIAG(II) NENS=NENS+1 LIG2.IMMM(KHG)=NENS if(LL.ittr(ii).NE.0) NENSLX=NENSLX+1 12 CONTINUE IMASQ(NNS/MASDIM+1)=1 IMASQI(NN/MASDIM+1)=1 DIAG(II)=LIG2.VAL(NN) IF(DIAG(II).NE.0.) GO TO 41 KQ1=1+NNM1S KQN=N+NNM1S DO 16 LFG=KQ1,KQN IF(LIG1.VAL(LFG).NE.0.) GO TO 17 16 CONTINUE KQ1=1+NNM1 KQN=N+NNM1 DO 160 LFG=KQ1,KQN IF(LIG2.VAL(LFG).NE.0.) GO TO 17 160 CONTINUE DIAG(II)=1. LIG2.VAL(NN)=1. 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(LL.ITTR(II).NE.0) NBLAG=NBLAG+1 LIG1.VAL(NNS)=LIG2.VAL(NN) 156 CONTINUE NA=IMMM(/1) * * RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE) * d'abord la triangulaire superieure * CALL COMPAC(LIG1.VAL(1),NBPAR,KIVPO(1),KIVLO(1), # NVALL,LIG1.IPPVV(1),IZROSF,NA,PREC,imasq(1)) IF (NBPAR.NE.NA+1.OR.NVALL.NE.LIG1.VAL(/1)) SEGADJ LIG1 NVSTOR=NVSTOR+NVALL nvstrm=max(nvstrm,nvall) DO 143 LHG=1,NBPAR LIG1.IVPO(2*LHG-1)=KIVPO(LHG) LIG1.IVPO(2*LHG)=KIVLO(LHG) 143 CONTINUE * IF(I.EQ.910) THEN * WRITE (6,*) ' NOEU ',I,' DIAG INVERSE ',LIG1.VAL(VAL(/1)) * WRITE (6,*) ' LIGNE APRES COMPACTION VAL ' * WRITE (6,*) (LIG1.VAL(IK),IK=1,LIG1.VAL(/1)) * WRITE (6,*) 'KIVPO ' * WRITE (6,*) (KIVPO(IK),IK=1,NBPAR) * WRITE (6,*) 'KIVLO ' * WRITE (6,*) (KIVLO(IK),IK=1,NBPAR) * write(6,*) ( LIG1.IPPVV(IU),IU=1,LIG1.IPPVV(/1)) * ENDIF * * RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE) * puis la triangulaire superieure * NBPAR=NA+1 CALL COMPAC(LIG2.VAL(1),NBPAR,KIVPO(1),KIVLO(1), # NVILL,LIG2.IPPVV(1),IZROSF,NA,PREC,imasqi(1)) NVALL=NVILL IF (NBPAR.NE.NA+1.OR.NVILL.NE.LIG2.VAL(/1)) SEGADJ LIG2 NVSTIR=NVSTIR+NVILL nvstrm=max(nvstrm,nvIll) DO 1430 LHG=1,NBPAR LIG2.IVPO(2*LHG-1)=KIVPO(LHG) LIG2.IVPO(2*LHG)=KIVLO(LHG) 1430 CONTINUE segsup imasq imasq=itmasq(i) * IF(I.EQ.910) THEN * WRITE (6,*) ' NOEU ',I,' DIAG INVERSE ',LIG2.VAL(VAL(/1)) * WRITE (6,*) ' LIGNE APRES COMPACTION VAL ' * WRITE (6,*) (LIG2.VAL(IK),IK=1,LIG2.VAL(/1)) * WRITE (6,*) 'KIVPO ' * WRITE (6,*) (KIVPO(IK),IK=1,NBPAR) * WRITE (6,*) 'KIVLO ' * WRITE (6,*) (KIVLO(IK),IK=1,NBPAR) * write(6,*) ( LIG2.IPPVV(IU),IU=1,LIG2.IPPVV(/1)) * ENDIF C **** ON TRIANGULARISE LES AUTRES LIGNES C GOTO 5 7 CONTINUE IF(I.GE.IMIN) THEN SEGACT,LIG2 IPPR=LIG2.IPREL IDDR=LIG2.IDERL C C **** BOUCLE *10* SUR LES LIGNES DES NOEUDS INFERIEUR JUSQU'A C **** IL2 ON CALCULE UNE COLONNE C DO 10 JBI=IL1,IL2 LIGN=MM.ILIGN(JBI) IMASQ=ITMASQ(JBI) CALL CHOLE3(IPREL,IDERL,IPPVV(1),IPPR,IDDR,IVPO(1), # LIG2.IPPVV(1),VAL(1),LIG2.VAL(1),LIG2.IVPO(1), > imasq(1),prec) IF(IMM.GT.I) GO TO 106 10 CONTINUE 106 CONTINUE * IF(I.LT.IL1) SEGDES,LIG1 SEGDES,LIG2 ENDIF IF(I.GE.IMINII) THEN SEGACT,LIG1 IPPR=LIG1.IPREL IDDR=LIG1.IDERL C C **** BOUCLE *10* SUR LES LIGNES DES NOEUDS INFERIEUR JUSQU'A C **** IL2 ON CALCULE UNE COLONNE C DO 100 JBI=IL1,IL2 LIGN=LL.ILIGN(JBI) IMASQ=ITMASI(JBI) CALL CHOLE3(IPREL,IDERL,IPPVV(1),IPPR,IDDR,IVPO(1), # LIG1.IPPVV(1),VAL(1),LIG1.VAL(1),LIG1.IVPO(1), > imasq(1),prec) IF(IMM.GT.I) GO TO 1060 100 CONTINUE 1060 CONTINUE * IF(I.LT.IL1) SEGDES,LIG2 SEGDES,LIG1 ENDIF 5 CONTINUE DO 11 I=IL1,IL2 LIGN=LL.ILIGN(I) SEGDES,LIGN LIGN=MM.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 c 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+NVSTIC 4821 FORMAT( ' NOMBRE DE VALEURS DANS LE PROFIL',I9) IF (IIMPI.EQ.1)WRITE(IOIMP,4822) NVSTOR+NVSTIR 4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I9) INTERR(1)=NVSTOR+NVSTIR IF (IPASV.EQ.0) CALL ERREUR(-278) IPASV=1 SEGDES,MINCPO SEGDES,MIMIK SEGDES,MMATRI SEGDES,LL,MM SEGDES,MDIAG MMATRX=MMATRI SEGSUP KIVPO,KIVLO,ITMASQ RETURN END