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

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