Télécharger tens2.eso

Retour à la liste

Numérotation des lignes :

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

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