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   26/01/16    21:15:18     12450          
  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. *              14
  52. *     $     ,'VALP'/
  53. *
  54. *
  55. * Executable statements
  56. *
  57.       lverif=.false.
  58.       ldbg=.false.
  59.       if (ldbg) then
  60.          WRITE(IOIMP,*) 'TENS2 : IOTENS,IKAS=',IOTENS,IKAS
  61.          IF (IOTENS.NE.12) THEN
  62.             write(ioimp,*) 'Input  A,I,J=',IDIM,IDIM
  63.             write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
  64.          ELSE
  65.             write(ioimp,*) 'Input  D,J=',IDIM
  66.             write(ioimp,*) (D(J),J=1,IDIM)
  67.             write(ioimp,*) 'Input  R,I,J=',IDIM,IDIM
  68.             write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  69.          ENDIF
  70.       endif
  71. *
  72.       IF (IOTENS.GT.6.AND.IOTENS.LT.12) THEN
  73. *
  74. *     Décomposition spectrale (positive pour les cas scalaires et
  75. *     vecteurs)
  76. *
  77.          IF (IKAS.EQ.1.OR.(IKAS.EQ.3.AND.IDIM.EQ.1)) THEN
  78. * Difficile de donner un sens...
  79. *            D(1)=ABS(A(1,1))
  80. *            S(1,1)=SIGN(1.D0,A(1,1))
  81.             D(1)=A(1,1)
  82.             S(1,1)=1.D0
  83.             JD=1
  84.          ELSEIF (IKAS.EQ.2) THEN
  85. * Pas de sens pour les vecteurs
  86. c$$$            XNO=0.D0
  87. c$$$            DO I=1,IDIM
  88. c$$$               XNO=XNO+A(I,1)**2
  89. c$$$            ENDDO
  90. c$$$            XNA=SQRT(XNO)
  91. c$$$            D(1)=XNA
  92. c$$$            DO I=1,IDIM
  93. c$$$               S(I,1)=A(I,1)/XNA
  94. c$$$            ENDDO
  95.             CALL ERREUR(980)
  96.             RETURN
  97. c$$$         write(ioimp,*) 'A,I,J=',IDIM,1
  98. c$$$         write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,1)
  99. c$$$         write(ioimp,*) 'D,J=',JDIM
  100. c$$$         write(ioimp,*) (D(J),J=1,JDIM)
  101. c$$$         write(ioimp,*) 'S,I,J=',IDIM,1
  102. c$$$         write(ioimp,*) ((S(I,J),I=1,IDIM),J=1,1)
  103.  
  104.          ELSEIF (IKAS.EQ.3) THEN
  105.             CALL JACOB3(A,IDIM,D,S)
  106.             IF (IERR.NE.0) RETURN
  107.             JD=IDIM
  108.             if (lverif) then
  109.                write(ioimp,*) 'A,I,J=',IDIM,IDIM
  110.                write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
  111.                write(ioimp,*) 'D,J=',IDIM
  112.                write(ioimp,*) (D(J),J=1,IDIM)
  113.                write(ioimp,*) 'S,I,J=',IDIM,IDIM
  114.                write(ioimp,*) ((S(I,J),I=1,IDIM),J=1,IDIM)
  115. *     Verif
  116.                DO K=1,IDIM
  117. *                  DO I=K,IDIM
  118.                   DO I=1,IDIM
  119.                      R(I,K)=0.D0
  120.                      DO J=1,IDIM
  121. *                        R(I,K)=R(I,K)+D(J)*S(J,I)*S(J,K)
  122.                         R(I,K)=R(I,K)+D(J)*S(I,J)*S(K,J)
  123.                      ENDDO
  124.                   ENDDO
  125.                ENDDO
  126.                write(ioimp,*) 'A2,I,J=',IDIM,IDIM
  127.                write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  128. c$$$               DO J=1,IDIM
  129. c$$$                  DO I=1,IDIM
  130. c$$$                     V(I)=0.D0
  131. c$$$                     DO K=1,IDIM
  132. c$$$                        V(I)=V(I)+A(I,K)*S(K,J)
  133. c$$$                     ENDDO
  134. c$$$                  ENDDO
  135. c$$$                  write(ioimp,*) 'VJ(I)=AIK.SKJ,J',J
  136. c$$$                  write(ioimp,*) (V(I),I=1,IDIM)
  137. c$$$                  write(ioimp,*) 'VJ(I)=D(J).SIJ,J',J
  138. c$$$                  write(ioimp,*) (S(I,J)*D(J),I=1,IDIM)
  139. c$$$               ENDDO
  140.             endif
  141.          ELSE
  142. *     D'abord faire la décomposition polaire avant la décomposition
  143. *     spectrale dans le cas non-symetrique (mais inversible).
  144. *     Pas encore fait
  145.             CALL ERREUR(969)
  146.             RETURN
  147.          ENDIF
  148.       ENDIF
  149. *
  150. * Les opérations tensorielles proprement dites
  151. *
  152.       IF (IOTENS.EQ.1) THEN
  153.          IF (IKAS.EQ.1) THEN
  154.             XNORM=ABS(A(1,1))
  155.          ELSE
  156.             IF (IKAS.EQ.2) THEN
  157.                ID=IDIM
  158.                JD=1
  159.             ELSE
  160.                ID=IDIM
  161.                JD=IDIM
  162.             ENDIF
  163.             XNORM=0.D0
  164.             DO J=1,JD
  165.                DO I=1,ID
  166.                   XNORM=XNORM+A(I,J)**2
  167.                ENDDO
  168.             ENDDO
  169.             XNORM=SQRT(XNORM)
  170.          ENDIF
  171.          R(1,1)=XNORM
  172.       ELSEIF (IOTENS.EQ.2) THEN
  173.          IF (IKAS.EQ.1) THEN
  174.             XNORM=ABS(A(1,1))
  175.          ELSE
  176.             IF (IKAS.EQ.2) THEN
  177.                ID=IDIM
  178.                JD=1
  179.             ELSE
  180.                ID=IDIM
  181.                JD=IDIM
  182.             ENDIF
  183.             XNORM=0.D0
  184.             DO J=1,JD
  185.                DO I=1,ID
  186.                   XNORM=max(xnorm,abs(A(I,J)))
  187.                ENDDO
  188.             ENDDO
  189.          ENDIF
  190.          R(1,1)=XNORM
  191.       ELSEIF (IOTENS.EQ.3.OR.IOTENS.EQ.5) THEN
  192. *
  193.          IF (IKAS.EQ.1) THEN
  194.             JD=1
  195.          ELSE
  196.             JD=IDIM
  197.          ENDIF
  198.          DO J=1,JD
  199.             DO I=1,JD
  200.                R(I,J)=A(I,J)
  201.             ENDDO
  202.          ENDDO
  203. * D sert de vecteur de travail ici
  204.          CALL INVERE(R,JD,3,D,DET)
  205.          IF(IOTENS.EQ.3) THEN
  206.             R(1,1)=DET
  207.          ELSE
  208.             IF (DET.EQ.XZERO) THEN
  209. * Matrice singuliere. Numero de ligne =%i1
  210.                INTERR(1)=0
  211.                CALL ERREUR(49)
  212. *            write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  213.                RETURN
  214.             ENDIF
  215.          ENDIF
  216.       ELSEIF (IOTENS.EQ.4) THEN
  217.          IF (IKAS.EQ.1) THEN
  218.             XTRACE=A(1,1)
  219.          ELSE
  220.             XTRACE=0.D0
  221.             DO I=1,IDIM
  222.                XTRACE=XTRACE+A(I,I)
  223.             ENDDO
  224.          ENDIF
  225.          R(1,1)=XTRACE
  226.       ELSEIF (IOTENS.EQ.6) THEN
  227.          DO i=1,idim
  228.             DO j=1,idim
  229.                R(J,I)=A(J,I)
  230.            ENDDO
  231.          ENDDO
  232.       ELSEIF (IOTENS.EQ.7) THEN
  233.          DO J=1,JD
  234.             IF (D(J).GE.XPETIT) THEN
  235.                D(J)=LOG(D(J))
  236.             ELSE
  237.                REAERR(1)=D(J)
  238.                REAERR(2)=XPETIT
  239.                CALL ERREUR(1063)
  240.                RETURN
  241.             ENDIF
  242.          ENDDO
  243.       ELSEIF (IOTENS.EQ.8) THEN
  244. *         write(ioimp,*) 'TENS2 D(J) avant:',(D(idx),idx=1,jd)
  245.          DO J=1,JD
  246.             D(J)=EXP(D(J))
  247.          ENDDO
  248. *         write(ioimp,*) 'TENS2 D(J) apres:',(D(idx),idx=1,jd)
  249.       ELSEIF (IOTENS.EQ.9) THEN
  250.          DO J=1,JD
  251.             IF (ABS(D(J)).GE.XPETIT) THEN
  252.                D(J)=1.D0/D(J)
  253.             ELSE
  254.                CALL ERREUR(835)
  255.                RETURN
  256.             ENDIF
  257.          ENDDO
  258.       ELSEIF (IOTENS.EQ.10) THEN
  259.          DO J=1,JD
  260.             D(J)=ABS(D(J))
  261.          ENDDO
  262.       ELSEIF (IOTENS.EQ.11) THEN
  263.             DO J=1,JD
  264.                DO I=1,JD
  265.                R(I,J)=S(I,J)
  266.             ENDDO
  267.          ENDDO
  268.       ELSEIF (IOTENS.EQ.12) THEN
  269.          DO i=1,idim
  270.             DO j=1,idim
  271.                S(J,I)=R(J,I)
  272.            ENDDO
  273.          ENDDO
  274.       ELSEIF (IOTENS.EQ.13) THEN
  275.          DO i=1,idim
  276.             DO j=1,idim
  277.                R(I,J)=A(J,I)
  278.            ENDDO
  279.          ENDDO
  280.       ELSEIF (IOTENS.EQ.14) THEN
  281.          IF (IDIM.EQ.1) THEN
  282.             D(1)=A(1,1)
  283.          ELSEIF (IDIM.EQ.2) THEN
  284.             CALL JACOD2(A,D)
  285.          ELSEIF (IDIM.EQ.3) THEN
  286.             CALL JACOD3(A,3,D)
  287.          ELSE
  288.             write(ioimp,*) 'idim=',idim,' not implemented'
  289.             MOTERR(1:8)='TENS2'
  290.             call erreur(1039)
  291.             return
  292.          ENDIF
  293.       ELSE
  294.          write(ioimp,*) 'iotens=',iotens,' not implemented'
  295.          MOTERR(1:8)='TENS2'
  296.          call erreur(1039)
  297.          return
  298.       ENDIF
  299. * Recomposition spectrale
  300.       IF (IOTENS.GT.6.AND.IOTENS.NE.11.AND.IOTENS.LT.13) THEN
  301. c$$$            write(ioimp,*) 'D,J=',JDIM
  302. c$$$            write(ioimp,*) (D(J),J=1,JDIM)
  303. c$$$            write(ioimp,*) 'IDIM,IKAS=',IDIM,IKAS
  304.          IF (IKAS.EQ.1) THEN
  305.             R(1,1)=D(1)*S(1,1)
  306.          ELSEIF (IKAS.EQ.2) THEN
  307. c$$$            DO I=1,IDIM
  308. c$$$               R(I,1)=D(1)*S(I,1)
  309. c$$$            ENDDO
  310.             CALL ERREUR(980)
  311.             RETURN
  312.          ELSEIF (IKAS.EQ.3) THEN
  313.             DO K=1,IDIM
  314. * Tenir compte de la symétrie               DO I=1,IDIM
  315.                DO I=K,IDIM
  316.                   R(I,K)=0.D0
  317.                   DO J=1,IDIM
  318. *     faux!                        R(I,K)=R(I,K)+D(J)*S(J,I)*S(J,K)
  319.                      R(I,K)=R(I,K)+D(J)*S(I,J)*S(K,J)
  320.                   ENDDO
  321.                   IF (I.NE.K) R(K,I)=R(I,K)
  322.                ENDDO
  323.             ENDDO
  324.          ELSEIF (IKAS.EQ.5) THEN
  325.             DO K=1,IDIM
  326. * Tenir compte de la symétrie               DO I=1,IDIM
  327.                DO I=K,IDIM
  328.                   A(I,K)=0.D0
  329.                   DO J=1,IDIM
  330.                      A(I,K)=A(I,K)+D(J)*S(I,J)*S(K,J)
  331.                   ENDDO
  332.                   IF (I.NE.K) A(K,I)=A(I,K)
  333.                ENDDO
  334.             ENDDO
  335.          ELSE
  336. *     Faire la recomposition polaire après la recomposition
  337. *     spectrale dans le cas non-symetrique (mais inversible).
  338. *     Pas encore fait
  339.             CALL ERREUR(969)
  340.             RETURN
  341.          ENDIF
  342.       ENDIF
  343.       if (ldbg) then
  344.          IF (IOTENS.NE.12) THEN
  345.             write(ioimp,*) 'Output R,I,J=',IDIM,IDIM
  346.             write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
  347.          ELSE
  348.             write(ioimp,*) 'Output A,I,J=',IDIM,IDIM
  349.             write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
  350.          ENDIF
  351.          IF (IOTENS.EQ.11) THEN
  352.             write(ioimp,*) 'Output D,J=',IDIM
  353.             write(ioimp,*) (D(J),J=1,IDIM)
  354.          ENDIF
  355.       endif
  356. *
  357. * Normal termination
  358. *
  359.       RETURN
  360. *
  361. * Format handling
  362. *
  363. *
  364. * Error handling
  365. *
  366. *
  367. * End of subroutine TENS2
  368. *
  369.       END
  370.  
  371.  

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