Sources Esope | Cast3M

vectp2
C VECTP2    SOURCE    GOUNAND   25/10/23    21:15:10     12385          
      subroutine vectp2(s,d,igr,x)
*
*     Recherche un vecteur propre x(j,igr) associe a une valeur propre
*     simple d(igr) (important) d'une matrice 3x3 S
*     Les deux autres vecteurs de x sont tels que les trois vecteurs
*     forment une base orthonormale
*
*     On utilise la méthode de deflation de Scherzinger et Dohrmann
*     CMAME'08
*     L'image de S-d(igr)I est de dimension 2 car d(igr) est valeur
*     propre simple. On determine d'abord une base B de cet espace.
*     Le complementaire est l'espace de dimension 1 associe au vecteur
*     propre x(j,igr) obtenu par produit vectoriel des 2 vecteurs de la
*     base B
*
*
      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 r(3,3),xnr2(3)
*
* Impressions
      ldbg=.false.
*
      vp=d(igr)
C Appliquer la methode ds Scherzinger et Dohrmann CMAME'08
      r(1,1)=s(1,1)-vp
      r(2,1)=s(1,2)
      r(3,1)=s(1,3)
      r(1,2)=r(2,1)
      r(2,2)=s(2,2)-vp
      r(3,2)=s(2,3)
      r(1,3)=r(3,1)
      r(2,3)=r(3,2)
      r(3,3)=s(3,3)-vp
      if (ldbg) then
         write(6,*) 'igr,d(igr)=',igr,d(igr)
         write(6,*) 'Matrice R = S-d(igr) I='
         do ii=1,idim
            write (6,*) (r(ii,jj),jj=1,idim)
         enddo
      endif
*
*     Les colonnes de R sont les vecteurs obtenus par produit de R avec
*     les vecteurs de la base canonique
*     Celle de norme la plus grande sera notre 1er vecteur vb1 de la
*     base B
*
      imax=0
      xmax=0.d0
      do i=1,3
         xnr2(i)=xzero
         do j=1,3
            xnr2(i)=xnr2(i)+r(j,i)**2
         enddo
         if (xnr2(i).gt.xmax) then
            xmax=xnr2(i)
            imax=i
         endif
      enddo
      jgr=imax
      if (ldbg) then
         write(6,*)  'Normes des colonnes de R'
         write(6,*) (xnr2(jj),jj=1,idim)
         write(6,*)  'indice du max des normes, jgr=',jgr
      endif
* Normalement pas possible si hypothèses vérifiées
      if (jgr.eq.0) then
         jgr=1
         do i=1,3
            if (i.eq.jgr) then
               r(i,jgr)=1
            else
               r(i,jgr)=0
            endif
         enddo
      endif
* On normalise vb1
      xnorm=sqrt(xnr2(jgr))
      do i=1,3
         r(i,jgr)=r(i,jgr)/xnorm
      enddo
* On orthogonalise les deux autres vecteurs par rapport a vb1
      do j=1,3
         if (j.ne.jgr) then
            xpsca=xzero
            do i=1,3
               xpsca=xpsca+r(i,jgr)*r(i,j)
            enddo
            do i=1,3
               r(i,j)=r(i,j)-xpsca*r(i,jgr)
            enddo
         endif
      enddo
      if (ldbg) then
         write(6,*) 'Matrice R orthogonalise vs colonne jgr=',jgr
         do ii=1,idim
            write (6,*) (r(ii,jj),jj=1,idim)
         enddo
      endif
* On calcule leur norme au carre
      xmax=xzero
      imax=0
      do i=1,3
         if (i.ne.jgr) then
            xnr2(i)=xzero
            do j=1,3
               xnr2(i)=xnr2(i)+r(j,i)**2
            enddo
            if (xnr2(i).gt.xmax) then
               xmax=xnr2(i)
               imax=i
            endif
         endif
      enddo
      jgr2=imax
      if (ldbg) then
         write(6,*)  'Normes des colonnes de R ortho'
         write(6,*) (xnr2(jj),jj=1,idim)
         write(6,*)  'indice du max des normes, jgr2=',jgr2
      endif
*  Le plus grand est notre vecteur vb2, qu'on normalise
      if (jgr2.eq.0) then
         do i=1,3
            if (i.ne.jgr) then
               jgr2=i
               goto 10
            endif
         enddo
 10      continue
         do i=1,3
            if (i.eq.jgr2) then
               r(i,jgr2)=1
            else
               r(i,jgr2)=0
            endif
         enddo
      else
         xnorm=sqrt(xnr2(jgr2))
         do i=1,3
            r(i,jgr2)=r(i,jgr2)/xnorm
         enddo
      endif
*     Le vecteur propre associe a d(igr) est le produit vectoriel de vb1
*     et vb2
      ioth=modulo(igr,3)+1
      do i=1,3
         x(i,ioth)=r(i,jgr)
      enddo
      ioth2=modulo(igr+1,3)+1
      do i=1,3
         x(i,ioth2)=r(i,jgr2)
      enddo
      x(1,igr)=x(2,ioth)*x(3,ioth2)-x(3,ioth)*x(2,ioth2)
      x(2,igr)=x(3,ioth)*x(1,ioth2)-x(1,ioth)*x(3,ioth2)
      x(3,igr)=x(1,ioth)*x(2,ioth2)-x(2,ioth)*x(1,ioth2)
      return
      end