Télécharger la plaquette Cast3M


fr en

Télécharger jacoba.eso

Retour à la liste

Numérotation des lignes : 

jacoba

  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