Sources Esope | Cast3M

tens2
C TENS2     SOURCE    GOUNAND   25/07/28    21:15:12     12338          
      SUBROUTINE TENS2(IOTENS,IKAS,A,D,R)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : TENS2
C DESCRIPTION : Opérations sur des tenseurs (unaires pour l'instant)
C               Ici on gère les opérations à l'échelle du noeud
C
C
C LANGAGE     : ESOPE
C AUTEUR      : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA)
C               mel : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          :
C APPELES (E/S)    :
C APPELES (BLAS)   :
C APPELES (CALCUL) :
C APPELE PAR       :
C***********************************************************************
C SYNTAXE GIBIANE    :
C ENTREES            :
C ENTREES/SORTIES    :
C SORTIES            :
C***********************************************************************
C VERSION    : v1, 28/08/2024, version initiale
C HISTORIQUE : v1, 28/08/2024, creation
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
*
C   ENTREEES
C        A(3,3) = MATRICE A DIAGONALISER
C   SORTIES
C        D(3)   =  VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3)
C        S(3,3) =  VECTEURS PROPRES
C        R(3,3) =  RESULTAT DES OPERATIONS TENSORIELLES
      DIMENSION A(3,3),D(3),S(3,3),R(3,3)
      DATA ((S(I,J),I=1,3),J=1,3) /9*0.D0/
*
*      DATA MOTENS/'NORM2','NORMINF','DET','TRACE','INVERSE','IDEN','LOG'
*     $     ,'EXP','INVS','ABSOLU','PRINCIPA'/
*
*
* Executable statements
*
      IF (IOTENS.GT.5) THEN
*
*     Décomposition spectrale (positive pour les cas scalaires et
*     vecteurs)
*
         IF (IKAS.EQ.1.OR.(IKAS.EQ.3.AND.IDIM.EQ.1)) THEN
* Difficile de donner un sens...
*            D(1)=ABS(A(1,1))
*            S(1,1)=SIGN(1.D0,A(1,1))
            D(1)=A(1,1)
            S(1,1)=1.D0
            JD=1
         ELSEIF (IKAS.EQ.2) THEN
* Pas de sens pour les vecteurs
c$$$            XNO=0.D0
c$$$            DO I=1,IDIM
c$$$               XNO=XNO+A(I,1)**2
c$$$            ENDDO
c$$$            XNA=SQRT(XNO)
c$$$            D(1)=XNA
c$$$            DO I=1,IDIM
c$$$               S(I,1)=A(I,1)/XNA
c$$$            ENDDO
            CALL ERREUR(980)
            RETURN
c$$$         write(ioimp,*) 'A,I,J=',IDIM,1
c$$$         write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,1)
c$$$         write(ioimp,*) 'D,J=',JDIM
c$$$         write(ioimp,*) (D(J),J=1,JDIM)
c$$$         write(ioimp,*) 'S,I,J=',IDIM,1
c$$$         write(ioimp,*) ((S(I,J),I=1,IDIM),J=1,1)
 
         ELSEIF (IKAS.EQ.3) THEN
            CALL JACOB3(A,IDIM,D,S)
            JD=IDIM
c$$$               write(ioimp,*) 'A,I,J=',IDIM,IDIM
c$$$               write(ioimp,*) ((A(I,J),I=1,IDIM),J=1,IDIM)
c$$$               write(ioimp,*) 'D,J=',IDIM
c$$$               write(ioimp,*) (D(J),J=1,IDIM)
c$$$               write(ioimp,*) 'S,I,J=',IDIM,IDIM
c$$$               write(ioimp,*) ((S(I,J),I=1,IDIM),J=1,IDIM)
c$$$*     Verif
c$$$               DO K=1,IDIM
c$$$*                  DO I=K,IDIM
c$$$                  DO I=1,IDIM
c$$$                     R(I,K)=0.D0
c$$$                     DO J=1,IDIM
c$$$*                        R(I,K)=R(I,K)+D(J)*S(J,I)*S(J,K)
c$$$                        R(I,K)=R(I,K)+D(J)*S(I,J)*S(K,J)
c$$$                     ENDDO
c$$$                  ENDDO
c$$$               ENDDO
c$$$               write(ioimp,*) 'A2,I,J=',IDIM,JDIM
c$$$               write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,JDIM)
c$$$               DO J=1,IDIM
c$$$                  DO I=1,IDIM
c$$$                     V(I)=0.D0
c$$$                     DO K=1,IDIM
c$$$                        V(I)=V(I)+A(I,K)*S(K,J)
c$$$                     ENDDO
c$$$                  ENDDO
c$$$                  write(ioimp,*) 'VJ(I)=AIK.SKJ,J',J
c$$$                  write(ioimp,*) (V(I),I=1,IDIM)
c$$$                  write(ioimp,*) 'VJ(I)=D(J).SIJ,J',J
c$$$                  write(ioimp,*) (S(I,J)*D(J),I=1,IDIM)
c$$$               ENDDO
         ELSE
*     D'abord faire la décomposition polaire avant la décomposition
*     spectrale dans le cas non-symetrique (mais inversible).
*     Pas encore fait
            CALL ERREUR(969)
            RETURN
         ENDIF
      ENDIF
*
* Les opérations tensorielles proprement dites
*
      IF (IOTENS.EQ.1) THEN
         IF (IKAS.EQ.1) THEN
            XNORM=ABS(A(1,1))
         ELSE
            IF (IKAS.EQ.2) THEN
               ID=IDIM
               JD=1
            ELSE
               ID=IDIM
               JD=IDIM
            ENDIF
            XNORM=0.D0
            DO J=1,JD
               DO I=1,ID
                  XNORM=XNORM+A(I,J)**2
               ENDDO
            ENDDO
            XNORM=SQRT(XNORM)
         ENDIF
         R(1,1)=XNORM
      ELSEIF (IOTENS.EQ.2) THEN
         IF (IKAS.EQ.1) THEN
            XNORM=ABS(A(1,1))
         ELSE
            IF (IKAS.EQ.2) THEN
               ID=IDIM
               JD=1
            ELSE
               ID=IDIM
               JD=IDIM
            ENDIF
            XNORM=0.D0
            DO J=1,JD
               DO I=1,ID
                  XNORM=max(xnorm,abs(A(I,J)))
               ENDDO
            ENDDO
         ENDIF
         R(1,1)=XNORM
      ELSEIF (IOTENS.EQ.3.OR.IOTENS.EQ.5) THEN
*
         IF (IKAS.EQ.1) THEN
            JD=1
         ELSE
            JD=IDIM
         ENDIF
         DO J=1,JD
            DO I=1,JD
               R(I,J)=A(I,J)
            ENDDO
         ENDDO
* D sert de vecteur de travail ici
         CALL INVERE(R,JD,3,D,DET)
         IF(IOTENS.EQ.3) THEN
            R(1,1)=DET
         ELSE
            IF (DET.EQ.XZERO) THEN
* Matrice singuliere. Numero de ligne =%i1
               INTERR(1)=0
               CALL ERREUR(49)
*            write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,IDIM)
               RETURN
            ENDIF
         ENDIF
      ELSEIF (IOTENS.EQ.4) THEN
         IF (IKAS.EQ.1) THEN
            XTRACE=A(1,1)
         ELSE
            XTRACE=0.D0
            DO I=1,IDIM
               XTRACE=XTRACE+A(I,I)
            ENDDO
         ENDIF
         R(1,1)=XTRACE
      ELSEIF (IOTENS.EQ.6) THEN
      ELSEIF (IOTENS.EQ.7) THEN
         DO J=1,JD
            IF (D(J).GE.XPETIT) THEN
               D(J)=LOG(D(J))
            ELSE
               REAERR(1)=D(J)
               REAERR(2)=XPETIT
               CALL ERREUR(1063)
               RETURN
            ENDIF
         ENDDO
      ELSEIF (IOTENS.EQ.8) THEN
         DO J=1,JD
            D(J)=EXP(D(J))
         ENDDO
      ELSEIF (IOTENS.EQ.9) THEN
         DO J=1,JD
            IF (ABS(D(J)).GE.XPETIT) THEN
               D(J)=1.D0/D(J)
            ELSE
               CALL ERREUR(835)
               RETURN
            ENDIF
         ENDDO
      ELSEIF (IOTENS.EQ.10) THEN
         DO J=1,JD
            D(J)=ABS(D(J))
         ENDDO
      ELSEIF (IOTENS.EQ.11) THEN
            DO J=1,JD
               DO I=1,JD
               R(I,J)=S(I,J)
            ENDDO
         ENDDO
      ELSE
         write(ioimp,*) 'iotens=',iotens,' not implemented'
         MOTERR(1:8)='TENS2'
         call erreur(1039)
         return
      ENDIF
* Recomposition spectrale
      IF (IOTENS.GT.5.AND.IOTENS.LT.11) THEN
c$$$            write(ioimp,*) 'D,J=',JDIM
c$$$            write(ioimp,*) (D(J),J=1,JDIM)
c$$$            write(ioimp,*) 'IDIM,IKAS=',IDIM,IKAS
         IF (IKAS.EQ.1.OR.(IKAS.EQ.3.AND.IDIM.EQ.1)) THEN
            R(1,1)=D(1)*S(1,1)
         ELSEIF (IKAS.EQ.2) THEN
c$$$            DO I=1,IDIM
c$$$               R(I,1)=D(1)*S(I,1)
c$$$            ENDDO
            CALL ERREUR(980)
            RETURN
         ELSEIF (IKAS.EQ.3) THEN
            DO K=1,IDIM
* Tenir compte de la symétrie               DO I=1,IDIM
               DO I=K,IDIM
                  R(I,K)=0.D0
                  DO J=1,IDIM
*faux!                        R(I,K)=R(I,K)+D(J)*S(J,I)*S(J,K)
                     R(I,K)=R(I,K)+D(J)*S(I,J)*S(K,J)
                  ENDDO
               ENDDO
            ENDDO
         ELSE
*     Faire la recomposition polaire après la recomposition
*     spectrale dans le cas non-symetrique (mais inversible).
*     Pas encore fait
            CALL ERREUR(969)
            RETURN
         ENDIF
c$$$            write(ioimp,*) 'R,I,J=',IDIM,JDIM
c$$$            write(ioimp,*) ((R(I,J),I=1,IDIM),J=1,JDIM)
      ENDIF
*
* Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
*
* End of subroutine TENS2
*
      END