Télécharger jacob3.eso

Retour à la liste

Numérotation des lignes :

jacob3
  1. C JACOB3 SOURCE AM 19/12/05 21:15:02 10416
  2. subroutine jacob3(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.  
  25. if (idim.ne.3) then
  26. call jacob2(a,d,x)
  27. return
  28. endif
  29.  
  30. c2=-a(1,1)-a(2,2)-a(3,3)
  31. c1= (a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1))
  32. & - a(1,3)**2 - a(1,2)**2 - a(2,3)**2
  33. c0=-2.*a(1,2)*a(1,3)*a(2,3) + a(1,1)*a(2,3)**2
  34. & + a(2,2)*a(1,3)**2 + a(3,3)*a(1,2)**2
  35. & - a(1,1)*a(2,2)*a(3,3)
  36. call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3)
  37. d(1)=max(d1,d2,d3)
  38. d(3)=min(d1,d2,d3)
  39. d(2)=d1+d2+d3-(d(1)+d(3))
  40. c CLB dans le cas ou toutes les VP sont negatives,
  41. c il faut prendre la valeur absolue
  42. deps=abs (d(1))*1.e-6
  43. C CLB si l'ordre a ete inverse par le calcul (cf comm precedent)
  44. c on prend la valeur absolue
  45. if (abs(d(1)-d(2)).le.deps) then
  46. * valeur propre double
  47. if (abs(d(2)-d(3)).le.deps) then
  48. * valeur propre triple
  49. call vectp(a,d(1),x(1,1),3)
  50. else
  51. call vectp(a,d(1),x(1,1),2)
  52. call vectp(a,d(3),x(1,3),1)
  53. endif
  54. else
  55. deps=d(2)*1.e-6
  56. if (abs(d(2)-d(3)).le.deps) then
  57. * valeur propre double
  58. call vectp(a,d(1),x(1,1),1)
  59. call vectp(a,d(2),x(1,2),2)
  60. else
  61. * cas normal
  62. call vectp(a,d(1),x(1,1),1)
  63. call vectp(a,d(2),x(1,2),1)
  64. call vectp(a,d(3),x(1,3),1)
  65. endif
  66. endif
  67.  
  68. return
  69. end
  70.  
  71.  
  72.  
  73.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales