tens2
C TENS2 SOURCE GOUNAND 24/09/27 21:15:23 12020 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 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 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=',JDIM c$$$ write(ioimp,*) (D(J),J=1,JDIM) 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 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 XNORM=SQRT(XNORM) 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 IF(IOTENS.EQ.3) THEN ELSE * Matrice singuliere. Numero de ligne =%i1 INTERR(1)=0 * 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 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 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' 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 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 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
© Cast3M 2003 - Tous droits réservés.
Mentions légales