Sources Esope | Cast3M

ellp58
C ELLP58    SOURCE    CHAT      05/01/12    23:35:55     5004
      SUBROUTINE ELLP58(ZA,ZB,ZU,N,N2,NNPOI,ICHAR,NIT,
     *                  A,B,U,R,R1,P,P1,CH,EPS)
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Y)
      IMPLICIT COMPLEX*16(Z)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C    OPERATEUR ...
C
C    RESOLUTION D'UN SYSTEME LINEAIRE COMPLEXE   ZA1 * ZX = ZV
C    PAR LA METHODE ITERATIVE DU GRADIANT CONJUGUE
C
C    PARAMETRES :
C
C    ZA(N*N)  : MATRICE (N*N) REPRESENTANT LE SYSTEME LINEAIRE
C    ZB(N)    : SECOND MEMBRE
C    ZU(N)    : TABLEAU SOLUTION
C
C    INTERMEDIAIRES :
C
C    A(2N*2N) :
C    B(2N)    :--->  MATRICES REELLES DU SYSTEME REEL CORRESPONDANT
C    U(2N)    :
C
C    AUTRES MATRICES  :  MATRICES DE TRAVAIL
C
C    SORTIE :
C
C    ZX       : SOLUTION
C
C    AUTEUR : SAINT - DIZIER
C    DATE   : 16 FEV 1990
C
C  ====================================================================
C
      COMPLEX*16  ZA(N,*),ZB(*),ZU(*)
C     COMPLEX*16  ZR(*),ZDELTA(*),ZP(*),ZH(*),ZH1(*)
C
      REAL*8 A(N2,*)
      REAL*8 B(*)
      REAL*8 U(*)
      REAL*8 R(*)
      REAL*8 R1(*)
      REAL*8 P(*)
      REAL*8 P1(*)
      REAL*8 CH(*)
C
      ZI = (0.D0 , -1.D0)
      N2 = N*2
C
      DO 1 I = 1 , N
          B(I  ) = ZB(I)
          B(I+N) = ZB(I)*ZI
          U(I  ) = ZU(I)
          U(I+N) = ZU(I)*ZI
          DO 2 J = 1 , N
              A(I  ,J  ) = ZA(I,J)
              A(I+N,J+N) = A(I,J)
              A(I+N,J  ) = ZA(I,J)*ZI
              A(I  ,J+N) = -A(I+N,J)
2         CONTINUE
1     CONTINUE
C
C
C     EPS = 1.D-8
C
C
C  ...............................  R = B - A*U
C
        DO 10 I = 1 , N2
            R(I) = B(I)
            DO 11 J = 1 , N2
                R(I) = R(I) - A(I,J)*U(J)
11          CONTINUE
10      CONTINUE
C
C  ...............................  P1 = P = R1 = R
C
        DO 12 I = 1 , N2
            R1(I) = R(I)
            P (I) = R(I)
            P1(I) = R(I)
12      CONTINUE
C
        NIT = 0
100     NIT = NIT + 1
C  ...............................  CH = (A) * P
C
        DO 20 I = 1 , N2
            CH (I) = 0.D0
            DO 21 J = 1 , N2
              CH (I) = CH(I) + A(I,J)*P(J)
21          CONTINUE
20      CONTINUE
C
C  ...............................  DEL1 = (R)T * R1
C
        DEL1 = 0.D0
        DO 30 I = 1 , N2
            DEL1 = DEL1 + R(I)*R1(I)
30      CONTINUE
C
C  ...............................  ALFA = DEL1/(CH)T*P1
C
        D = 0.D0
        DO 40 I = 1 , N2
            D = D + CH(I)*P1(I)
40      CONTINUE
C
        ALFA = DEL1 / D
C
C  ...............................  ZU = ZU + ZLANDA*ZP
C  ...............................  ZDELTA = ZDELTA + ZLANDA*ZH1
C  ...............................  DELTA1 = (ZDELTA)T * ZDELTA
C
        DEL2 = 0.D0
        DO 80 I = 1 , N2
            U(I) = U(I) + ALFA * P(I)
            R(I) = R(I) - ALFA * CH(I)
80      CONTINUE
C
        DO 81 I = 1 , N2
            CH (I) = 0.D0
            DO 82 J = 1 , N2
              CH (I) = CH(I) + A(J,I)*P1(J)
82          CONTINUE
C
            R1(I) = R1(I) - ALFA*CH(I)
            DEL2 = DEL2 + R(I)*R1(I)
C
81      CONTINUE
C
        BETA = DEL2 / DEL1
C
        XNU = 0.D0
        XDU = 0.D0
        XNR = 0.D0
        XDR = 0.D0
        XNF = 0.D0
        XDF = 0.D0
        XNM = 0.D0
        XDM = 0.D0
        ND12 = N / 12
        DO 85 I = 1 , ND12
            II = (I-1)*12
            DO 86 J = 1 , 3
              XNU = XNU + ALFA*ALFA*(P(II+J  )**2 + P(II+J  +N)**2)
              XNR = XNR + ALFA*ALFA*(P(II+J+3)**2 + P(II+J+3+N)**2)
              XNF = XNF + ALFA*ALFA*(P(II+J+6)**2 + P(II+J+6+N)**2)
              XNM = XNM + ALFA*ALFA*(P(II+J+9)**2 + P(II+J+9+N)**2)
C
              XDU = XDU + U(II+J  )**2 + U(II+J  +N)**2
              XDR = XDR + U(II+J+3)**2 + U(II+J+3+N)**2
              XDF = XDF + U(II+J+6)**2 + U(II+J+6+N)**2
              XDM = XDM + U(II+J+9)**2 + U(II+J+9+N)**2
86          CONTINUE
85      CONTINUE
C
C       NN = (NNPOI-1)*12 + ICHAR
C       XNUM = P(NN)**2 +P(NN+N)**2
C       XDEN = U(NN)**2 +U(NN+N)**2
C
        XST1 = SQRT(XNU/XDU)
        XST2 = SQRT(XNR/XDR)
        XST3 = SQRT(XNF/XDF)
        XST4 = SQRT(XNM/XDM)
        XSTOP = SQRT(XST1*XST1+XST2*XST2+XST3*XST3+XST4*XST4)
C       WRITE(6,*)'XST1 :',XST1
C       WRITE(6,*)'XST2 :',XST2
C       WRITE(6,*)'XST3 :',XST3
C       WRITE(6,*)'XST4 :',XST4
C
        IF (XSTOP.LE.EPS) THEN
          DO 88 I = 1 , N
              ZU(I) = U(I) - ZI*U(I+N)
88        CONTINUE
          WRITE(6,*) 'VECTEUR SOLUTION :',NNPOI,ICHAR
C         DO 89 I = 1 , N
C             WRITE(6,*) I , ZU(I)
C9        CONTINUE
C
          RETURN
        END IF
C
C
C ************************************************************
C
        DO 90 I = 1 , N2
            P (I) = BETA*P (I) + R (I)
            P1(I) = BETA*P1(I) + R1(I)
90      CONTINUE
C
        GO TO 100
C
      END