1 : C ELLP54 SOURCE CHAT 05/01/12 23:35:51 5004 2 : SUBROUTINE ELLP54(ZA,ZB,ZU,N,N2,NNPOI,ICHAR,XUR,XFM,NIT, 3 : * A,B,U,R,DELTA,P,CH,CH1,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 DELTA(*) 47 : REAL*8 P(*) 48 : REAL*8 CH(*) 49 : REAL*8 CH1(*) 50 : C 51 : ZI = (0.D0 , -1.D0) 52 : N2 = N*2 53 : NP4= N/6 54 : C 55 : DO 1 I = 1 , N 56 : B(I ) = ZB(I) 57 : B(I+N) = ZB(I)*ZI 58 : U(I ) = ZU(I) 59 : U(I+N) = ZU(I)*ZI 60 : DO 2 J = 1 , N 61 : A(I ,J ) = ZA(I,J) 62 : A(I+N,J+N) = A(I,J) 63 : A(I+N,J ) = ZA(I,J)*ZI 64 : A(I ,J+N) = -A(I+N,J) 65 : 2 CONTINUE 66 : 1 CONTINUE 67 : 68 : DO 6 I = 1 , NP4 69 : DO 7 J = 1 , 6 70 : U((I-1)*12+J ) = U((I-1)*12+J ) / XUR 71 : U((I-1)*12+J+6) = U((I-1)*12+J+6) / XFM 72 : 7 CONTINUE 73 : 6 CONTINUE 74 : C 75 : DO 3 J = 1 , NP4 76 : DO 5 I = 1 , N 77 : DO 4 JJ = 1 , 6 78 : A(I ,(J-1)*12+JJ ) = A(I,(J-1)*12+JJ )*XUR 79 : A(I ,(J-1)*12+JJ+6 ) = A(I,(J-1)*12+JJ+6)*XFM 80 : 4 CONTINUE 81 : 5 CONTINUE 82 : 3 CONTINUE 83 : C 84 : C 85 : C EPS = 1.D-8 86 : C 87 : C 88 : C ............................... R = B - A*U 89 : C 90 : DO 10 I = 1 , N2 91 : R(I) = B(I) 92 : DO 11 J = 1 , N2 93 : R(I) = R(I) - A(I,J)*U(J) 94 : 11 CONTINUE 95 : 10 CONTINUE 96 : C 97 : C 98 : C ............................... ZDELTA = (ZA)T * ZR 99 : C 100 : DO 20 I = 1 , N2 101 : DELTA (I) = 0.D0 102 : DO 21 J = 1 , N2 103 : DELTA (I) = DELTA(I) + A(J,I)*R(J) 104 : 21 CONTINUE 105 : 20 CONTINUE 106 : C 107 : C ............................... DELTA = (ZDELTA)T * ZDELTA 108 : C 109 : DEL = 0.D0 110 : NIT = 0 111 : C 112 : DO 30 I = 1 , N2 113 : DEL = DEL + DELTA(I)**2 114 : 30 CONTINUE 115 : C ............................... ZP = ZDELTA 116 : C 117 : DO 40 I = 1 , N2 118 : P(I) = DELTA(I) 119 : 40 CONTINUE 120 : C 121 : 100 NIT = NIT + 1 122 : WRITE(6,*)'********************* ',NIT 123 : DO 41 I = 1 , N 124 : WRITE(6,*)I,' ',U(I),' ',U(I+N) 125 : 41 CONTINUE 126 : C 127 : C 128 : C ............................... ZH = ZA * ZP 129 : C 130 : DO 50 I = 1 , N2 131 : CH(I) = 0.D0 132 : DO 51 J = 1 , N2 133 : CH(I) = CH(I) + A(I,J)*P(J) 134 : 51 CONTINUE 135 : 50 CONTINUE 136 : C 137 : C ............................... ZH1 = (ZA)T * ZH 138 : C 139 : DO 60 I = 1 , N2 140 : CH1(I) = 0.D0 141 : DO 61 J = 1 , N2 142 : CH1(I) = CH1(I) + A(J,I)*CH(J) 143 : 61 CONTINUE 144 : 60 CONTINUE 145 : C 146 : C ............................... ZLANDA = DELTA / ((ZP)T * ZH1) 147 : C 148 : D = 0.D0 149 : DO 70 I = 1 , N2 150 : D = D + CH1(I)*P(I) 151 : 70 CONTINUE 152 : C 153 : CLANDA = DEL / D 154 : C 155 : C ............................... ZU = ZU + ZLANDA*ZP 156 : C ............................... ZDELTA = ZDELTA + ZLANDA*ZH1 157 : C ............................... DELTA1 = (ZDELTA)T * ZDELTA 158 : C 159 : DEL1 = 0.D0 160 : DO 80 I = 1 , N2 161 : U(I) = U(I) + CLANDA * P(I) 162 : DELTA(I) = DELTA(I) - CLANDA * CH1(I) 163 : DEL1 = DEL1 + DELTA(I)**2 164 : 80 CONTINUE 165 : C 166 : C ************** TEST DE CONVERGENCE ***************** 167 : C 168 : XSTOP = 0.D0 169 : XNU = 0.D0 170 : XDU = 0.D0 171 : XNR = 0.D0 172 : XDR = 0.D0 173 : XNF = 0.D0 174 : XDF = 0.D0 175 : XNM = 0.D0 176 : XDM = 0.D0 177 : C 178 : ND12 = N / 12 179 : DO 85 I = 1 , ND12 180 : II = (I-1) * 12 181 : DO 86 J = 1 , 3 182 : XNU = XNU + CLANDA*CLANDA*(P(II+J )**2+P(II+J +N)**2) 183 : XNR = XNR + CLANDA*CLANDA*(P(II+J+3)**2+P(II+J+3+N)**2) 184 : XNF = XNF + CLANDA*CLANDA*(P(II+J+6)**2+P(II+J+6+N)**2) 185 : XNM = XNM + CLANDA*CLANDA*(P(II+J+9)**2+P(II+J+9+N)**2) 186 : XDU = XDU + U(II+J )**2 + U(II+J +N)**2 187 : XDR = XDR + U(II+J+3)**2 + U(II+J+3+N)**2 188 : XDF = XDF + U(II+J+6)**2 + U(II+J+6+N)**2 189 : XDM = XDM + U(II+J+9)**2 + U(II+J+9+N)**2 190 : 86 CONTINUE 191 : 85 CONTINUE 192 : C 193 : C NN = (NNPOI-1)*12+ICHAR 194 : C XNUM = (CLANDA*P(NN))**2+(CLANDA*P(NN+N))**2 195 : C XDEN = U(NN)**2 + U(NN+N)**2 196 : XST1 = SQRT(XNU/XDU) 197 : XST2 = SQRT(XNR/XDR) 198 : XST3 = SQRT(XNF/XDF) 199 : XST4 = SQRT(XNM/XDM) 200 : XSTOP = SQRT(XST1*XST1+XST2*XST2+XST3*XST3+XST4*XST4) 201 : C WRITE(6,*)'XST1 :',XST1 202 : C WRITE(6,*)'XST2 :',XST2 203 : C WRITE(6,*)'XST3 :',XST3 204 : C WRITE(6,*)'XST4 :',XST4 205 : C 206 : IF (XSTOP.LE.EPS) THEN 207 : WRITE(6 ,*)'SOLUTION DU SYSTEME :',NNPOI,ICHAR 208 : DO 8 I = 1 , NP4 209 : DO 9 J = 1 , 6 210 : U((I-1)*12+J ) = U((I-1)*12+J ) * XUR 211 : U((I-1)*12+J+6) = U((I-1)*12+J+6) * XFM 212 : 9 CONTINUE 213 : 8 CONTINUE 214 : C 215 : DO 84 I = 1 , N 216 : ZU(I) = U(I) - ZI*U(I+N) 217 : 84 CONTINUE 218 : C DO 89 I = 1 , N 219 : C WRITE(21 ,*) I,ZU(I) 220 : C9 CONTINUE 221 : C 222 : C WRITE(6,*)'--------------------------' 223 : RETURN 224 : END IF 225 : C 226 : C ************************************************************ 227 : ALFA = DEL1 / DEL 228 : DEL = DEL1 229 : C 230 : DO 90 I = 1 , N2 231 : P(I) = ALFA*P(I) + DELTA(I) 232 : 90 CONTINUE 233 : C 234 : GO TO 100 235 : C 236 : END 237 : 238 :
© Cast3M 2003 - All rights reserved.
Disclaimer