C JACOB3    SOURCE    GOUNAND   25/10/23    21:15:05     12385          
      subroutine jacob3(a,idima,d,x)
C======================================================================
C     OBJET
C     -----
C     DIAGONALISATION D UNE MATRICE 3*3 SYMETRIQUE
C
C     ENTREES
C     -------
C     A(3,3) = MATRICE SYMETRIQUE
C     IDIMA   = 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     X(3,3) = VECTEURS PROPRES ( X(IP,2) EST LE VECTEUR
C                                 ASSOCIE A D(2) )
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),x(3,3),s(3,3),as(3,3)
      logical ldbg,lverif
      dimension r(3,3),v(3)
*
* Verifications
      lverif=.false.
* Impressions
      ldbg=.false.
      xpre=xzprec*10.D0
      xpet=xpetit*10.D0
      xprev=xzprec*10.D0
*      xpre=xzprec*30.D0
*      xpet=xpetit*30.D0
      if (idima.ne.3) then
        call jacob2(a,d,x)
        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
            do j=1,idima
               if (i.eq.j) then
                  x(j,i)=1.d0
               else
                  x(j,i)=0.d0
               endif
            enddo
         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
      if (ldbg) then
         write(6,*) 'Matrice S='
         do ii=1,idima
            write (6,*) (S(ii,jj),jj=1,idima)
         enddo
      endif
*     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)
      if (ldbg) then
* Impression des invariants de S
         c2=-s(1,1)-s(2,2)-s(3,3)
         c1= (s(1,1)*s(2,2)+s(2,2)*s(3,3)+s(3,3)*s(1,1))
     &        - s(1,3)**2 - s(1,2)**2 - s(2,3)**2
         c0=-2.d0*s(1,2)*s(1,3)*s(2,3) + s(1,1)*s(2,3)**2
     &        + s(2,2)*s(1,3)**2 + s(3,3)*s(1,2)**2
     &        - s(1,1)*s(2,2)*s(3,3)
         xj1b=-c2+xj1
         xj2b=-c1
         xj3b=-c0
         xj4b=4*xj2b**3-27*xj3b**2
         write(6,*) 'Invariants du deviateur S : new old  abs(old-new)'
         write(6,*) ' trace= ',xj1,xj1b,abs(xj1b-xj1)
         write(6,*) ' J2=    ',xj2,xj2b,abs(xj2b-xj2)
         write(6,*) ' det=   ',xj3,xj3b,abs(xj3b-xj3)
         write(6,*) ' discr= ',xj4,xj4b,abs(xj4b-xj4)
      endif
*
*     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
* Impression des valeurs propres de A par les methodes nouvelles et anciennes
      if (ldbg) then
         d(2)=das2
         if (sdet.gt.0) then
            d(1)=das1
            d(3)=das3
         else
            d(3)=das1
            d(1)=das3
         endif
*   Methode ancienne
         c2=-as(1,1)-as(2,2)-as(3,3)
         c1= (as(1,1)*as(2,2)+as(2,2)*as(3,3)+as(3,3)*as(1,1))
     &        - as(1,3)**2 - as(1,2)**2 - as(2,3)**2
         c0=-2.d0*as(1,2)*as(1,3)*as(2,3) + as(1,1)*as(2,3)**2
     &        + as(2,2)*as(1,3)**2 + as(3,3)*as(1,2)**2
     &        - as(1,1)*as(2,2)*as(3,3)
         write(6,*) 'c0,c1,c2=',c0,c1,c2
         call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3)
         write(6,*) 'Som vp=',d1+d2+d3
         write(6,*) 'Prod vp=',d1*d2*d3
         if (d1.lt.d2) then
            dd=d1
            d1=d2
            d2=dd
         endif
         if (d1.lt.d3) then
            dd=d1
            d1=d3
            d3=dd
         endif
         if (d2.lt.d3) then
            dd=d2
            d2=d3
            d3=dd
         endif
         db1=d1
         db2=d2
         db3=d3
         write(6,*) 'Val. propres de A scalee : new old  abs(old-new)'
         write(6,*) ' vp1= ',d(1),db1,abs(db1-d(1))
         write(6,*) ' vp2= ',d(2),db2,abs(db2-d(2))
         write(6,*) ' vp3= ',d(3),db3,abs(db3-d(3))
      endif
*
*     Faire le calcul des vecteurs propres avec S plutôt qu'avec A
*     (moins d'erreur de cancellation)
*     On utilise la méthode de deflation de Scherzinger et Dohrmann
*     CMAME'08 car sinon, on n'a pas reussi a avoir une bonne
*     orthogonalite des vecteurs propres lorsque les valeurs propres
*     sont proches
*
      d(2)=ds2
* Si sdet est positif, d(1) est la valeur propre la plus isolee
* Si sdet est negatif, d(3) est la valeur propre la plus isolee
* On recherche d'abord le vecteur propre associe a cette valeur
*     propre
      if (sdet.gt.0) then
         igr=1
      else
         igr=3
      endif
* igr vaut 1 ou 3 quand ipe vaut 3 ou 1
      ipe=4-igr
      d(igr)=ds1
      d(ipe)=ds3
      if (ldbg) then
         write(6,*) 'Val. propres de S'
         write(6,*) ' vp1= ',d(1)
         write(6,*) ' vp2= ',d(2)
         write(6,*) ' vp3= ',d(3)
      endif
*
*     CALCUL DES VECTEURS PROPRES
*
*     Raccourci lorsque les trois valeurs propres sont egales
*     Dans ce cas S=0 et J2(S)=0
*
      xsca=max(xpre,xpet)
      if (ldbg) then
         write(6,*) 'xpre,xpet,xsca,xsj2=',xpre,xpet,xsca,xsj2
      endif
*
      if (xsj2.lt.xsca) then
         if (ldbg) write(6,*) '3 valeurs propres egales detectees'
         do i=1,idima
            do j=1,idima
               if (i.eq.j) then
                  x(j,i)=1.d0
               else
                  x(j,i)=0.d0
               endif
            enddo
         enddo
         goto 1000
      endif
*
*     Faire le calcul des vecteurs propres avec S plutôt qu'avec A
*     (moins d'erreur de cancellation)
*     On utilise la méthode de deflation de Scherzinger et Dohrmann
*     CMAME'08 car sinon, on n'a pas reussi a avoir une bonne
*     orthogonalite des vecteurs propres lorsque les valeurs propres
*     sont proches
*
*     RECHERCHE DU 1er VECTEUR PROPRE associe a la valeur propre isolee
*     d(igr)
*
      call vectp2(s,d,igr,x)
      if (lverif) then
*     Verif
         call vmorth(x,idima,3,xpre,r)
         if (ierr.ne.0) then
            write(6,*) '!!! X non orthonormale'
            if (ldbg) then
               write(6,*) 'Matrice des VPs X='
               do ii=1,idima
                  write (6,*) (X(ii,jj),jj=1,idima)
               enddo
               write(6,*) 'Matrice des Pscals R='
               do ii=1,idima
                  write (6,*) (R(ii,jj),jj=1,idima)
               enddo
            endif
            goto 9999
         endif
*     vp(igr) est-il bien le vp associe a la valeur propre d(igr)
         call vvectp(s,idima,3,d(igr),x(1,igr),xprev,v)
         if (ierr.ne.0) then
            write(6,*) '!!!   vp(igr) nest pas vecteur propre '
            if (ldbg) then
               write(6,*) ' v= (S - D(igr)) x(j,igr)  ='
               write(6,*) (v(jj),jj=1,idima)
            endif
            goto 9999
         endif
      endif
*
*     FIN DE RECHERCHE DU 1er VECTEUR PROPRE
*
*
*     Raccourci lorsque les deux autres valeurs propres sont egales
*     Dans ce cas le discriminant est nul J4(S)=0
*     Le developpement limite de lambda_2 et lambda_3 en J4 tendant vers
*     0 est :
*     lambda_2,3 \propto sqrt(J2) ( 1 + sqrt(J4) / (J2^3/2)  )
*

      xsca=max((xsj2**3)*xpre,xpet)
      if (ldbg) then
         write(6,*) 'xsj2,xsj23,xpre,xsca,xsj4=',xsj2,xsj2**3,xpre,xsca
     $        ,xsj4
      endif
      if (xsj4.lt.xsca) then
         if (ldbg) write(6,*) '2 valeurs propres egales detectees'
*     La matrice X contient une base donc on peut quitter prematurement
         goto 1000
      endif
*
*    RECHERCHE des 2emes et 3emes VECTEURS PROPRES
*    les valeurs propres sont normalement numeriquement distinctes
*
*     Recherche un vecteur propre x(j,2) associe à la valeur propre
*     simple d(2)  d'une matrice 3x3 S
*
*     L'image de S-d(2)I est de dimension 2 car d(2) est valeur
*     propre simple. On cherche uniquement la partie de l'image
*     orthogonale à vp(igr).
*     Pour ce faire, on prend la plus grande des images des deux
*     vecteurs de la base X différents de vp(igr). Ces deux images sont
*     colinéaires car à la fois orthogonales à vp(igr) et vp(2) : w1
*
*     On choisit ensuite :  vp(2) = w1 \pvec vp(1)
*                  puis  :  vp(3) = vp(1) \pvec vp(2)
*
*
      call vectp3(s,d,igr,x)
*
* Verif
      if (lverif) then
* Orthonormalite des vps
         call vmorth(x,idima,3,xpre,r)
         if (ierr.ne.0) then
            write(6,*) '!!! X non orthonormale'
            if (ldbg) then
               write(6,*) 'Matrice des VPs X='
               do ii=1,idima
                  write (6,*) (X(ii,jj),jj=1,idima)
               enddo
               write(6,*) 'Matrice des Pscals R='
               do ii=1,idima
                  write (6,*) (R(ii,jj),jj=1,idima)
               enddo
            endif
            goto 9999
         endif
         do i=1,idima
*     vp(i) est-il bien le vp associe a la valeur propre d(i)
            call vvectp(s,idima,3,d(i),x(1,i),xprev,v)
            if (ierr.ne.0) then
               write(6,*) '!!!   vp(i) nest pas vecteur propre i=',i
               if (ldbg) then
                  write(6,*) ' v= (S - D(i)) x(j,i)  ='
                  write(6,*) (v(jj),jj=1,idima)
               endif
               goto 9999
            endif
         enddo
      endif
*
      if (.false.) then
* Methode ancienne
         dmaxa=max(abs(d(1)),abs(d(2)),abs(d(3)))
         deps=max(dmaxa*xpre,xpetit*10.d0)
         if (abs(d(1)-d(2)).le.deps) then
*     valeur propre double
            if (abs(d(2)-d(3)).le.deps) then
*     valeur propre triple
               if (ldbg) write(6,*) 'JACOB3 Valeur propre triple'
               call vectp(s,d(1),x(1,1),3)
            else
               if (ldbg) write(6,*) 'JACOB3 Valeur propre double 1'
               call vectp(a,d(1),x(1,1),2)
               x(1,3)=x(2,1)*x(3,2)-x(3,1)*x(2,2)
               x(2,3)=x(3,1)*x(1,2)-x(1,1)*x(3,2)
               x(3,3)=x(1,1)*x(2,2)-x(2,1)*x(1,2)
*old     call vectp(a,d(3),x(1,3),1)
            endif
         else
            if (abs(d(2)-d(3)).le.deps) then
               if (ldbg) write(6,*) 'JACOB3 Valeur propre double 2'
*       valeur propre double
               call vectp(a,d(2),x(1,2),2)
               x(1,1)=x(2,2)*x(3,3)-x(3,2)*x(2,3)
               x(2,1)=x(3,2)*x(1,3)-x(1,2)*x(3,3)
               x(3,1)=x(1,2)*x(2,3)-x(2,2)*x(1,3)
*old            call vectp(a,d(1),x(1,1),1)
            else
*       cas normal
               if (ldbg) write(6,*) 'JACOB3 Cas normal'
               call vectp(s,d(1),x(1,1),1)
               call vectp(s,d(2),x(1,2),1)
               call vectp(s,d(3),x(1,3),1)
            endif
         endif
      endif
*
* On remet dans D les valeurs propres de A a la place de celles de S
*
 1000 continue
      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
*
* Error handling
*

 9999 CONTINUE
      MOTERR(1:8)='JACOB3'
      call erreur(1127)
      return
      end