1 : C POLA2 SOURCE CHAT 05/01/13 02:17:10 5004 2 : SUBROUTINE POLA2(F,R,U,N) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : -INC CCOPTIO 6 : DIMENSION F(*),R(*),U(*) 7 : DIMENSION C(3,3),C2(3,3),UI(9) 8 : * 9 : JEBOUC=0 10 : IIMPI0=IIMPI 11 : 2020 JEBOUC=JEBOUC+1 12 : KERRE=0 13 : * 14 : * CALCUL DE C = U2 = FT*F 15 : * 16 : DO 1 I=1,N 17 : DO 1 J=1,N 18 : C(I,J)=0.D0 19 : DO 1 K=1,N 20 : KI= (K-1)*N+I 21 : KJ= (K-1)*N+J 22 : C(I,J)=C(I,J)+F(KI)*F(KJ) 23 : 1 CONTINUE 24 : * 25 : IF(IIMPI.EQ.199) THEN 26 : WRITE(6,77771) N 27 : 77771 FORMAT(2X,'POLA2 - N=',I3/) 28 : N2=N*N 29 : WRITE(6,77772) (F(I),I=1,N2) 30 : 77772 FORMAT(2X,'F '/(3(1X,1PE12.5))) 31 : ENDIF 32 : * 33 : IF(N.EQ.2) THEN 34 : * 35 : * CAS 2D 36 : * 37 : TRC=C(1,1)+C(2,2) 38 : RDTC=SQRT(C(1,1)*C(2,2)-C(1,2)*C(2,1)) 39 : DTU=RDTC 40 : TRU=SQRT(TRC+2.D0*RDTC) 41 : * 42 : IF(IIMPI.EQ.199) THEN 43 : WRITE(6,77773) TRU,DTU 44 : 77773 FORMAT(2X,'POLA2 - TRU= ',1PE12.5,2X,'DTU= ',1PE12.5/) 45 : ENDIF 46 : IF(TRU.EQ.0.D0.OR.DTU.EQ.0.D0) THEN 47 : WRITE(6,77883) TRU,DTU 48 : 77883 FORMAT(2X,'TRU=',1PE12.5,2X,'DTU=',1PE12.5/) 49 : KERRE=26 50 : GO TO 2021 51 : ENDIF 52 : * 53 : U(1)=(DTU+C(1,1))/TRU 54 : U(2)=C(1,2)/TRU 55 : U(3)=C(2,1)/TRU 56 : U(4)=(DTU+C(2,2))/TRU 57 : * 58 : UI(1)=(TRU-U(1))/DTU 59 : UI(2)=-U(2)/DTU 60 : UI(3)=-U(3)/DTU 61 : UI(4)=(TRU-U(4))/DTU 62 : * 63 : ELSE IF(N.EQ.3) THEN 64 : * 65 : * CAS 3D 66 : * 67 : * CAS FAUX 3D 68 : * 69 : IF(C(1,3).EQ.0.D0.AND.C(2,3).EQ.0.D0.AND. 70 : & C(3,1).EQ.0.D0.AND.C(3,2).EQ.0.D0) THEN 71 : 72 : TRC=C(1,1)+C(2,2) 73 : RDTC=SQRT(C(1,1)*C(2,2)-C(1,2)*C(2,1)) 74 : DTU=RDTC 75 : TRU=SQRT(TRC+2.D0*RDTC) 76 : * 77 : IF(TRU.EQ.0.D0.OR.DTU.EQ.0.D0.OR.F(9).EQ.0.D0) THEN 78 : KERRE=26 79 : GO TO 2021 80 : ENDIF 81 : * 82 : U(1)=(DTU+C(1,1))/TRU 83 : U(2)=C(1,2)/TRU 84 : U(3)=0.D0 85 : U(4)=C(2,1)/TRU 86 : U(5)=(DTU+C(2,2))/TRU 87 : U(6)=0.D0 88 : U(7)=0.D0 89 : U(8)=0.D0 90 : U(9)=F(9) 91 : * 92 : UI(1)=(TRU-U(1))/DTU 93 : UI(2)=-U(2)/DTU 94 : UI(3)=0.D0 95 : UI(4)=-U(4)/DTU 96 : UI(5)=(TRU-U(5))/DTU 97 : UI(6)=0.D0 98 : UI(7)=0.D0 99 : UI(8)=0.D0 100 : UI(9)=1.D0/F(9) 101 : * 102 : ELSE 103 : * 104 : * CAS VRAI 3D 105 : * 106 : DO 3 I=1,N 107 : DO 3 J=1,N 108 : C2(I,J)=0.D0 109 : DO 3 K=1,N 110 : C2(I,J)=C2(I,J)+C(I,K)*C(K,J) 111 : 3 CONTINUE 112 : TRC = C(1,1)+C(2,2)+C(3,3) 113 : TRC2 = C2(1,1)+C2(2,2)+C2(3,3) 114 : AUX=TRC**2 115 : P2C = (AUX - TRC2)/2.D0 116 : DTC = C(1,1)*(C(2,2)*C(3,3)-C(2,3)*C(3,2)) 117 : . -C(1,2)*(C(2,1)*C(3,3)-C(2,3)*C(3,1)) 118 : . +C(1,3)*(C(2,1)*C(3,2)-C(2,2)*C(3,1)) 119 : XK=AUX-3.D0*P2C 120 : TOL=1.D-6 121 : * 122 : IF(IIMPI.EQ.199) THEN 123 : WRITE(6,77764) TRC,TRC2,AUX,P2C 124 : 77764 FORMAT(2X,'POLA2 - TRC= ',1PE12.5,2X,'TRC2= ',1PE12.5/ 125 : . 2X,'AUX= ',1PE12.5,2X,'P2C=',1PE12.5/) 126 : WRITE(6,77774) XK,TOL 127 : 77774 FORMAT(2X,'POLA2 - XK= ',1PE12.5,2X,'TOL= ',1PE12.5/) 128 : ENDIF 129 : * 130 : IF(XK.LT.TOL) THEN 131 : XLAM=SQRT(TRC/3.D0) 131 : CALL ZERO(U,9,1) 131 : CALL ZERO(UI,9,1) 132 : U(1)=XLAM 133 : U(5)=XLAM 134 : U(9)=XLAM 135 : * 136 : UNXLAM=1.D0/XLAM 137 : UI(1)=UNXLAM 138 : UI(5)=UNXLAM 139 : UI(9)=UNXLAM 140 : * 141 : ELSE 142 : XL=TRC*(AUX-4.5D0*P2C)+ 13.5D0*DTC 143 : AUX2=XL/SQRT(XK**3) 144 : IF(IIMPI.EQ.199) THEN 145 : WRITE(6,77777) XL,TRC,AUX2,DTC 146 : 77777 FORMAT(2X,'POLA2 - XL =',1PE12.5,2X,'TRC=',1PE12.5, 147 : . 2X,'AUX2=',1PE12.5,2X,'DTC=',1PE12.5/) 148 : ENDIF 149 : IF(ABS(AUX2).GT.1.D0) THEN 150 : TOLEPS=1.D-10 151 : IF((ABS(AUX2)-1.D0).GE.TOLEPS) THEN 152 : ZZZ = ABS(AUX2)-1.D0 153 : WRITE(6,77884) ZZZ ,TOLEPS 154 : 77884 FORMAT(2X,'ZOB =',1PE12.5,2X,'TOLEPS=',1PE12.5/) 155 : KERRE=26 156 : GO TO 2021 157 : ELSE 158 : IF(AUX2.GT.0.D0) THEN 159 : AUX2 = AUX2 - TOLEPS 160 : ELSE 161 : AUX2 = AUX2 + TOLEPS 162 : ENDIF 163 : ENDIF 164 : ENDIF 165 : PHI=ACOS(AUX2) 166 : XLAM2=(TRC+2.D0*SQRT(XK)*COS(PHI/3.D0))/3.D0 167 : XLAM=SQRT(XLAM2) 168 : DTU=SQRT(DTC) 169 : TRU=XLAM+SQRT(TRC-XLAM2+(2.D0*DTU/XLAM)) 170 : P2U=(TRU**2-TRC)*0.5D0 171 : * 172 : FAC1=TRU*P2U-DTU 173 : FAC2=TRU*DTU 174 : FAC3=TRU**2-P2U 175 : * 176 : IF(IIMPI.EQ.199) THEN 177 : WRITE(6,77775) FAC1,DTU,TRU,P2U 178 : 77775 FORMAT(2X,'POLA2 - FAC1 =',1PE12.5,2X,'DTU=',1PE12.5, 179 : . 2X,'TRU=',1PE12.5,2X,'P2U=',1PE12.5/) 180 : ENDIF 181 : IF(FAC1.EQ.0.D0.OR.DTU.EQ.0.D0) THEN 182 : WRITE(6,77885) FAC1,DTU 183 : 77885 FORMAT(2X,'FAC1=',1PE12.5,2X,'DTU=',1PE12.5/) 184 : KERRE=26 185 : GO TO 2021 186 : ENDIF 187 : * 188 : DO 4 I=1,N 189 : IN=(I-1)*N 190 : DO 4 J=1,N 191 : IJ= IN+J 192 : IF(J.EQ.I) THEN 193 : U(IJ)=(FAC2+FAC3*C(I,I)-C2(I,I))/FAC1 194 : UI(IJ)=(P2U-TRU*U(IJ)+C(I,I))/DTU 195 : ELSE 196 : U(IJ)=(FAC3*C(I,J)-C2(I,J))/FAC1 197 : UI(IJ)=(-TRU*U(IJ)+C(I,J))/DTU 198 : ENDIF 199 : 4 CONTINUE 200 : * 201 : ENDIF 202 : ENDIF 203 : * 204 : ELSE 205 : KERRE=19 206 : GO TO 2021 207 : ENDIF 208 : * 209 : DO 2 I=1,N 210 : IN=(I-1)*N 211 : DO 2 J=1,N 212 : IJ= IN+J 213 : R(IJ)=0.D0 214 : DO 2 K=1,N 215 : IK= IN+K 216 : KJ= (K-1)*N+J 217 : R(IJ)=R(IJ)+F(IK)*UI(KJ) 218 : 2 CONTINUE 219 : * 220 : * PROOF 221 : * 222 : IF(IIMPI.EQ.199) THEN 223 : DO 7 I=1,N 224 : IN=(I-1)*N 225 : DO 7 J=1,N 226 : IJ= IN+J 227 : UI(IJ)=F(IJ) 228 : DO 7 K=1,N 229 : IK= IN+K 230 : KJ= (K-1)*N+J 231 : UI(IJ)=UI(IJ)-R(IK)*U(KJ) 232 : 7 CONTINUE 233 : * 234 : WRITE(6,77830) (UI(K),K=1,N2) 235 : 77830 FORMAT(2X,'POLA2- PROOF '/(6(1X,1PE12.5))) 236 : WRITE(6,77831) (R(K),K=1,N2) 237 : 77831 FORMAT(2X,'POLA2- R '/(3(1X,1PE12.5))) 238 : WRITE(6,77832) (U(K),K=1,N2) 239 : 77832 FORMAT(2X,'POLA2- U '/(3(1X,1PE12.5))) 240 : ENDIF 241 : * 242 : 2021 CONTINUE 243 : IF(KERRE.EQ.0) GO TO 9999 244 : * 245 : IF(JEBOUC.EQ.1.AND.IIMPI.EQ.1199) THEN 246 : IIMPI=199 247 : GO TO 2020 248 : ENDIF 249 : * 250 : 9999 CONTINUE 251 : IIMPI=IIMPI0 252 : * 252 : IF(KERRE.NE.0) CALL ERREUR(KERRE) 253 : * 254 : RETURN 255 : END 256 : 257 : 258 :
© Cast3M 2003 - All rights reserved.
Disclaimer