Sources Esope | Cast3M

vectp3
C VECTP3    SOURCE    GOUNAND   25/10/23    21:15:11     12385          
      subroutine vectp3(s,d,igr,x)
*
*    RECHERCHE des 2emes et 3emes VECTEURS PROPRES
*    les valeurs propres doivent être numeriquement distinctes
*
*     X en entrée est une matrice orthonormee :
*        x(j,igr)           : 1er vecteur propre vp(1) associe a d(igr)
*        x(j,i) ; i.ne.igr    2 vecteurs orthonormé orthogonaux a vp(1)
*
*     S est une matrice symetrique
*     igr vaut 1 ou 3 et correspond a la valeur propre isolee
*     d(i) : valeurs propres
 
*     Recherche un vecteur propre x(j,2) associe à la valeur propre
*     simple d(2) (important) 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)
*
*
*     X en sortie contient les vecteurs propres correspondant a chaque
*     valeur propre
*
*     On utilise la méthode de deflation de Scherzinger et Dohrmann
*     CMAME'08
*
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
      dimension s(3,3),d(3),x(3,3)
      logical ldbg
      dimension sr(3,3),u(3,3),xnu2(3)
*
* Impressions
      ldbg=.false.
* igr vaut 1 ou 3 quand ipe vaut 3 ou 1
      ipe=4-igr
      vp=d(2)
*
      do i=1,3
         do j=1,3
            sr(j,i)=s(j,i)
         enddo
      enddo
      do i=1,3
         sr(i,i)=sr(i,i)-vp
      enddo
*
*     Recherche du plus grand vecteur de l'image de Sred - d(2) I par
*     les vecteurs orthogonaux à vp1 :
*     w1
*
      if (ldbg) then
         write(6,*) 'Matrice S-d(2) I='
         do ii=1,3
            write (6,*) (sr(ii,jj),jj=1,3)
         enddo
*inu         do j=1,3
*inu            X1=SR(1,1)*x(1,j)+SR(1,2)*x(2,j)+SR(1,3)*x(3,j)
*inu            X2=SR(2,1)*x(1,j)+SR(2,2)*x(2,j)+SR(2,3)*x(3,j)
*inu            X3=SR(3,1)*x(1,j)+SR(3,2)*x(2,j)+SR(3,3)*x(3,j)
*inu            write(6,*) 'vpj d(igr) j=',j,x1,x2,x3
*inu         enddo
      endif
*
*     Calcul des images i par S-d(2)I des deux vecteurs de X orthogonaux a
*     VP1 = U (j,i) # et de leur norme xnu(i)
*
      do i=1,3
         do j=1,3
            u(j,i)=0.d0
            if (i.ne.igr) then
               do k=1,3
                  u(j,i)=u(j,i)+sr(j,k)*x(k,i)
               enddo
            endif
         enddo
      enddo
      if (ldbg) then
         write(6,*) 'Matrice U des images par S-d(2)I'
         do ii=1,3
            write (6,*) (u(ii,jj),jj=1,3)
         enddo
      endif
*
* Normalement, les deux vecteurs dans u sont orthogonaux a vp(1)
* Toutefois, lorsque les deux valeurs propres hors d(igr) sont proches,
* on peut perdre en precision sur l'orthogonalite. Donc on
* reorthogonalise systematiquement...
*
      do i=1,3
         if (i.ne.igr) then
            xpsca=xzero
            do j=1,3
               xpsca=xpsca+u(j,i)*x(j,igr)
            enddo
            do j=1,3
               u(j,i)=u(j,i)-xpsca*x(j,igr)
            enddo
         endif
      enddo
      if (ldbg) then
         write(6,*) 'Matrice U des images reorthog vs vp(1)'
         do ii=1,3
            write (6,*) (u(ii,jj),jj=1,3)
         enddo
      endif
* Calcul des normes
      imax=0
      xmax=0.d0
      do i=1,3
         xnu2(i)=xzero
         if (i.ne.igr) then
            do j=1,3
               xnu2(i)=xnu2(i)+u(j,i)**2
            enddo
            if (xnu2(i).gt.xmax) then
               xmax=xnu2(i)
               imax=i
            endif
         endif
      enddo
*
      igr2=imax
      if (ldbg) then
         write(6,*)  'Normes des colonnes de U reortho'
         write(6,*) (xnu2(jj),jj=1,3)
         write(6,*)  'indice du max des normes, igr2=',igr2
      endif
* On garde le plus grand et on le norme
* Normalement pas possible si hypothèses vérifiées (valeurs propres distinctes)
      if (igr2.eq.0) then
         do i=1,3
            if (i.ne.igr) then
               igr2=i
               goto 10
            endif
         enddo
 10      continue
         do i=1,3
            if (i.eq.igr2) then
               u(i,igr2)=1
            else
               u(i,igr2)=0
            endif
         enddo
      else
         xnorm=sqrt(xnu2(igr2))
         do i=1,3
            u(i,igr2)=u(i,igr2)/xnorm
         enddo
      endif
*
* Cela permet de construire notre 2eme et 3eme vecteur propre :
*   v2 = w1 \pvect v1
*   v3 = v1 \pvect v2
 
      x(1,2)=u(2,igr2)*x(3,igr)-u(3,igr2)*x(2,igr)
      x(2,2)=u(3,igr2)*x(1,igr)-u(1,igr2)*x(3,igr)
      x(3,2)=u(1,igr2)*x(2,igr)-u(2,igr2)*x(1,igr)
      x(1,ipe)=x(2,igr)*x(3,2)-x(3,igr)*x(2,2)
      x(2,ipe)=x(3,igr)*x(1,2)-x(1,igr)*x(3,2)
      x(3,ipe)=x(1,igr)*x(2,2)-x(2,igr)*x(1,2)
      if (ipe.eq.1) then
         x(1,ipe)=-x(1,ipe)
         x(2,ipe)=-x(2,ipe)
         x(3,ipe)=-x(3,ipe)
      endif
      return
      end