Télécharger la plaquette Cast3M


fr en

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