Télécharger la plaquette Cast3M


fr en

Télécharger jacob3.eso

Retour à la liste

Numérotation des lignes : 

jacob3

  1. C JACOB3    SOURCE    PV090527  23/02/15    21:15:04     11592          
  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.       if (d1.lt.d2) then
  38.            dd=d1
  39.            d1=d2
  40.            d2=dd
  41.       endif
  42.       if (d1.lt.d3) then
  43.            dd=d1
  44.            d1=d3
  45.            d3=dd
  46.       endif
  47.       if (d2.lt.d3) then
  48.            dd=d2
  49.            d2=d3
  50.            d3=dd
  51.       endif
  52.       d(1)=d1                      
  53.       d(2)=d2                      
  54.       d(3)=d3                      
  55. c CLB dans le cas ou toutes les VP sont negatives,
  56. c     il faut prendre la valeur absolue
  57.       deps=abs (d(1))*1.e-6
  58. C CLB si l'ordre a ete inverse par le calcul (cf comm precedent) 
  59. c     on prend la valeur absolue
  60.       if (abs(d(1)-d(2)).le.deps) then
  61. *      valeur propre double
  62.        if (abs(d(2)-d(3)).le.deps) then
  63. *       valeur propre triple
  64.         call vectp(a,d(1),x(1,1),3)
  65.        else
  66.         call vectp(a,d(1),x(1,1),2)
  67.         call vectp(a,d(3),x(1,3),1)
  68.        endif
  69.       else
  70.        deps=d(2)*1.e-6
  71.        if (abs(d(2)-d(3)).le.deps) then
  72. *       valeur propre double
  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.         call vectp(a,d(1),x(1,1),1)
  78.         call vectp(a,d(2),x(1,2),1)
  79.         call vectp(a,d(3),x(1,3),1)
  80.        endif
  81.       endif
  82.  
  83.       return
  84.       end
  85.  
  86.  
  87.  
  88.  
  89.  

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