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

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