Télécharger la plaquette Cast3M


fr en

Télécharger tens2.eso

Retour à la liste

Numérotation des lignes : 

tens2

  1. C TENS2     SOURCE    GOUNAND   25/10/23    21:15:09     12386          
  2.       SUBROUTINE TENS2(IOTENS,IKAS,A,D,R)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : TENS2
  7. C DESCRIPTION : Opérations sur des tenseurs (unaires pour l'instant)
  8. C               Ici on gère les opérations à l'échelle du noeud
  9. C
  10. C
  11. C LANGAGE     : ESOPE
  12. C AUTEUR      : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA)
  13. C               mel : gounand@semt2.smts.cea.fr
  14. C***********************************************************************
  15. C APPELES          :
  16. C APPELES (E/S)    :
  17. C APPELES (BLAS)   :
  18. C APPELES (CALCUL) :
  19. C APPELE PAR       :
  20. C***********************************************************************
  21. C SYNTAXE GIBIANE    :
  22. C ENTREES            :  A si IOTENS.NE.12
  23. C                       D et R si IOTENS.EQ.12
  24. C ENTREES/SORTIES    :
  25. C SORTIES            :  R si IOTENS.NE.12
  26. C                       D si IOTENS.EQ.11
  27. C                       A si IOTENS.EQ.12
  28. C***********************************************************************
  29. C VERSION    : v1, 28/08/2024, version initiale
  30. C HISTORIQUE : v1, 28/08/2024, creation
  31. C HISTORIQUE :
  32. C HISTORIQUE :
  33. C***********************************************************************
  34. -INC PPARAM
  35. -INC CCOPTIO
  36. -INC CCREEL
  37. *
  38. C   ENTREEES
  39. C        A(3,3) = MATRICE A DIAGONALISER
  40. C   SORTIES
  41. C        D(3)   =  VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3)
  42. C        S(3,3) =  VECTEURS PROPRES
  43. C        R(3,3) =  RESULTAT DES OPERATIONS TENSORIELLES
  44.       DIMENSION A(3,3),D(3),S(3,3),R(3,3)
  45.       LOGICAL LVERIF,LDBG
  46.       DATA ((S(I,J),I=1,3),J=1,3) /9*0.D0/
  47. *                     1       2        3     4        5        6      7
  48. *     DATA MOTENS/'NORM2','NORMINF','DET','TRACE','INVERSE','IDEN' ,'LOG'
  49. *              8     9      10        11         12         13
  50. *     $     ,'EXP','INVS','ABSOLU','PRINCIPA','RECOMPOS','TRANSPOS'/
  51. *
  52. *
  53. * Executable statements
  54. *
  55.       lverif=.false.
  56.       ldbg=.false.
  57.       if (ldbg) then
  58.          WRITE(IOIMP,*) 'TENS2 : IOTENS,IKAS=',IOTENS,IKAS
  59.          IF (IOTENS.NE.12) THEN
  60.             write(ioimp,*) 'Input  A,I,J=',IDIM,IDIM
  61.             write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
  62.          ELSE
  63.             write(ioimp,*) 'Input  D,J=',IDIM
  64.             write(ioimp,*) (D(J),J=1,IDIM)
  65.             write(ioimp,*) 'Input  R,I,J=',IDIM,IDIM
  66.             write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  67.          ENDIF
  68.       endif
  69. *
  70.       IF (IOTENS.GT.6.AND.IOTENS.LT.12) THEN
  71. *
  72. *     Décomposition spectrale (positive pour les cas scalaires et
  73. *     vecteurs)
  74. *
  75.          IF (IKAS.EQ.1.OR.(IKAS.EQ.3.AND.IDIM.EQ.1)) THEN
  76. * Difficile de donner un sens...
  77. *            D(1)=ABS(A(1,1))
  78. *            S(1,1)=SIGN(1.D0,A(1,1))
  79.             D(1)=A(1,1)
  80.             S(1,1)=1.D0
  81.             JD=1
  82.          ELSEIF (IKAS.EQ.2) THEN
  83. * Pas de sens pour les vecteurs
  84. c$$$            XNO=0.D0
  85. c$$$            DO I=1,IDIM
  86. c$$$               XNO=XNO+A(I,1)**2
  87. c$$$            ENDDO
  88. c$$$            XNA=SQRT(XNO)
  89. c$$$            D(1)=XNA
  90. c$$$            DO I=1,IDIM
  91. c$$$               S(I,1)=A(I,1)/XNA
  92. c$$$            ENDDO
  93.             CALL ERREUR(980)
  94.             RETURN
  95. c$$$         write(ioimp,*) 'A,I,J=',IDIM,1
  96. c$$$         write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,1)
  97. c$$$         write(ioimp,*) 'D,J=',JDIM
  98. c$$$         write(ioimp,*) (D(J),J=1,JDIM)
  99. c$$$         write(ioimp,*) 'S,I,J=',IDIM,1
  100. c$$$         write(ioimp,*) ((S(I,J),I=1,IDIM),J=1,1)
  101.  
  102.          ELSEIF (IKAS.EQ.3) THEN
  103.             CALL JACOB3(A,IDIM,D,S)
  104.             IF (IERR.NE.0) RETURN
  105.             JD=IDIM
  106.             if (lverif) then
  107.                write(ioimp,*) 'A,I,J=',IDIM,IDIM
  108.                write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
  109.                write(ioimp,*) 'D,J=',IDIM
  110.                write(ioimp,*) (D(J),J=1,IDIM)
  111.                write(ioimp,*) 'S,I,J=',IDIM,IDIM
  112.                write(ioimp,*) ((S(I,J),I=1,IDIM),J=1,IDIM)
  113. *     Verif
  114.                DO K=1,IDIM
  115. *                  DO I=K,IDIM
  116.                   DO I=1,IDIM
  117.                      R(I,K)=0.D0
  118.                      DO J=1,IDIM
  119. *                        R(I,K)=R(I,K)+D(J)*S(J,I)*S(J,K)
  120.                         R(I,K)=R(I,K)+D(J)*S(I,J)*S(K,J)
  121.                      ENDDO
  122.                   ENDDO
  123.                ENDDO
  124.                write(ioimp,*) 'A2,I,J=',IDIM,IDIM
  125.                write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  126. c$$$               DO J=1,IDIM
  127. c$$$                  DO I=1,IDIM
  128. c$$$                     V(I)=0.D0
  129. c$$$                     DO K=1,IDIM
  130. c$$$                        V(I)=V(I)+A(I,K)*S(K,J)
  131. c$$$                     ENDDO
  132. c$$$                  ENDDO
  133. c$$$                  write(ioimp,*) 'VJ(I)=AIK.SKJ,J',J
  134. c$$$                  write(ioimp,*) (V(I),I=1,IDIM)
  135. c$$$                  write(ioimp,*) 'VJ(I)=D(J).SIJ,J',J
  136. c$$$                  write(ioimp,*) (S(I,J)*D(J),I=1,IDIM)
  137. c$$$               ENDDO
  138.             endif
  139.          ELSE
  140. *     D'abord faire la décomposition polaire avant la décomposition
  141. *     spectrale dans le cas non-symetrique (mais inversible).
  142. *     Pas encore fait
  143.             CALL ERREUR(969)
  144.             RETURN
  145.          ENDIF
  146.       ENDIF
  147. *
  148. * Les opérations tensorielles proprement dites
  149. *
  150.       IF (IOTENS.EQ.1) THEN
  151.          IF (IKAS.EQ.1) THEN
  152.             XNORM=ABS(A(1,1))
  153.          ELSE
  154.             IF (IKAS.EQ.2) THEN
  155.                ID=IDIM
  156.                JD=1
  157.             ELSE
  158.                ID=IDIM
  159.                JD=IDIM
  160.             ENDIF
  161.             XNORM=0.D0
  162.             DO J=1,JD
  163.                DO I=1,ID
  164.                   XNORM=XNORM+A(I,J)**2
  165.                ENDDO
  166.             ENDDO
  167.             XNORM=SQRT(XNORM)
  168.          ENDIF
  169.          R(1,1)=XNORM
  170.       ELSEIF (IOTENS.EQ.2) THEN
  171.          IF (IKAS.EQ.1) THEN
  172.             XNORM=ABS(A(1,1))
  173.          ELSE
  174.             IF (IKAS.EQ.2) THEN
  175.                ID=IDIM
  176.                JD=1
  177.             ELSE
  178.                ID=IDIM
  179.                JD=IDIM
  180.             ENDIF
  181.             XNORM=0.D0
  182.             DO J=1,JD
  183.                DO I=1,ID
  184.                   XNORM=max(xnorm,abs(A(I,J)))
  185.                ENDDO
  186.             ENDDO
  187.          ENDIF
  188.          R(1,1)=XNORM
  189.       ELSEIF (IOTENS.EQ.3.OR.IOTENS.EQ.5) THEN
  190. *
  191.          IF (IKAS.EQ.1) THEN
  192.             JD=1
  193.          ELSE
  194.             JD=IDIM
  195.          ENDIF
  196.          DO J=1,JD
  197.             DO I=1,JD
  198.                R(I,J)=A(I,J)
  199.             ENDDO
  200.          ENDDO
  201. * D sert de vecteur de travail ici
  202.          CALL INVERE(R,JD,3,D,DET)
  203.          IF(IOTENS.EQ.3) THEN
  204.             R(1,1)=DET
  205.          ELSE
  206.             IF (DET.EQ.XZERO) THEN
  207. * Matrice singuliere. Numero de ligne =%i1
  208.                INTERR(1)=0
  209.                CALL ERREUR(49)
  210. *            write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  211.                RETURN
  212.             ENDIF
  213.          ENDIF
  214.       ELSEIF (IOTENS.EQ.4) THEN
  215.          IF (IKAS.EQ.1) THEN
  216.             XTRACE=A(1,1)
  217.          ELSE
  218.             XTRACE=0.D0
  219.             DO I=1,IDIM
  220.                XTRACE=XTRACE+A(I,I)
  221.             ENDDO
  222.          ENDIF
  223.          R(1,1)=XTRACE
  224.       ELSEIF (IOTENS.EQ.6) THEN
  225.          DO i=1,idim
  226.             DO j=1,idim
  227.                R(J,I)=A(J,I)
  228.            ENDDO
  229.          ENDDO
  230.       ELSEIF (IOTENS.EQ.7) THEN
  231.          DO J=1,JD
  232.             IF (D(J).GE.XPETIT) THEN
  233.                D(J)=LOG(D(J))
  234.             ELSE
  235.                REAERR(1)=D(J)
  236.                REAERR(2)=XPETIT
  237.                CALL ERREUR(1063)
  238.                RETURN
  239.             ENDIF
  240.          ENDDO
  241.       ELSEIF (IOTENS.EQ.8) THEN
  242. *         write(ioimp,*) 'TENS2 D(J) avant:',(D(idx),idx=1,jd)
  243.          DO J=1,JD
  244.             D(J)=EXP(D(J))
  245.          ENDDO
  246. *         write(ioimp,*) 'TENS2 D(J) apres:',(D(idx),idx=1,jd)
  247.       ELSEIF (IOTENS.EQ.9) THEN
  248.          DO J=1,JD
  249.             IF (ABS(D(J)).GE.XPETIT) THEN
  250.                D(J)=1.D0/D(J)
  251.             ELSE
  252.                CALL ERREUR(835)
  253.                RETURN
  254.             ENDIF
  255.          ENDDO
  256.       ELSEIF (IOTENS.EQ.10) THEN
  257.          DO J=1,JD
  258.             D(J)=ABS(D(J))
  259.          ENDDO
  260.       ELSEIF (IOTENS.EQ.11) THEN
  261.             DO J=1,JD
  262.                DO I=1,JD
  263.                R(I,J)=S(I,J)
  264.             ENDDO
  265.          ENDDO
  266.       ELSEIF (IOTENS.EQ.12) THEN
  267.          DO i=1,idim
  268.             DO j=1,idim
  269.                S(J,I)=R(J,I)
  270.            ENDDO
  271.          ENDDO
  272.       ELSEIF (IOTENS.EQ.13) THEN
  273.          DO i=1,idim
  274.             DO j=1,idim
  275.                R(I,J)=A(J,I)
  276.            ENDDO
  277.          ENDDO
  278.       ELSE
  279.          write(ioimp,*) 'iotens=',iotens,' not implemented'
  280.          MOTERR(1:8)='TENS2'
  281.          call erreur(1039)
  282.          return
  283.       ENDIF
  284. * Recomposition spectrale
  285.       IF (IOTENS.GT.6.AND.IOTENS.NE.11.AND.IOTENS.LT.13) THEN
  286. c$$$            write(ioimp,*) 'D,J=',JDIM
  287. c$$$            write(ioimp,*) (D(J),J=1,JDIM)
  288. c$$$            write(ioimp,*) 'IDIM,IKAS=',IDIM,IKAS
  289.          IF (IKAS.EQ.1.OR.((IKAS.EQ.3.OR.IKAS.EQ.5).AND.IDIM.EQ.1)) THEN
  290.             R(1,1)=D(1)*S(1,1)
  291.          ELSEIF (IKAS.EQ.2) THEN
  292. c$$$            DO I=1,IDIM
  293. c$$$               R(I,1)=D(1)*S(I,1)
  294. c$$$            ENDDO
  295.             CALL ERREUR(980)
  296.             RETURN
  297.          ELSEIF (IKAS.EQ.3) THEN
  298.             DO K=1,IDIM
  299. * Tenir compte de la symétrie               DO I=1,IDIM
  300.                DO I=K,IDIM
  301.                   R(I,K)=0.D0
  302.                   DO J=1,IDIM
  303. *     faux!                        R(I,K)=R(I,K)+D(J)*S(J,I)*S(J,K)
  304.                      R(I,K)=R(I,K)+D(J)*S(I,J)*S(K,J)
  305.                   ENDDO
  306.                   IF (I.NE.K) R(K,I)=R(I,K)
  307.                ENDDO
  308.             ENDDO
  309.          ELSEIF (IKAS.EQ.5) THEN
  310.             DO K=1,IDIM
  311. * Tenir compte de la symétrie               DO I=1,IDIM
  312.                DO I=K,IDIM
  313.                   A(I,K)=0.D0
  314.                   DO J=1,IDIM
  315.                      A(I,K)=A(I,K)+D(J)*S(I,J)*S(K,J)
  316.                   ENDDO
  317.                   IF (I.NE.K) A(K,I)=A(I,K)
  318.                ENDDO
  319.             ENDDO
  320.          ELSE
  321. *     Faire la recomposition polaire après la recomposition
  322. *     spectrale dans le cas non-symetrique (mais inversible).
  323. *     Pas encore fait
  324.             CALL ERREUR(969)
  325.             RETURN
  326.          ENDIF
  327.       ENDIF
  328.       if (ldbg) then
  329.          IF (IOTENS.NE.12) THEN
  330.             write(ioimp,*) 'R,I,J=',IDIM,IDIM
  331.             write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  332.          ELSE
  333.             write(ioimp,*) 'A,I,J=',IDIM,IDIM
  334.             write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
  335.          ENDIF
  336.          IF (IOTENS.EQ.11) THEN
  337.             write(ioimp,*) 'D,J=',IDIM
  338.             write(ioimp,*) (D(J),J=1,IDIM)
  339.          ENDIF
  340.       endif
  341. *
  342. * Normal termination
  343. *
  344.       RETURN
  345. *
  346. * Format handling
  347. *
  348. *
  349. * Error handling
  350. *
  351. *
  352. * End of subroutine TENS2
  353. *
  354.       END
  355.  
  356.  

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