Télécharger cgresd.eso

Retour à la liste

cgresd
C CGRESD    SOURCE    CHAT      05/01/12    21:53:31     5004
      SUBROUTINE CGRESD(KSTO,
     1                  NL,ILG,
     2                  IMAT,IPRC,IA,JA,KA,
     3                  IZB,IZP,ICOLD,
     4                  NPT,NPITE,NEFF,ICONV,EPI,IPOU,VPOU,
     5                  NIMPR,IPAT)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                            C
C          RESOLUTION ITERATIVE D'UN SYSTEME                 C
C          SYMETRIQUE CREUX PAR LA METHODE DU                C
C          GRADIENT CONJUGUE PRECONDITIONNE PAR LA           C
C          DIAGONALE.                                        C
C                                                            C
C***** INFOS SUR LA METHODE (LIGNE 0)                        C
C                                                            C
C    DEUX MODES DE STOCKAGE SONT PROPOSES                    C
C    MORSE OU COMPRESSE.                                     C
C                                                            C
C          MORSE       (KSTO=0)                              C
C          COMPRESSE   (KSTO=1)                              C
C                                                            C
C    POUR LE CAS DCG, ON RESOUD L'EQUATION SUIVANTE          C
C                                                            C
C          ( (K-1)*A*(K-1) )*xi = (K-1)*b                    C
C                                                            C
C          AVEC K = D**(1/2)                                 C
C               pi est la pseudo solution                    C
C               La solution x est obtenue en faisant         C
C               x = (K-1)*xi                                 C
C                                                            C
C***** INFOS SUR LES DIMENSIONS (LIGNE 1)                    C
C                                                            C
C          NL (L'ORDRE DU SYSTEME)                           C
C                                                            C
C          ILG (LONGUEUR LA POUR LE MORSE)                   C
C              (NOMBRE MAXI DE VALEURS PAR LIGNE DE LA       C
C               MATRICE POUR LE COMPRESSE)                   C
C                                                            C
C***** INFOS SUR LA MATRICE (LIGNE 2)                        C
C                                                            C
C          IMAT (LA MATRICE DANS UN LISTREEL)                C
C                                                            C
C          IPRC (LE PRECONDITIONNEMENT DANS UN LISTREEL)     C
C                                                            C
C               DCG : IPRC EST LE POINTEUR SUR UN LISTREEL   C
C                     CONTENANT D**(-1/2) CE LISTREEL A POUR C
C                     DIMENSION NL AU MOINS                  C
C                                                            C
C          IA   (TABLEAU DE POINTEURS DU STOCKAGE MORSE)     C
C                                                            C
C               IL DOIT ETRE DIMENSIONNE A NL+1 AU MINIMUM   C
C                                                            C
C          JA   (TABLEAU DE CONNECTIVITES DE LA MATRICE)     C
C                                                            C
C               IL EST STOCKE EN MORSE ET DE LONGUEUR LA     C
C                                                            C
C          KA   (TABLEAU DE CONNECTIVITES EN MODE COMPRESSE) C
C                                                            C
C               IL EST DE LA FORME KA(NL,NNZ) AVEC NNZ       C
C               LE NOMBRE MAXI DE TERMES NON NULS POUR UNE   C
C               LIGNE DE LA MATRICE                          C
C                                                            C
C  NB : EN MORSE, KA N'EST PAS UTILISE                       C
C       EN COMPRESSE, IA ET JA NE SONT PAS UTILISES          C
C                                                            C
C  IA,JA ET KA SONT DES LISTENTI.                            C
C                                                            C
C***** INFOS SUR LE SECOND MEMBRE, LA SOLUTION INITIALE      C
C***** ET LA SOLUTION FINALE.   (LIGNE 3)                    C
C                                                            C
C          IZB    SECOND MEMBRE EN ENTREE                    C
C                 PAS MODIFIE.                               C
C                                                            C
C          IZP    ESTIMATION DE LA SOLUTION EN ENTREE        C
C                 LA SOLUTION EN SORTIE, IL FAUT LA          C
C                 CONSERVER POUR LA RESOLUTION SUIVANTE      C
C                                                            C
C  IZB ET IZP SONT DES LISTREEL DIMENSIONNES A NL            C
C                                                            C
C          ICOLD =1 (DEMARRAGE FROID)                        C
C          ICOLD!=1 (DEMARRAGE CHAUD A PARTIR DE LA          C
C                    SOLUTION INITIALE FOURNIE DANS IZP)     C
C                                                            C
C***** INFOS SUR LA CONVERGENCE (LIGNE 4)                    C
C                                                            C
C          NPT    NOMBRE MAXI D'ITERATIONS EN ENTREE         C
C                 SI IL EST < 10, IL EST MIS A 10            C
C          NPITE  FREQUENCE DES TESTS DE CONVERGENCE         C
C          NEFF   NOMBRE REEL EN SORTIE                      C
C          ICONV  0 NON CONVERGE, 1 CONVERGE                 C
C          EPI    PRECISION DEMANDEE                         C
C          IPOU=1 LA PRECISION VA ETRE POUSSE D'UN FACTEUR   C
C                 VPOU. CECI PEUT S'AVERER UTILE LORS D'UN   C
C                 CALCUL ITERATIF OU ON VEUT DES LE DEBUT    C
C                 UNE SOLUTION PRECISE POUR EVITER DE PARTIR C
C                 SUR UNE MAUVAISE PISTE. VPOU A UNE VALEUR  C
C                 PAR DEFAUT DE 100.                         C
C                                                            C
C***** INFOS SUR LE CONTROLE DES IMPRESSIONS (LIGNE 5)       C
C                                                            C
C          NIMPR =1 (INFOS SUR LA CONVERGENCE)               C
C          NIMPR!=1 (SILENCIEUX SAUF EN CAS DE DIVERGENCE)   C
C                                                            C
C          IPAT  =  (FACULTATIF, PAR EXEMPLE UN PAS DE TEMPS)C
C                                                            C
C RSETD(B,A)          B <- A                                C
C ADIVEC(A,B,C,alpha)  C <- A + alpha*B                      C
C PIMDV(A,B,C)         C <- A*B                              C
C                                                            C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC SMLENTI
      POINTEUR IA.MLENTI
      POINTEUR JA.MLENTI
      POINTEUR KA.MLENTI
-INC SMLREEL
      POINTEUR IMAT.MLREEL
      POINTEUR IPRC.MLREEL
      POINTEUR IZB.MLREEL
      POINTEUR IZBB.MLREEL
      POINTEUR IZP.MLREEL
      POINTEUR IZPS.MLREEL
      POINTEUR IZZPR.MLREEL
      POINTEUR IZPMV.MLREEL
      POINTEUR IZTMP.MLREEL
      POINTEUR IZDIR.MLREEL
 
      DATA VPOU0/100./
 
      IF(KSTO.EQ.0) LA=ILG
      IF(KSTO.EQ.1) NNZ=ILG
      IF(IPOU.EQ.1) THEN
         VPOUE=VPOU0
         IF(VPOU.NE.VPOU0) VPOUE=VPOU
      ENDIF
 
      CALL OOOMRU(1)
 
C On sauve le second membre en le dupliquant dans IZBB
C qui lui sera ecrase.
      SEGACT IZB
      JG=IZB.PROG(/1)
      SEGINI IZBB
      CALL RSETD(IZBB.PROG,IZB.PROG,JG)
      SEGDES IZB
 
      SEGACT IZP
 
      JG=NL
      SEGINI IZPS
      SEGINI IZZPR
      SEGINI IZPMV
      SEGINI IZTMP
      SEGINI IZDIR
 
      SEGACT IMAT
      SEGACT IPRC
      IF(KSTO.EQ.0) SEGACT IA
      IF(KSTO.EQ.0) SEGACT JA
      IF(KSTO.EQ.1) SEGACT KA
 
      IF(NPT.LT.10) NPT=10
      IF(NPT.GT.NL) NPT=NL
 
C Mettre le second membre et l'estimation initiale
C a l'echelle, B <- D-(1/2)*B, PIM <- PIM/(D**(-1/2))
      CALL PIMDV(IZBB.PROG,IPRC.PROG,IZBB.PROG,NL)
      CALL DIVISE(NL,IZP.PROG,IZP.PROG,IPRC.PROG)
      IF(ICOLD.EQ.1) CALL INITD (IZP.PROG,NL,0.D0)
      CALL RSETD(IZPS.PROG,IZP.PROG,NL)
      CALL PIMDV(IZP.PROG,IPRC.PROG,IZTMP.PROG,NL)
      IF(KSTO.EQ.0) THEN
C      CALL PMVM(A,IA,JA,X,Y,NL,LA)
       CALL PMVM(IMAT.PROG,IA.LECT,JA.LECT,
     1           IZTMP.PROG,IZPMV.PROG,NL,LA)
      ELSE
C      CALL PMVC(A,KA,X,Y,NL,NNZ)
       CALL PMVC(IMAT.PROG,KA.LECT,
     1           IZTMP.PROG,IZPMV.PROG,NL,NNZ)
      ENDIF
      CALL PIMDV(IZPMV.PROG,IPRC.PROG,IZPMV.PROG,NL)
      CALL ADIVEC (IZBB.PROG,IZPMV.PROG,IZBB.PROG,-1.d0,NL)
      ANORB=SQRT(PISCAL(IZBB.PROG,IZBB.PROG,NL))
      CALL RSETD(IZZPR.PROG,IZBB.PROG,NL)
      CALL RSETD(IZDIR.PROG,IZZPR.PROG,NL)
      SCAR1=PISCAL(IZBB.PROG,IZDIR.PROG,NL)
C
C===========DEBUT DE LA BOUCLE ITERATIVE=====600========
C
      NEFF=0
      NCOMP=0
      ICONV=0
      DO 600 I=1,NPT
       NEFF=NEFF+1
       NCOMP=NCOMP+1
       CALL PIMDV(IZDIR.PROG,IPRC.PROG,IZTMP.PROG,NL)
       IF(KSTO.EQ.0) THEN
C       CALL PMVM(A,IA,JA,X,Y,NL,LA)
        CALL PMVM(IMAT.PROG,IA.LECT,JA.LECT,
     1            IZTMP.PROG,IZPMV.PROG,NL,LA)
       ELSE
C       CALL PMVC(A,KA,X,Y,NL,NNZ)
        CALL PMVC(IMAT.PROG,KA.LECT,
     1            IZTMP.PROG,IZPMV.PROG,NL,NNZ)
       ENDIF
       CALL PIMDV(IZPMV.PROG,IPRC.PROG,IZPMV.PROG,NL)
       SCADIR=PISCAL(IZDIR.PROG,IZPMV.PROG,NL)
       ALPHA=SCAR1/SCADIR
       CALL ADIVEC(IZP.PROG,IZDIR.PROG,IZP.PROG,ALPHA,NL)
       ALP=-1.D0*ALPHA
       CALL ADIVEC (IZBB.PROG,IZPMV.PROG,IZBB.PROG,ALP,NL)
       CALL RSETD(IZZPR.PROG,IZBB.PROG,NL)
       SCAR2=PISCAL(IZBB.PROG,IZZPR.PROG,NL)
       BETA=SCAR2/SCAR1
       CALL ADIVEC (IZZPR.PROG,IZDIR.PROG,IZDIR.PROG,BETA,NL)
       IF (NCOMP.EQ.NPITE.OR.I.EQ.NPT) THEN
        INDIC=0
        IF (I.EQ.NPT) INDIC=1
C On pousse la precision si necessaire
        ZEPS=EPI
        IF(IPOU.EQ.1) ZEPS=EPI/VPOUE
        CALL COTEST (IZP.PROG,IZPS.PROG,
     1               ZEPS,DELTA,DELTAP,
     2               INDIC,ICONV,NL)
        SEPA=DELTA
        RSEPA=DELTA/DELTAP
        CALL RSETD(IZPS.PROG,IZP.PROG,NL)
        NCOMP=0
       ENDIF
       IF (ICONV.EQ.1) GOTO 601
       SCAR1=SCAR2
 600  CONTINUE
 601  CONTINUE
      CALL MTABL(IZBB.PROG,RES,IMA,NL)
      RRES=RES/ANORB
C Pour le DCG la solution est D**(-1/2)*IZP
C Pour le DCG on remet dans PIM, la solution
      CALL PIMDV(IZP.PROG,IPRC.PROG,IZP.PROG,NL)
      IF (ICONV.NE.1) THEN
       WRITE(6,1111)IPAT,NPT,RSEPA,DELTAP,RES,IZBB.PROG(IMA)
      ELSE
       IF(NIMPR.EQ.1) WRITE(6,2222)IPAT,EPI,
     1                             RES,SEPA,I,
     2                             RRES,RSEPA,IMA
      ENDIF
 
      SEGDES IZP
 
      SEGSUP IZBB
      SEGSUP IZPS
      SEGSUP IZZPR
      SEGSUP IZDIR
      SEGSUP IZTMP
      SEGSUP IZPMV
 
      SEGDES IMAT*(NOMOD,MRU)
      SEGDES IPRC*(NOMOD,MRU)
      IF(KSTO.EQ.0) SEGDES IA*(NOMOD,MRU)
      IF(KSTO.EQ.0) SEGDES JA*(NOMOD,MRU)
      IF(KSTO.EQ.1) SEGDES KA*(NOMOD,MRU)
 
      CALL OOOMRU(0)
 
C
 1111 FORMAT(1X,'*CGRESD* ',I6,' NON CONV EN ',I5,
     1' ITERA, RSEPA = ',D8.2,1X,'DELTAP = ',D8.2,1X,
     2' RES = ',D8.2,' P = ',D8.2)
 2222 FORMAT(1X,'*CGRESD* ',I6,' *PREC ',D8.2,
     1          ' * RESIDU-SEPARA @ NITE ',D8.2,1X,D8.2,1X,I5,
     2          ' * RRES-RSEP ',D8.2,1X,D8.2,
     3          ' * ELEMENT ',I7)
C
C  FIN RESOLUTION
C
      RETURN
      END
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales