Esope source files

   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
  27 :         call jacob2(a,d,x)
  28 :         return
  29 :       endif
  30 :       c3=-1.e0
  31 :       c2=a(1,1)+a(2,2)+a(3,3)
  32 :       c1=-(a(1,1)*a(2,2)+a(2,2)*a(3,3)+a(3,3)*a(1,1))+
  33 :      >  a(1,3)**2+a(1,2)**2+a(2,3)**2
  34 :       c0=2*a(1,2)*a(1,3)*a(2,3)-a(1,1)*a(2,3)**2-
  35 :      >  a(2,2)*a(1,3)**2-a(3,3)*a(1,2)**2+
  36 :      >  a(1,1)*a(2,2)*a(3,3)
  37 :       c0=c0/c3
  38 :       c1=c1/c3
  39 :       c2=c2/c3
  39 :       call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3)
  40 :       d(1)=max(d1,d2,d3)
  41 :       d(3)=min(d1,d2,d3)
  42 :       d(2)=d1+d2+d3-d(1)-d(3)
  43 :       deps=max(abs(d1),abs(d2),abs(d3))*1e-4
  44 : * A ttention j'ai du mettre des val abs 
  45 : * car ça ne marchait pas trop bien lorsque les
  46 : * d(1) tend vers 0+ 
  47 : * d(2) tend vers 0- 
  48 : * d(3) négatif 
  49 : * comme on avait deps tends vers 0+ 
  50 : * d(1) - d(2) = 2.*0+
  51 : * il considérait ce cas comme normal 
  52 : * alors qu'on a une valeur propre double  
  53 : *        PRINT *,'d1,d2,d3',d1,d2,d3
  54 : *        PRINT *,'d(1),d(2),d(3),deps',d(1),d(2),d(3),deps
  55 : *        PRINT *,'c0,c1,c2',c0,c1,c2
  56 :       if ((abs(d(1)-d(2))).le.deps) then
  57 : *      valeur propre double
  58 :        if ((abs(d(2)-d(3))).le.deps) then
  59 : *       valeur propre triple
  59 :         call vectp(a,d(1),x(1,1),3)
  60 :        else
  60 :         call vectp(a,d(1),x(1,1),2)
  60 :         call vectp(a,d(3),x(1,3),1)
  61 :        endif
  62 :       else
  63 :        if ((abs(d(2)-d(3))).le.deps) then
  64 : 
  65 : *       valeur propre double 
  66 : *        PRINT *,'JACOBA VP BOUBLE' 
  67 : 
  67 :         call vectp(a,d(1),x(1,1),1)
  67 :         call vectp(a,d(2),x(1,2),2)
  68 :        else
  69 : *       cas normal
  70 : *        PRINT *,'JACOBA CAS normal'  
  70 :         call vectp(a,d(1),x(1,1),1)
  70 :         call vectp(a,d(2),x(1,2),1)
  70 :         call vectp(a,d(3),x(1,3),1)
  71 :        endif
  72 :       endif
  73 :       end
  74 : *
  75 :  
  76 :  
  77 :  
  78 :