Télécharger pola2.eso

Retour à la liste

Numérotation des lignes :

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

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