Sources Esope | Cast3M

jacob2
C JACOB2    SOURCE    GOUNAND   25/10/23    21:15:04     12385          
      SUBROUTINE JACOB2(A,D,S)
C=====================================================================
C
C       DIAGONALISATION DE A      SYMETRIQUE
C   ENTREEES
C        A(3,3) = MATRICE A DIAGONALISER  A(1,3)=A(2,3)=0.
C                                         A(3,1)=A(3,2)=0.
C   SORTIES
C        D(3)   =  LES 3 VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3)
C        S(3,3) =  VECTEURS PROPRES    S(I,3)=0. 0. 1.
C
C        RECUPERATION INCA FEVRIER 85    EBERSOLT
C====================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
C
      DIMENSION A(3,*),D(*),S(3,*)
      LOGICAL lverif,ldbg
      DIMENSION r(3,3),xnr(3)
C
C Valeur particuliere : SQRT(2.)*0.5 = 1/SQRT(2.)
      PARAMETER ( X1sr2 = 0.70710678118654752440084436210484904D0 )
* Verifications
      lverif=.false.
* Impressions
      ldbg=.false.
      xpre=xzprec*3.D0
      xpet=xpetit*10.D0
      xprev=xzprec*10.D0
*
      D(3)=A(3,3)
      S(3,1)=0.D0
      S(3,2)=0.D0
      S(1,3)=0.D0
      S(2,3)=0.D0
      S(3,3)=1.D0
* Scaling de A
      amax=0.d0
      do i=1,2
         do j=1,2
            amax=max(amax,abs(a(j,i)))
         enddo
      enddo
* Matrice nulle
      if (amax.lt.xpet) then
         D(1)=XZERO
         D(2)=XZERO
         S(1,1)=1.D0
         S(2,1)=XZERO
         S(1,2)=XZERO
         S(2,2)=1.D0
         return
      endif
C
      amax1=1.D0/amax
      do i=1,2
         do j=1,2
            R(j,i)=a(j,i)*amax1
         enddo
      enddo
      X1  =2.D0*R(1,2)
      X2  =R(1,1)-R(2,2)
      X3  =SQRT(X2*X2+X1*X1)
      D(1)=0.5D0*(R(1,1)+R(2,2)+X3)
      D(2)=D(1)-X3
*     Raccourci lorsque les deux valeurs propres sont egales
*
      xsca=max(xpre,xpet)
      if (ldbg) then
         write(6,*) 'xpre,xpet,xsca,x3=',xpre,xpet,xsca,x3
      endif
*
      if (x3.lt.xsca) then
         if (ldbg) write(6,*) '2 valeurs propres egales detectees'
         do i=1,2
            do j=1,2
               if (i.eq.j) then
                  s(j,i)=1.d0
               else
                  s(j,i)=0.d0
               endif
            enddo
         enddo
         goto 1000
      endif
C
C Appliquer la methode de Scherzinger et Dohrmann CMAME'08
* Recherche du plus grand vecteur de l'image de A - d(1) I :
C
      do i=1,2
         R(i,i)=R(i,i)-d(1)
      enddo
      imax=0
      xmax=0.d0
      do i=1,2
         xnr(i)=0.d0
         do j=1,2
            xnr(i)=xnr(i)+r(j,i)**2
         enddo
         if (xnr(i).gt.xmax) then
            xmax=xnr(i)
            imax=i
         endif
      enddo
      igr=imax
* Normalement pas possible si valeurs propres non egales
      if (igr.eq.0) then
         igr=1
         r(1,1)=1.d0
         r(2,1)=0.d0
      else
         xnorm=sqrt(xnr(igr))
         do i=1,2
            r(i,igr)=r(i,igr)/xnorm
         enddo
      endif
      S(1,2)= r(1,igr)
      S(2,2)= r(2,igr)
      S(1,1)= S(2,2)
      S(2,1)=-S(1,2)
 
      if (.false.) then
*       Methode ancienne
         R11=A(1,1)-D(1)
         R21=A(1,2)
         XNR1=SQRT(R11*R11+R21*R21)
         R12=A(2,1)
         R22=A(2,2)-D(1)
         XNR2=SQRT(R12*R12+R22*R22)
         IF (XNR1.GE.XNR2) THEN
            I=1
            T1=R11
            T2=R21
            XT=XNR1
         ELSE
            I=2
            T1=R12
            T2=R22
            XT=XNR2
         ENDIF
C
         XSCA=MAX(MAX(ABS(D(1)),ABS(D(2)))*XZPREC,XPETIT)*2
         IF (XT.LE.XSCA) THEN
*         write(6,*) 'Cas 0'
            X5    = X1sr2
            S(1,1)=X5
            S(2,1)=SIGN(X5,X1)
            S(1,2)=-S(2,1)
            S(2,2)=X5
         ELSE
*
            S1=T1/XT
            S2=T2/XT
*
*       IF (I.EQ.1) THEN
*            write(6,*) 'Cas 1'
            S(1,1)= S2
            S(2,1)=-S1
            S(1,2)= S1
            S(2,2)= S2
*         ELSE
*            write(6,*) 'Cas 2'
*            S(1,1)= S1
*            S(2,1)= S2
*            S(1,2)=-S2
*            S(2,2)= S1
*         ENDIF
         ENDIF
      endif
 1000 continue
* Scaling des valeurs propres
      d(1)=d(1)*amax
      d(2)=d(2)*amax
      if (lverif) then
* Orthonormalite
         call vmorth(s,2,3,xpre,r)
         if (ierr.ne.0) then
            write(6,*) '!!! X non orthonormale'
            if (ldbg) then
               write(6,*) 'Matrice des VPs s='
               do ii=1,2
                  write (6,*) (s(ii,jj),jj=1,2)
               enddo
               write(6,*) 'Matrice des Pscals R='
               do ii=1,2
                  write (6,*) (R(ii,jj),jj=1,2)
               enddo
            endif
            goto 9999
         endif
* Vecteur propre
         xprevp=amax*xprev
         do i=1,2
*     vp(i) est-il bien le vp associe a la valeur propre d(i)
*     xnr espace de stockage ici
            call vvectp(a,2,3,d(i),s(1,i),xprevp,xnr)
            if (ierr.ne.0) then
               write(6,*) '!!!   vp(i) nest pas vecteur propre i=',i
               if (ldbg) then
                  write(6,*) ' v= (S - D(i)) s(j,i)  ='
                  write(6,*) (xnr(jj),jj=1,2)
               endif
               goto 9999
            endif
         enddo
      endif
      RETURN
*
* Error handling
*
 
 9999 CONTINUE
      MOTERR(1:8)='JACOB2'
      call erreur(1127)
      RETURN
      END