1 : C JACOB4 SOURCE AM 08/12/19 21:16:48 6237 2 : subroutine jacob4(a,idim,d,x) 3 : C====================================================================== 4 : C OBJET 5 : C ----- 6 : C DIAGONALISATION D UNE MATRICE 3*3 SYMETRIQUE 7 : C 8 : C ENTREES 9 : C ------- 10 : C A(3,3) = MATRICE SYMETRIQUE 11 : C IDIM = 2 OU 3 SI 2 ON NE S OCCUPE QUE DE A(2,2) 12 : C SI 3 DE A(3,3) 13 : C SORTIES 14 : C ------- 15 : C D(3) = VALEURS PROPRES ORDONNEES D(1)>D(2)>D(3) 16 : C 17 : C S(3,3) = VECTEURS PROPRES ( S(IP,2) EST LE VECTEUR 18 : C ASSOCIE A D(2) ) 19 : C 20 : C=============================================================== 21 : IMPLICIT INTEGER(I-N) 22 : IMPLICIT REAL*8 (A-H,O-Z) 23 : dimension a(3,3),d(3),x(3,3) 24 : if (idim.ne.3) then 24 : call jacob2(a,d,x) 25 : return 26 : endif 27 : c2=-a(1,1)-a(2,2)-a(3,3) 28 : c1= (a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1)) 29 : & - a(1,3)**2 - a(1,2)**2 - a(2,3)**2 30 : c0=-2.*a(1,2)*a(1,3)*a(2,3) + a(1,1)*a(2,3)**2 31 : & + a(2,2)*a(1,3)**2 + a(3,3)*a(1,2)**2 32 : & - a(1,1)*a(2,2)*a(3,3) 32 : call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3) 33 : d(1)=max(d1,d2,d3) 34 : d(3)=min(d1,d2,d3) 35 : d(2)=d1+d2+d3-d(1)-d(3) 36 : deps=d(1)*1.e-4 37 : if (d(1)-d(2).le.deps) then 38 : * valeur propre double 39 : if (d(2)-d(3).le.deps) then 40 : * valeur propre triple 40 : call vectp(a,d(1),x(1,1),3) 41 : else 41 : call vectp(a,d(1),x(1,1),2) 41 : call vectp(a,d(3),x(1,3),1) 42 : endif 43 : else 44 : deps=d(2)*1.e-4 45 : if (d(2)-d(3).le.deps) then 46 : * valeur propre double 46 : call vectp(a,d(1),x(1,1),1) 46 : call vectp(a,d(2),x(1,2),2) 47 : else 48 : * cas normal 48 : call vectp(a,d(1),x(1,1),1) 48 : call vectp(a,d(2),x(1,2),1) 48 : call vectp(a,d(3),x(1,3),1) 49 : endif 50 : endif 51 : end 52 : * 53 : 54 :
© Cast3M 2003 - All rights reserved.
Disclaimer