Sources Esope | Cast3M

vectp
C VECTP     SOURCE    PV        22/04/15    13:20:15     11344          
      subroutine vectp(aa,vp,x,n)
*  recherche vecteur propre associe a une valeur propre (matrice 3x3)
*  n : multiplicite de la valeur propre
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCREEL
      dimension aa(3,3),x(*)
      a=aa(1,1)-vp
      b=aa(1,2)
      c=aa(1,3)
      d=aa(2,2)-vp
      e=aa(2,3)
      f=aa(3,3)-vp
      if (n.eq.1) then
      det1=d*f-e*e
      det2=a*f-c*c
      det3=a*d-b*b
      adet1=abs(det1)
      adet2=abs(det2)
      adet3=abs(det3)
      if (adet1.gt.xpetit) then
        det1i=1.d0/det1
      else
        det1i=1.d0
      endif
      if (adet2.gt.xpetit) then
        det2i=1.d0/det2
      else
        det2i=1.d0
      endif
      if (adet3.gt.xpetit) then
        det3i=1.d0/det3
      else
        det3i=1.d0
      endif
      if (adet3.ge.adet1.and.adet3.ge.adet2) then
        if (adet3.gt.xpetit) then
          x(1)=(b*e-c*d)*det3i
          x(2)=(b*c-a*e)*det3i
        else
          x(1)=0.D0
          x(2)=0.D0
        endif
          x(3)=+1.d0
      elseif (adet1.ge.adet2.and.adet1.ge.adet3) then
        if (adet1.gt.0.D0) then
         x(2)=(c*e-b*f)*det1i
         x(3)=(b*e-c*d)*det1i
        else
          x(2)=0.D0
          x(3)=0.D0
        endif
          x(1)=+1.d0
 
      elseif (adet2.ge.adet1.and.adet2.ge.adet3) then
        if (adet2.gt.0.D0) then
         x(1)=(c*e-b*f)*det2i
         x(3)=(b*c-a*e)*det2i
        else
          x(1)=0.D0
          x(3)=0.D0
        endif
          x(2)=+1.d0
      endif
      xn=sqrt(x(1)**2+x(2)**2+x(3)**2)
      x(1)=x(1)/xn
      x(2)=x(2)/xn
      x(3)=x(3)/xn
      elseif (n.eq.2) then
       if (abs(a).ge.abs(d).and.abs(a).ge.abs(f)) then
       x(1)=-b/a
       x(2)=1
       x(3)=0
       x(4)=-c/a
       x(5)=0
       x(6)=1
       elseif (abs(d).ge.abs(a).and.abs(d).ge.abs(f)) then
       x(1)=1
       x(2)=-b/d
       x(3)=0
       x(4)=0
       x(5)=-e/d
       x(6)=1
       elseif (abs(f).ge.abs(a).and.abs(f).ge.abs(d)) then
       x(1)=1
       x(2)=0
       x(3)=-c/f
       x(4)=0
       x(5)=1
       x(6)=-e/f
       endif
       xn=sqrt(x(1)**2+x(2)**2+x(3)**2)
       x(1)=x(1)/xn
       x(2)=x(2)/xn
       x(3)=x(3)/xn
       xsc=x(1)*x(4)+x(2)*x(5)+x(3)*x(6)
       x(4)=x(4)-xsc*x(1)
       x(5)=x(5)-xsc*x(2)
       x(6)=x(6)-xsc*x(3)
       xn=sqrt(x(4)**2+x(5)**2+x(6)**2)
       x(4)=x(4)/xn
       x(5)=x(5)/xn
       x(6)=x(6)/xn
      elseif (n.eq.3) then
       x(1)=1
       x(2)=0
       x(3)=0
       x(4)=0
       x(5)=1
       x(6)=0
       x(7)=0
       x(8)=0
       x(9)=1
      endif
      end