Sources Esope | Cast3M

jacod3
C JACOD3    SOURCE    GOUNAND   25/10/27    21:15:04     12387          
      subroutine jacod3(a,idima,d)
C======================================================================
C     OBJET
C     -----
C     DIAGONALISATION D UNE MATRICE 3*3 SYMETRIQUE
C     sans sortir les vecteurs propres
C
C     2025/10 Gounand : on reprend la methode robuste et precise de
C                       jacob3.eso
C
C     ENTREES
C     -------
C     A(3,3) = MATRICE SYMETRIQUE
C     IDIM   = 2 OU 3  SI 2 ON NE S OCCUPE QUE DE A(2,2)
C                      SI 3                    DE A(3,3)
C     SORTIES
C     -------
C     D(3)   = VALEURS PROPRES ORDONNEES D(1)>D(2)>D(3)
C
C===============================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
      dimension a(3,3),d(3),s(3,3),as(3,3)
*
* Impressions
      xpet=xpetit*10.D0
      if (idima.ne.3) then
        call jacod2(a,d)
        return
      endif
*
*     SCALING DE A
*
      amax=0.d0
      do i=1,idima
         do j=1,idima
            amax=max(amax,abs(a(j,i)))
         enddo
      enddo
* Matrice nulle
      if (amax.lt.xpet) then
         do i=1,idima
            d(i)=0.D0
         enddo
         return
      endif
C
      amax1=1.D0/amax
      do i=1,idima
         do j=1,idima
            as(j,i)=a(j,i)*amax1
         enddo
      enddo
*
*     CALCUL DES INVARIANTS
*
* On applique la methode de Harari et Albocher IJNME'2023
* pour calculer les invariants et les valeurs propres
* Invariants J1=TRACE, J2, J3=DET
      xj1=as(1,1)+as(2,2)+as(3,3)
      do i=1,idima
         do j=1,idima
            if (i.eq.j) then
               s(j,i)=as(j,i)-xj1/3
            else
               s(j,i)=as(j,i)
            endif
         enddo
      enddo
*     Dans la methode, on veut les invariants J2(S) et J3(S)
*     où S est le deviateur de A (A = S + J1(A)/3 I)
*     On a J1(S)=0 et lambda_A = lambda_S + J1(A)/3
*
*     A et S ont les memes vecteurs propres
*
*     On exprime les invariants de S de maniere stable en fonction des
*     termes de A en particulier les differences diagonales
*     En particulier J2 et J4 sont exprimes sous forme de sommes de
*     carres donc positif. On peut donc en prendre les racines carrees
*     sans fremir
      d12=as(1,1)-as(2,2)
      d23=as(2,2)-as(3,3)
      d31=as(3,3)-as(1,1)
* Expression de J2(S)
      xj2d=(d12**2+d23**2+d31**2)/6
      xj2o=as(1,2)**2+as(2,3)**2+as(3,1)**2
      xj2=xj2d+xj2o
* Expression du determinant J3(S)
      xj3d=(d12-d31)*(d23-d12)*(d31-d23)/27
      xj3m=(d12-d31)*as(2,3)**2+(d23-d12)*as(3,1)**2
     $     +(d31-d23)*as(1,2)**2
      xj3m=-xj3m/3
      xj3o=2.d0*as(1,2)*as(2,3)*as(3,1)
      xj3=xj3d+xj3m+xj3o
* Expression du discriminant Delta(S)=J4(S)=4 J2(S)^3 - 27 J3(S)^2
      xj4x=d12*d23*d31
     $     + as(1,2)**2*d12 + as(2,3)**2*d23 + as(1,3)**2*d31
      xj4y1=as(2,3)*(2*as(2,3)**2 - as(1,3)**2 - as(1,2)**2 + 2*d12*d31)
     $     +as(1,2)*as(1,3)*(d12-d31)
      xj4y2=as(1,3)*(2*as(1,3)**2 - as(2,3)**2 - as(1,2)**2 + 2*d23*d12)
     $     +as(1,2)*as(2,3)*(d23-d12)
      xj4y3=as(1,2)*(2*as(1,2)**2 - as(2,3)**2 - as(1,3)**2 + 2*d31*d23)
     $     +as(1,3)*as(2,3)*(d31-d23)
      xj4z1=as(2,3)*(as(1,3)**2 - as(1,2)**2) + as(1,2)*as(1,3)*d23
      xj4z2=as(1,3)*(as(1,2)**2 - as(2,3)**2) + as(2,3)*as(1,2)*d31
      xj4z3=as(1,2)*(as(2,3)**2 - as(1,3)**2) + as(1,3)*as(2,3)*d12
*
      xj4=xj4x**2+xj4y1**2+xj4y2**2+xj4y3**2
     $     + 15* (xj4z1**2+xj4z2**2+xj4z3**2)
*
*     CALCUL DES VALEURS PROPRES
*
      sdet=sign(1.d0,xj3)
      xsj4=sqrt(xj4)
      xsj2=sqrt(xj2)
      xdenom=2*xj2*xsj2+3*sqrt(3.d0)*sdet*xj3
      ang=2.d0/3.d0*atan2(xsj4,xdenom)
      xf1=xsj2*sdet
      xf2=cos(ang)/sqrt(3.d0)
      xf3=sin(ang)
* Si sdet=1   ds1 .GE. ds2 .GE. ds3
* Si sdet=-1  ds3 .GE. ds2 .GE. ds1
* dsi sont les valeurs propres de S
* dAi sont les valeurs propres de A
      ds1=2*xf1*xf2
      ds2=xf1*(xf3-xf2)
      ds3=xf1*(-xf3-xf2)
      xj13=xj1/3
      das1=ds1+xj13
      das2=ds2+xj13
      das3=ds3+xj13
      d(2)=das2*amax
      if (sdet.gt.0) then
         d(1)=das1*amax
         d(3)=das3*amax
      else
         d(3)=das1*amax
         d(1)=das3*amax
      endif
      return
      end