C CHOLE SOURCE PV 13/10/31 21:17:33 7850 SUBROUTINE CHOLE(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 C C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX C -INC CCOPTIO -INC SMMATRI -INC CCASSIS -INC CCHOLE SEGMENT KIVPO(IIMAX) SEGMENT KIVLO(IIMAX) SEGMENT ITMASQ(NBLIG) SEGMENT IMASQ(LMASQ) external chole3i SAVE IPASV DATA IPASV/0/ * character*8 zen * equivalence (zen,izen) iposm=0 * zen='CPU'//char(0) * le=4 * nbthro=0 nvaor=0 nbthro=max(nbesc,1) if (nbthro.gt.1) call threadii nbthr=nbthro * 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/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 INEG=0 NBLAG=0 NENSLX=0 NVSTOC=0 NVSTOR=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 IL1=IL2+1 * write (6,*) ' il1 il2 apres 1 ',il1,il2 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 nvaor = nvaor + xxva(/1) C C **** DECOMPACTAGE C LMASQ=NVALL/MASDIM+2 SEGINI IMASQ 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) IMASQ(IPLA/MASDIM+1)=1 122 CONTINUE IPA=IPA+IMMMM(JPA)-LPA+1 IMMM(JPA)=IPNO(LPA) IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA) 121 CONTINUE **** **** **** ITMASQ(I)=IMASQ IPREL=IMMMM(1) IDERL=IMMMM(NA) IPPVV(NA+1)=IPA-1 SEGSUP LLIGN ILIGN(I)=LIGN nbths=nbthr * nb de ligne multiple du nb de threads if (mod(i+1-il1,nbths).eq.0 > .and.mactic.gt.NGMAXZ*nbths) then * write (6,*) ' chole il1 il2 ',il1,i 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) 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 ******* 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)+ # 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) if (imasq(nn/masdim+1).le.0) imasq(nn/masdim+1)=1 8 CONTINUE IF(ITTR(II).EQ.0.AND. & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.D-10) GO TO 12 IF(ITTR(II).NE.0.AND. & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.D-6 ) GO TO 12 * il faut mettre une valeur plus grande sur les LX car on a un probleme de conditionnement * sur le calcul des reactions en cas de 2 relations presque identique 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,*) ' chole mode d ensemble ittr ligne ', * > ittr(ii),ii,diag(ii),val(nn) * on garde le signe car il fau un moins sur les ML 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 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)=-1.d0 VAL(NN)=DIAG(II) GO TO 41 17 CONTINUE * write (6,*) ' chole apres 17 ',val(lfg) 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 SINGUIERE 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 AUTRE LIGNES C il1=il1+1 if (il1.gt.il2) goto 5 LIG1=ILIGN(I) lign=ilign(il1) IL11=IPREL goto 7 71 continue if (i.ne.il1-1) then write (6,*) ' passage au grand ngmax ',i,il1,il2 ngmaxz=2 000 000 000 endif * 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 (nbthr.ge.1) then do ith=2,nbthr call threadid(ith,chole3i) enddo call chole3i(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 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 (iper.ne.il1) goto 7 goto 5 7 CONTINUE SEGACT/err=71/LIG1 ipos=ipos+1 ider=i if (i.eq.il1-1) goto 71 5 CONTINUE * write (6,*) ' il1 il2 apres 5 ',il1,il2 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( ' 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) if (iimpi.eq.1) WRITE(IOIMP,4823) NVaor 4823 FORMAT( ' 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 MMATRX=MMATRI SEGSUP KIVPO,KIVLO,ITMASQ * write (6,*) ' chole ipos max ',iposm if (nbthr.gt.1) call threadis RETURN END