1 : C ELLP58 SOURCE CHAT 05/01/12 23:35:55 5004 2 : SUBROUTINE ELLP58(ZA,ZB,ZU,N,N2,NNPOI,ICHAR,NIT, 3 : * A,B,U,R,R1,P,P1,CH,EPS) 4 : C 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8(A-H,O-Y) 7 : IMPLICIT COMPLEX*16(Z) 8 : C 9 : C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 10 : C 11 : C OPERATEUR ... 12 : C 13 : C RESOLUTION D'UN SYSTEME LINEAIRE COMPLEXE ZA1 * ZX = ZV 14 : C PAR LA METHODE ITERATIVE DU GRADIANT CONJUGUE 15 : C 16 : C PARAMETRES : 17 : C 18 : C ZA(N*N) : MATRICE (N*N) REPRESENTANT LE SYSTEME LINEAIRE 19 : C ZB(N) : SECOND MEMBRE 20 : C ZU(N) : TABLEAU SOLUTION 21 : C 22 : C INTERMEDIAIRES : 23 : C 24 : C A(2N*2N) : 25 : C B(2N) :---> MATRICES REELLES DU SYSTEME REEL CORRESPONDANT 26 : C U(2N) : 27 : C 28 : C AUTRES MATRICES : MATRICES DE TRAVAIL 29 : C 30 : C SORTIE : 31 : C 32 : C ZX : SOLUTION 33 : C 34 : C AUTEUR : SAINT - DIZIER 35 : C DATE : 16 FEV 1990 36 : C 37 : C ==================================================================== 38 : C 39 : COMPLEX*16 ZA(N,*),ZB(*),ZU(*) 40 : C COMPLEX*16 ZR(*),ZDELTA(*),ZP(*),ZH(*),ZH1(*) 41 : C 42 : REAL*8 A(N2,*) 43 : REAL*8 B(*) 44 : REAL*8 U(*) 45 : REAL*8 R(*) 46 : REAL*8 R1(*) 47 : REAL*8 P(*) 48 : REAL*8 P1(*) 49 : REAL*8 CH(*) 50 : C 51 : ZI = (0.D0 , -1.D0) 52 : N2 = N*2 53 : C 54 : DO 1 I = 1 , N 55 : B(I ) = ZB(I) 56 : B(I+N) = ZB(I)*ZI 57 : U(I ) = ZU(I) 58 : U(I+N) = ZU(I)*ZI 59 : DO 2 J = 1 , N 60 : A(I ,J ) = ZA(I,J) 61 : A(I+N,J+N) = A(I,J) 62 : A(I+N,J ) = ZA(I,J)*ZI 63 : A(I ,J+N) = -A(I+N,J) 64 : 2 CONTINUE 65 : 1 CONTINUE 66 : C 67 : C 68 : C EPS = 1.D-8 69 : C 70 : C 71 : C ............................... R = B - A*U 72 : C 73 : DO 10 I = 1 , N2 74 : R(I) = B(I) 75 : DO 11 J = 1 , N2 76 : R(I) = R(I) - A(I,J)*U(J) 77 : 11 CONTINUE 78 : 10 CONTINUE 79 : C 80 : C ............................... P1 = P = R1 = R 81 : C 82 : DO 12 I = 1 , N2 83 : R1(I) = R(I) 84 : P (I) = R(I) 85 : P1(I) = R(I) 86 : 12 CONTINUE 87 : C 88 : NIT = 0 89 : 100 NIT = NIT + 1 90 : C ............................... CH = (A) * P 91 : C 92 : DO 20 I = 1 , N2 93 : CH (I) = 0.D0 94 : DO 21 J = 1 , N2 95 : CH (I) = CH(I) + A(I,J)*P(J) 96 : 21 CONTINUE 97 : 20 CONTINUE 98 : C 99 : C ............................... DEL1 = (R)T * R1 100 : C 101 : DEL1 = 0.D0 102 : DO 30 I = 1 , N2 103 : DEL1 = DEL1 + R(I)*R1(I) 104 : 30 CONTINUE 105 : C 106 : C ............................... ALFA = DEL1/(CH)T*P1 107 : C 108 : D = 0.D0 109 : DO 40 I = 1 , N2 110 : D = D + CH(I)*P1(I) 111 : 40 CONTINUE 112 : C 113 : ALFA = DEL1 / D 114 : C 115 : C ............................... ZU = ZU + ZLANDA*ZP 116 : C ............................... ZDELTA = ZDELTA + ZLANDA*ZH1 117 : C ............................... DELTA1 = (ZDELTA)T * ZDELTA 118 : C 119 : DEL2 = 0.D0 120 : DO 80 I = 1 , N2 121 : U(I) = U(I) + ALFA * P(I) 122 : R(I) = R(I) - ALFA * CH(I) 123 : 80 CONTINUE 124 : C 125 : DO 81 I = 1 , N2 126 : CH (I) = 0.D0 127 : DO 82 J = 1 , N2 128 : CH (I) = CH(I) + A(J,I)*P1(J) 129 : 82 CONTINUE 130 : C 131 : R1(I) = R1(I) - ALFA*CH(I) 132 : DEL2 = DEL2 + R(I)*R1(I) 133 : C 134 : 81 CONTINUE 135 : C 136 : BETA = DEL2 / DEL1 137 : C 138 : XNU = 0.D0 139 : XDU = 0.D0 140 : XNR = 0.D0 141 : XDR = 0.D0 142 : XNF = 0.D0 143 : XDF = 0.D0 144 : XNM = 0.D0 145 : XDM = 0.D0 146 : ND12 = N / 12 147 : DO 85 I = 1 , ND12 148 : II = (I-1)*12 149 : DO 86 J = 1 , 3 150 : XNU = XNU + ALFA*ALFA*(P(II+J )**2 + P(II+J +N)**2) 151 : XNR = XNR + ALFA*ALFA*(P(II+J+3)**2 + P(II+J+3+N)**2) 152 : XNF = XNF + ALFA*ALFA*(P(II+J+6)**2 + P(II+J+6+N)**2) 153 : XNM = XNM + ALFA*ALFA*(P(II+J+9)**2 + P(II+J+9+N)**2) 154 : C 155 : XDU = XDU + U(II+J )**2 + U(II+J +N)**2 156 : XDR = XDR + U(II+J+3)**2 + U(II+J+3+N)**2 157 : XDF = XDF + U(II+J+6)**2 + U(II+J+6+N)**2 158 : XDM = XDM + U(II+J+9)**2 + U(II+J+9+N)**2 159 : 86 CONTINUE 160 : 85 CONTINUE 161 : C 162 : C NN = (NNPOI-1)*12 + ICHAR 163 : C XNUM = P(NN)**2 +P(NN+N)**2 164 : C XDEN = U(NN)**2 +U(NN+N)**2 165 : C 166 : XST1 = SQRT(XNU/XDU) 167 : XST2 = SQRT(XNR/XDR) 168 : XST3 = SQRT(XNF/XDF) 169 : XST4 = SQRT(XNM/XDM) 170 : XSTOP = SQRT(XST1*XST1+XST2*XST2+XST3*XST3+XST4*XST4) 171 : C WRITE(6,*)'XST1 :',XST1 172 : C WRITE(6,*)'XST2 :',XST2 173 : C WRITE(6,*)'XST3 :',XST3 174 : C WRITE(6,*)'XST4 :',XST4 175 : C 176 : IF (XSTOP.LE.EPS) THEN 177 : DO 88 I = 1 , N 178 : ZU(I) = U(I) - ZI*U(I+N) 179 : 88 CONTINUE 180 : WRITE(6,*) 'VECTEUR SOLUTION :',NNPOI,ICHAR 181 : C DO 89 I = 1 , N 182 : C WRITE(6,*) I , ZU(I) 183 : C9 CONTINUE 184 : C 185 : RETURN 186 : END IF 187 : C 188 : C 189 : C ************************************************************ 190 : C 191 : DO 90 I = 1 , N2 192 : P (I) = BETA*P (I) + R (I) 193 : P1(I) = BETA*P1(I) + R1(I) 194 : 90 CONTINUE 195 : C 196 : GO TO 100 197 : C 198 : END 199 : 200 :
© Cast3M 2003 - All rights reserved.
Disclaimer