Télécharger jacoba.eso

Retour à la liste

Numérotation des lignes :

  1. C JACOBA SOURCE CMAU 06/01/10 21:15:29 5284
  2. C JACOB3 SOURCE CHAT 05/01/13 00:48:17 5004
  3. subroutine jacoba(a,idim,d,x)
  4. C======================================================================
  5. C routine calquée sur jacob3.eso avec ajout de valeur absolue au moment
  6. C des test
  7. C OBJET
  8. C -----
  9. C DIAGONALISATION D UNE MATRICE 3*3 SYMETRIQUE
  10. C
  11. C ENTREES
  12. C -------
  13. C A(3,3) = MATRICE SYMETRIQUE
  14. C IDIM = 2 OU 3 SI 2 ON NE S OCCUPE QUE DE A(2,2)
  15. C SI 3 DE A(3,3)
  16. C SORTIES
  17. C -------
  18. C D(3) = VALEURS PROPRES ORDONNEES D(1)>D(2)>D(3)
  19. C
  20. C X(3,3) = VECTEURS PROPRES ( X(IP,2) EST LE VECTEUR
  21. C ASSOCIE A D(2) )
  22. C
  23. C===============================================================
  24. IMPLICIT INTEGER(I-N)
  25. IMPLICIT REAL*8 (A-H,O-Z)
  26. dimension a(3,3),d(3),x(3,3)
  27. if (idim.ne.3) then
  28. call jacob2(a,d,x)
  29. return
  30. endif
  31. c3=-1.e0
  32. c2=a(1,1)+a(2,2)+a(3,3)
  33. c1=-(a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1))+
  34. > a(1,3)**2+a(1,2)**2+a(2,3)**2
  35. c0=2*a(1,2)*a(1,3)*a(2,3)-a(1,1)*a(2,3)**2-
  36. > a(2,2)*a(1,3)**2-a(3,3)*a(1,2)**2+
  37. > a(1,1)*a(2,2)*a(3,3)
  38. c0=c0/c3
  39. c1=c1/c3
  40. c2=c2/c3
  41. call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3)
  42. d(1)=max(d1,d2,d3)
  43. d(3)=min(d1,d2,d3)
  44. d(2)=d1+d2+d3-d(1)-d(3)
  45. deps=max(abs(d1),abs(d2),abs(d3))*1e-4
  46. * A ttention j'ai du mettre des val abs
  47. * car ça ne marchait pas trop bien lorsque les
  48. * d(1) tend vers 0+
  49. * d(2) tend vers 0-
  50. * d(3) négatif
  51. * comme on avait deps tends vers 0+
  52. * d(1) - d(2) = 2.*0+
  53. * il considérait ce cas comme normal
  54. * alors qu'on a une valeur propre double
  55. * PRINT *,'d1,d2,d3',d1,d2,d3
  56. * PRINT *,'d(1),d(2),d(3),deps',d(1),d(2),d(3),deps
  57. * PRINT *,'c0,c1,c2',c0,c1,c2
  58. if ((abs(d(1)-d(2))).le.deps) then
  59. * valeur propre double
  60. if ((abs(d(2)-d(3))).le.deps) then
  61. * valeur propre triple
  62. call vectp(a,d(1),x(1,1),3)
  63. else
  64. call vectp(a,d(1),x(1,1),2)
  65. call vectp(a,d(3),x(1,3),1)
  66. endif
  67. else
  68. if ((abs(d(2)-d(3))).le.deps) then
  69.  
  70. * valeur propre double
  71. * PRINT *,'JACOBA VP BOUBLE'
  72.  
  73. call vectp(a,d(1),x(1,1),1)
  74. call vectp(a,d(2),x(1,2),2)
  75. else
  76. * cas normal
  77. * PRINT *,'JACOBA CAS normal'
  78. call vectp(a,d(1),x(1,1),1)
  79. call vectp(a,d(2),x(1,2),1)
  80. call vectp(a,d(3),x(1,3),1)
  81. endif
  82. endif
  83. end
  84. *
  85.  
  86.  
  87.  
  88.  
  89.  

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