C CHOMOD SOURCE PV 13/10/31 21:17:34 7850 C C DECOMPOSITION DE CROUT MODIFIEE POUR LE CALCUL DES MATRICES C CONDENSEES D'UN SUPER-ELEMENT C RESULTAT DANS XMATRI C SUBROUTINE CHOMOD(MMATRX,NBNNMA,SNTT,SNTO,XMATRI,NLIGRA) 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,ILIMITE SEGMENT ITMASQ(NBLIG) SEGMENT IMASQ(LMASQ) external chole3i C LOGICAL TEST1,LIMIT C C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX C -INC CCOPTIO -INC SMMATRI -INC SMRIGID -INC CCASSIS -INC CCHOLE C SEGMENT SNTO INTEGER NTOTMA(NN) ENDSEGMENT C SEGMENT SNTT INTEGER NTTMAI(NN) ENDSEGMENT C SEGMENT KIVPO(IIMAX) SEGMENT KIVLO(IIMAX) SAVE IPASV DATA IPASV/0/ DATA PREC /1D-18/ * character*8 zen * equivalence (zen,izen) * multithread iposm=0 * zen='CPU'//char(0) * le=4 * nbthro=0 nvaor=0 * call ooozen(izen,le) * if (le.eq.1) read(zen,fmt='(i1)') nbthro * if (le.eq.2) read(zen,fmt='(i2)') nbthro * if (nbthro.eq.1) nbthro=0 * nbthro=max(nbesc,2*nbthro,1) * nbthro=min(nbthro,64) nbthro=max( nbesc,1) if (nbthro.gt.1) call threadii nbthr=nbthro * nouvelle methode de gestion de l'espace memoire necessitee par la para * memoire vive totale MACTIT=OOOVAL(1,1) * un bloc de memoire fera au plus macti/2 nvstrm=0 MMATRI=MMATRX SEGACT,MMATRI*MOD PRCHLV=PREC MILIGN=IILIGN SEGACT,MILIGN*MOD INO=ILIGN(/1) MDIAG=IDIAG SEGACT,MDIAG*MOD NBLIG=INO SEGINI ITMASQ precc=prec INC=DIAG(/1) nvstrm=inc INCC=INC MIMIK=IIMIK MINCPO=IINCPO SEGACT,MINCPO,MIMIK IPLUMI=IMIK(/2)*2 +4 IL2=0 IIMAX=IJMAX+IPLUMI SEGINI KIVPO,KIVLO NLIGRP=NLIGRA NLIGRD=NLIGRA NGMAXZ=NGMAXY nelrig=1 SEGINI,XMATRI C *** ILIMITE = nb de noeuds qui ne sont pas maitres C SEGACT,SNTO SEGACT,SNTT ILIMITE = INO - (NTOTMA(/1) + NTTMAI(/1)) TEST1=.FALSE. IF(ILIMITE.EQ.0) THEN IL2=0 TEST1=.TRUE. ENDIF INEG=0 NBLAG=0 NENSLX=0 NVSTOC=0 NVSTOR=0 LIMIT=.FALSE. C C **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD, C **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUD C C **** LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX C 1 CONTINUE IF (TEST1) LIMIT=.TRUE. C IVALMA=IJMAX+IPLUMI IL1=IL2+1 IMIN=IL1 mactic=0 DO 2 I=IL1,INO VALREG=1D50 LLIGN=ILIGN(I) SEGACT /ERR=3/LLIGN NA=IMMMM(/1) NBPAR=NA+1 NVALL=NJMAX * IF(IVALMA+NVALL.GT.NGMAXZ.AND.I.GT.IL1 ) GO TO 3 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 nvaor=nvaor+xxva(/1) C C **** DECOMPACTAGE C IPA=1 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) 122 CONTINUE * * REDUCTION EVENTUELLE DE LA LIGNE A SA PARTIE SIGNIFICATIVE * > 1E-15*DIAGONALE * VALREF=ABS(VAL(IPA+IMMMM(JPA)-LPA))*prec VALREG=MIN(VALREG,VALREF) ITROP=0 IPAHAU=IPA+IMMMM(JPA)-LPA DO 123 ILD=IPA,IPAHAU-1 IF (ABS(VAL(ILD)).GT.VALREF) GOTO 124 ITROP=ITROP+1 123 CONTINUE 124 CONTINUE 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 (NVALL.NE.VAL(/1)) SEGADJ LIGN LMASQ=NVALL/MASDIM+1 SEGINI IMASQ DO ILD=1,VAL(/1) IF (ABS(VAL(ILD)).GT.VALREG) IMASQ(ILD/MASDIM+1)=1 ENDDO ITMASQ(I)=IMASQ **** **** **** IPREL=IMMMM(1) IDERL=IMMMM(NA) IPPVV(NA+1)=IPA-1 SEGSUP 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. > i-il1+1.ge.2*nbths.and.mod(i-il1+1,nbths).eq.0 > .and.mactic.gt.NGMAXZ*nbths) then il2=i GOTO 4 endif C IF ((I.EQ.(ILIMITE)).AND.(TEST1.EQV..FALSE.)) THEN TEST1=.TRUE. IL2 = I GOTO 4 ENDIF 2 CONTINUE IL2=INO GO TO 4 3 IL2=I-1 4 CONTINUE C IF(IL2.GE.IL1) GO TO 40 C C **** APPEL AUX ERREURS MESSAGE PAS ASSEZ DE PLACE MEMOIRE C ITYP=48 CALL ERREUR(48) if (nbthr.gt.1) call threadis RETURN 40 CONTINUE 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 lig1=ilign(imin) ipos=0 iper=imin ider=imin-1 ngmaxz=ngmaxy DO 5 I=IMIN,IL2 IMASQ=ITMASQ(I) LIG1=ILIGN(I) IF(I.LT.IL1) GO TO 7 C____________ 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) THEN IF(LIMIT) THEN RE(II-NBNNMA,II-NBNNMA,1)=VAL(NN) GO TO 41 ELSE GO TO 8 ENDIF ENDIF NMI=N-II KI1=N-1 KQ=-NMI IF(LIMIT) THEN C C ****** CAS NOEUD MAITRE C IDEP=MAX(IL11,1-NMI) KIDEP=IDEP+NMI VAL(NN)=VAL(NN)+ # SUPER1(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI), # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,NBNNMA,XMATRI, # IMASQ(1),1+IPPVV(KHG)) RE(II-NBNNMA,II-NBNNMA,1)=VAL(NN) GOTO 41 ENDIF C C ****** CAS NOEUD NON MAITRE C IDEP=MAX(IL11,2-NMI) KIDEP=IDEP+NMI VAL(NN)=VAL(NN)+ # CHOLE1(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) C 8 CONTINUE 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 * write (6,*) ' chomod nens ittr val diag ',ittr(ii),nens+1, * > val(nn),diag(ii),ii,ipno(ii),ilimite if (ittr(ii).eq.0) then VAL(NN)=DIAG(II)*1d0 else VAL(NN)=val(nn)-4.D0*DIAG(II) NENSLX=NENSLX+1 endif NENS=NENS+1 IMMM(KHG)=NENS 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(VAL(LFG).NE.0.d0) GO TO 17 16 CONTINUE DIAG(II)=1.d0 if (ittr(ii).ne.0) diag(ii)=-1.d0 VAL(NN)=diag(ii) GO TO 41 17 CONTINUE diag(ii)=val(lfg) if (ittr(ii).ne.0) diag(ii)=-abs(diag(ii)) val(nn)=diag(ii) goto 41 C C **** ENVOI ERREUR MATRICE SINGULIERE C ITYP=49 INTERR(1)=I CALL ERREUR(49) if (nbthr.gt.1) call threadis 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) * CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1), # NVALL,IPPVV(1),IZROSF,NA,PREC,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 C 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 72 continue * write (6,*) ' chomod segact ',i,ILIMITE,il1,il2 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 ' call erreur(5) endif IF(I.LE.ILIMITE) THEN if (nbthr.gt.1) then do ith=1,min(il2-il1+1,nbthr-1) call threadid(ith,chole3i) enddo call chole3i(nbthr) do ith=1,min(il2-il1+1,nbthr-1) 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 chole3i(ith) enddo endif endif iposm=max(iposm,ipos) 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 (i.eq.ilimite) goto 5 if (iper.ne.il1) goto 7 ENDIF goto 5 7 CONTINUE SEGACT/err=72/LIG1 ipos=ipos+1 ider=i if (i.eq.il1-1.and.i.le.ilimite) goto 71 IF(I.EQ.ILIMITE) GOTO 71 IF(I.GT.ILIMITE) THEN IPPR=LIG1.IPREL IDDR=LIG1.IDERL C C DO 10 JBI=IL1,IL2 LIGN=ILIGN(JBI) IMASQ=ITMASQ(JBI) C C ****** CAS NOEUD MAITRE C CALL SUPER3(IPREL,IDERL,IPPVV(1),IPPR,IDDR, # IVPO(1),LIG1.IPPVV(1),VAL(1), # LIG1.VAL(1),LIG1.IVPO(1),NBNNMA,XMATRI, # imasq(1)) IF(IMM.GT.I) GO TO 106 10 CONTINUE 106 CONTINUE IF(I.LT.IL1) then SEGDES,LIG1 endif ENDIF 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 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 4821 FORMAT( ' CHOMOD NOMBRE DE VALEURS DANS LE PROFIL',I9) IF (IIMPI.EQ.1)WRITE(IOIMP,4822) NVSTOR 4822 FORMAT( ' CHOMOD NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I9) IF (IIMPI.EQ.1)WRITE(IOIMP,4823) NVaor 4823 FORMAT( ' CHOMOD NOMBRE DE VALEURS initiales',I9) INTERR(1)=NVSTOR reaerr(1)=nvstor/inc**(4./3) if (ipasv.eq.0) CALL ERREUR(-278) IPASV=1 SEGDES,MINCPO SEGDES,MIMIK SEGDES,MMATRI SEGDES,MILIGN SEGDES,MDIAG SEGDES,SNTO SEGDES,SNTT SEGDES,XMATRI MMATRX=MMATRI SEGSUP KIVPO,KIVLO,ITMASQ if (nbthr.gt.1) call threadis RETURN END