Télécharger pola2.eso

Retour à la liste

Numérotation des lignes :

  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)
  132. CALL ZERO(U,9,1)
  133. CALL ZERO(UI,9,1)
  134. U(1)=XLAM
  135. U(5)=XLAM
  136. U(9)=XLAM
  137. *
  138. UNXLAM=1.D0/XLAM
  139. UI(1)=UNXLAM
  140. UI(5)=UNXLAM
  141. UI(9)=UNXLAM
  142. *
  143. ELSE
  144. XL=TRC*(AUX-4.5D0*P2C)+ 13.5D0*DTC
  145. AUX2=XL/SQRT(XK**3)
  146. IF(IIMPI.EQ.199) THEN
  147. WRITE(6,77777) XL,TRC,AUX2,DTC
  148. 77777 FORMAT(2X,'POLA2 - XL =',1PE12.5,2X,'TRC=',1PE12.5,
  149. . 2X,'AUX2=',1PE12.5,2X,'DTC=',1PE12.5/)
  150. ENDIF
  151. IF(ABS(AUX2).GT.1.D0) THEN
  152. TOLEPS=1.D-10
  153. IF((ABS(AUX2)-1.D0).GE.TOLEPS) THEN
  154. ZZZ = ABS(AUX2)-1.D0
  155. WRITE(6,77884) ZZZ ,TOLEPS
  156. 77884 FORMAT(2X,'ZOB =',1PE12.5,2X,'TOLEPS=',1PE12.5/)
  157. KERRE=26
  158. GO TO 2021
  159. ELSE
  160. IF(AUX2.GT.0.D0) THEN
  161. AUX2 = AUX2 - TOLEPS
  162. ELSE
  163. AUX2 = AUX2 + TOLEPS
  164. ENDIF
  165. ENDIF
  166. ENDIF
  167. PHI=ACOS(AUX2)
  168. XLAM2=(TRC+2.D0*SQRT(XK)*COS(PHI/3.D0))/3.D0
  169. XLAM=SQRT(XLAM2)
  170. DTU=SQRT(DTC)
  171. TRU=XLAM+SQRT(TRC-XLAM2+(2.D0*DTU/XLAM))
  172. P2U=(TRU**2-TRC)*0.5D0
  173. *
  174. FAC1=TRU*P2U-DTU
  175. FAC2=TRU*DTU
  176. FAC3=TRU**2-P2U
  177. *
  178. IF(IIMPI.EQ.199) THEN
  179. WRITE(6,77775) FAC1,DTU,TRU,P2U
  180. 77775 FORMAT(2X,'POLA2 - FAC1 =',1PE12.5,2X,'DTU=',1PE12.5,
  181. . 2X,'TRU=',1PE12.5,2X,'P2U=',1PE12.5/)
  182. ENDIF
  183. IF(FAC1.EQ.0.D0.OR.DTU.EQ.0.D0) THEN
  184. WRITE(6,77885) FAC1,DTU
  185. 77885 FORMAT(2X,'FAC1=',1PE12.5,2X,'DTU=',1PE12.5/)
  186. KERRE=26
  187. GO TO 2021
  188. ENDIF
  189. *
  190. DO 4 I=1,N
  191. IN=(I-1)*N
  192. DO 4 J=1,N
  193. IJ= IN+J
  194. IF(J.EQ.I) THEN
  195. U(IJ)=(FAC2+FAC3*C(I,I)-C2(I,I))/FAC1
  196. UI(IJ)=(P2U-TRU*U(IJ)+C(I,I))/DTU
  197. ELSE
  198. U(IJ)=(FAC3*C(I,J)-C2(I,J))/FAC1
  199. UI(IJ)=(-TRU*U(IJ)+C(I,J))/DTU
  200. ENDIF
  201. 4 CONTINUE
  202. *
  203. ENDIF
  204. ENDIF
  205. *
  206. ELSE
  207. KERRE=19
  208. GO TO 2021
  209. ENDIF
  210. *
  211. DO 2 I=1,N
  212. IN=(I-1)*N
  213. DO 2 J=1,N
  214. IJ= IN+J
  215. R(IJ)=0.D0
  216. DO 2 K=1,N
  217. IK= IN+K
  218. KJ= (K-1)*N+J
  219. R(IJ)=R(IJ)+F(IK)*UI(KJ)
  220. 2 CONTINUE
  221. *
  222. * PROOF
  223. *
  224. IF(IIMPI.EQ.199) THEN
  225. DO 7 I=1,N
  226. IN=(I-1)*N
  227. DO 7 J=1,N
  228. IJ= IN+J
  229. UI(IJ)=F(IJ)
  230. DO 7 K=1,N
  231. IK= IN+K
  232. KJ= (K-1)*N+J
  233. UI(IJ)=UI(IJ)-R(IK)*U(KJ)
  234. 7 CONTINUE
  235. *
  236. WRITE(6,77830) (UI(K),K=1,N2)
  237. 77830 FORMAT(2X,'POLA2- PROOF '/(6(1X,1PE12.5)))
  238. WRITE(6,77831) (R(K),K=1,N2)
  239. 77831 FORMAT(2X,'POLA2- R '/(3(1X,1PE12.5)))
  240. WRITE(6,77832) (U(K),K=1,N2)
  241. 77832 FORMAT(2X,'POLA2- U '/(3(1X,1PE12.5)))
  242. ENDIF
  243. *
  244. 2021 CONTINUE
  245. IF(KERRE.EQ.0) GO TO 9999
  246. *
  247. IF(JEBOUC.EQ.1.AND.IIMPI.EQ.1199) THEN
  248. IIMPI=199
  249. GO TO 2020
  250. ENDIF
  251. *
  252. 9999 CONTINUE
  253. IIMPI=IIMPI0
  254. *
  255. IF(KERRE.NE.0) CALL ERREUR(KERRE)
  256. *
  257. RETURN
  258. END
  259.  
  260.  
  261.  
  262.  

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