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    GOUNAND   25/10/27    21:15:02     12387          
  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 Gounand 2025/10 On renvoie sur jacob3.eso plus robuste et precise
  24. C===============================================================
  25.       IMPLICIT INTEGER(I-N)
  26.       IMPLICIT REAL*8 (A-H,O-Z)
  27.       dimension a(3,3),d(3),x(3,3)
  28.       call jacob3(a,idim,d,x)
  29.       return
  30. * Old source
  31.       if (idim.ne.3) then
  32.         call jacob2(a,d,x)
  33.         return
  34.       endif
  35.       c3=-1.e0
  36.       c2=a(1,1)+a(2,2)+a(3,3)
  37.       c1=-(a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1))+
  38.      >  a(1,3)**2+a(1,2)**2+a(2,3)**2
  39.       c0=2*a(1,2)*a(1,3)*a(2,3)-a(1,1)*a(2,3)**2-
  40.      >  a(2,2)*a(1,3)**2-a(3,3)*a(1,2)**2+
  41.      >  a(1,1)*a(2,2)*a(3,3)
  42.       c0=c0/c3
  43.       c1=c1/c3
  44.       c2=c2/c3
  45.       call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3)
  46.       d(1)=max(d1,d2,d3)
  47.       d(3)=min(d1,d2,d3)
  48.       d(2)=d1+d2+d3-d(1)-d(3)
  49.       deps=max(abs(d1),abs(d2),abs(d3))*1e-4
  50. * A ttention j'ai du mettre des val abs
  51. * car ça ne marchait pas trop bien lorsque les
  52. * d(1) tend vers 0+
  53. * d(2) tend vers 0-
  54. * d(3) négatif
  55. * comme on avait deps tends vers 0+
  56. * d(1) - d(2) = 2.*0+
  57. * il considérait ce cas comme normal
  58. * alors qu'on a une valeur propre double
  59. *        PRINT *,'d1,d2,d3',d1,d2,d3
  60. *        PRINT *,'d(1),d(2),d(3),deps',d(1),d(2),d(3),deps
  61. *        PRINT *,'c0,c1,c2',c0,c1,c2
  62.       if ((abs(d(1)-d(2))).le.deps) then
  63. *      valeur propre double
  64.        if ((abs(d(2)-d(3))).le.deps) then
  65. *       valeur propre triple
  66.         call vectp(a,d(1),x(1,1),3)
  67.        else
  68.         call vectp(a,d(1),x(1,1),2)
  69.         call vectp(a,d(3),x(1,3),1)
  70.        endif
  71.       else
  72.        if ((abs(d(2)-d(3))).le.deps) then
  73.  
  74. *       valeur propre double
  75. *        PRINT *,'JACOBA VP BOUBLE'
  76.  
  77.         call vectp(a,d(1),x(1,1),1)
  78.         call vectp(a,d(2),x(1,2),2)
  79.        else
  80. *       cas normal
  81. *        PRINT *,'JACOBA CAS normal'
  82.         call vectp(a,d(1),x(1,1),1)
  83.         call vectp(a,d(2),x(1,2),1)
  84.         call vectp(a,d(3),x(1,3),1)
  85.        endif
  86.       endif
  87.       end
  88. *
  89.  
  90.  

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