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    GOUNAND   25/10/23    21:15:05     12385          
  2.       subroutine jacob3(a,idima,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     IDIMA   = 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     X(3,3) = VECTEURS PROPRES ( X(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. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC CCREEL
  26.       dimension a(3,3),d(3),x(3,3),s(3,3),as(3,3)
  27.       logical ldbg,lverif
  28.       dimension r(3,3),v(3)
  29. *
  30. * Verifications
  31.       lverif=.false.
  32. * Impressions
  33.       ldbg=.false.
  34.       xpre=xzprec*10.D0
  35.       xpet=xpetit*10.D0
  36.       xprev=xzprec*10.D0
  37. *      xpre=xzprec*30.D0
  38. *      xpet=xpetit*30.D0
  39.       if (idima.ne.3) then
  40.         call jacob2(a,d,x)
  41.         return
  42.       endif
  43. *
  44. *     SCALING DE A
  45. *
  46.       amax=0.d0
  47.       do i=1,idima
  48.          do j=1,idima
  49.             amax=max(amax,abs(a(j,i)))
  50.          enddo
  51.       enddo
  52. * Matrice nulle
  53.       if (amax.lt.xpet) then
  54.          do i=1,idima
  55.             d(i)=0.D0
  56.             do j=1,idima
  57.                if (i.eq.j) then
  58.                   x(j,i)=1.d0
  59.                else
  60.                   x(j,i)=0.d0
  61.                endif
  62.             enddo
  63.          enddo
  64.          return
  65.       endif
  66. C
  67.       amax1=1.D0/amax
  68.       do i=1,idima
  69.          do j=1,idima
  70.             as(j,i)=a(j,i)*amax1
  71.          enddo
  72.       enddo
  73. *
  74. *     CALCUL DES INVARIANTS
  75. *
  76. * On applique la methode de Harari et Albocher IJNME'2023
  77. * pour calculer les invariants et les valeurs propres
  78. * Invariants J1=TRACE, J2, J3=DET
  79.       xj1=as(1,1)+as(2,2)+as(3,3)
  80.       do i=1,idima
  81.          do j=1,idima
  82.             if (i.eq.j) then
  83.                s(j,i)=as(j,i)-xj1/3
  84.             else
  85.                s(j,i)=as(j,i)
  86.             endif
  87.          enddo
  88.       enddo
  89.       if (ldbg) then
  90.          write(6,*) 'Matrice S='
  91.          do ii=1,idima
  92.             write (6,*) (S(ii,jj),jj=1,idima)
  93.          enddo
  94.       endif
  95. *     Dans la methode, on veut les invariants J2(S) et J3(S)
  96. *     où S est le deviateur de A (A = S + J1(A)/3 I)
  97. *     On a J1(S)=0 et lambda_A = lambda_S + J1(A)/3
  98. *
  99. *     A et S ont les memes vecteurs propres
  100. *
  101. *     On exprime les invariants de S de maniere stable en fonction des
  102. *     termes de A en particulier les differences diagonales
  103. *     En particulier J2 et J4 sont exprimes sous forme de sommes de
  104. *     carres donc positif. On peut donc en prendre les racines carrees
  105. *     sans fremir
  106.       d12=as(1,1)-as(2,2)
  107.       d23=as(2,2)-as(3,3)
  108.       d31=as(3,3)-as(1,1)
  109. * Expression de J2(S)
  110.       xj2d=(d12**2+d23**2+d31**2)/6
  111.       xj2o=as(1,2)**2+as(2,3)**2+as(3,1)**2
  112.       xj2=xj2d+xj2o
  113. * Expression du determinant J3(S)
  114.       xj3d=(d12-d31)*(d23-d12)*(d31-d23)/27
  115.       xj3m=(d12-d31)*as(2,3)**2+(d23-d12)*as(3,1)**2
  116.      $     +(d31-d23)*as(1,2)**2
  117.       xj3m=-xj3m/3
  118.       xj3o=2.d0*as(1,2)*as(2,3)*as(3,1)
  119.       xj3=xj3d+xj3m+xj3o
  120. * Expression du discriminant Delta(S)=J4(S)=4 J2(S)^3 - 27 J3(S)^2
  121.       xj4x=d12*d23*d31
  122.      $     + as(1,2)**2*d12 + as(2,3)**2*d23 + as(1,3)**2*d31
  123.       xj4y1=as(2,3)*(2*as(2,3)**2 - as(1,3)**2 - as(1,2)**2 + 2*d12*d31)
  124.      $     +as(1,2)*as(1,3)*(d12-d31)
  125.       xj4y2=as(1,3)*(2*as(1,3)**2 - as(2,3)**2 - as(1,2)**2 + 2*d23*d12)
  126.      $     +as(1,2)*as(2,3)*(d23-d12)
  127.       xj4y3=as(1,2)*(2*as(1,2)**2 - as(2,3)**2 - as(1,3)**2 + 2*d31*d23)
  128.      $     +as(1,3)*as(2,3)*(d31-d23)
  129.       xj4z1=as(2,3)*(as(1,3)**2 - as(1,2)**2) + as(1,2)*as(1,3)*d23
  130.       xj4z2=as(1,3)*(as(1,2)**2 - as(2,3)**2) + as(2,3)*as(1,2)*d31
  131.       xj4z3=as(1,2)*(as(2,3)**2 - as(1,3)**2) + as(1,3)*as(2,3)*d12
  132. *
  133.       xj4=xj4x**2+xj4y1**2+xj4y2**2+xj4y3**2
  134.      $     + 15* (xj4z1**2+xj4z2**2+xj4z3**2)
  135.       if (ldbg) then
  136. * Impression des invariants de S
  137.          c2=-s(1,1)-s(2,2)-s(3,3)
  138.          c1= (s(1,1)*s(2,2)+s(2,2)*s(3,3)+s(3,3)*s(1,1))
  139.      &        - s(1,3)**2 - s(1,2)**2 - s(2,3)**2
  140.          c0=-2.d0*s(1,2)*s(1,3)*s(2,3) + s(1,1)*s(2,3)**2
  141.      &        + s(2,2)*s(1,3)**2 + s(3,3)*s(1,2)**2
  142.      &        - s(1,1)*s(2,2)*s(3,3)
  143.          xj1b=-c2+xj1
  144.          xj2b=-c1
  145.          xj3b=-c0
  146.          xj4b=4*xj2b**3-27*xj3b**2
  147.          write(6,*) 'Invariants du deviateur S : new old  abs(old-new)'
  148.          write(6,*) ' trace= ',xj1,xj1b,abs(xj1b-xj1)
  149.          write(6,*) ' J2=    ',xj2,xj2b,abs(xj2b-xj2)
  150.          write(6,*) ' det=   ',xj3,xj3b,abs(xj3b-xj3)
  151.          write(6,*) ' discr= ',xj4,xj4b,abs(xj4b-xj4)
  152.       endif
  153. *
  154. *     CALCUL DES VALEURS PROPRES
  155. *
  156.       sdet=sign(1.d0,xj3)
  157.       xsj4=sqrt(xj4)
  158.       xsj2=sqrt(xj2)
  159.       xdenom=2*xj2*xsj2+3*sqrt(3.d0)*sdet*xj3
  160.       ang=2.d0/3.d0*atan2(xsj4,xdenom)
  161.       xf1=xsj2*sdet
  162.       xf2=cos(ang)/sqrt(3.d0)
  163.       xf3=sin(ang)
  164. * Si sdet=1   ds1 .GE. ds2 .GE. ds3
  165. * Si sdet=-1  ds3 .GE. ds2 .GE. ds1
  166. * dsi sont les valeurs propres de S
  167. * dAi sont les valeurs propres de A
  168.       ds1=2*xf1*xf2
  169.       ds2=xf1*(xf3-xf2)
  170.       ds3=xf1*(-xf3-xf2)
  171.       xj13=xj1/3
  172.       das1=ds1+xj13
  173.       das2=ds2+xj13
  174.       das3=ds3+xj13
  175. * Impression des valeurs propres de A par les methodes nouvelles et anciennes
  176.       if (ldbg) then
  177.          d(2)=das2
  178.          if (sdet.gt.0) then
  179.             d(1)=das1
  180.             d(3)=das3
  181.          else
  182.             d(3)=das1
  183.             d(1)=das3
  184.          endif
  185. *   Methode ancienne
  186.          c2=-as(1,1)-as(2,2)-as(3,3)
  187.          c1= (as(1,1)*as(2,2)+as(2,2)*as(3,3)+as(3,3)*as(1,1))
  188.      &        - as(1,3)**2 - as(1,2)**2 - as(2,3)**2
  189.          c0=-2.d0*as(1,2)*as(1,3)*as(2,3) + as(1,1)*as(2,3)**2
  190.      &        + as(2,2)*as(1,3)**2 + as(3,3)*as(1,2)**2
  191.      &        - as(1,1)*as(2,2)*as(3,3)
  192.          write(6,*) 'c0,c1,c2=',c0,c1,c2
  193.          call degre3(c0,c1,c2,d1,XI1,d2,XI2,d3,XI3)
  194.          write(6,*) 'Som vp=',d1+d2+d3
  195.          write(6,*) 'Prod vp=',d1*d2*d3
  196.          if (d1.lt.d2) then
  197.             dd=d1
  198.             d1=d2
  199.             d2=dd
  200.          endif
  201.          if (d1.lt.d3) then
  202.             dd=d1
  203.             d1=d3
  204.             d3=dd
  205.          endif
  206.          if (d2.lt.d3) then
  207.             dd=d2
  208.             d2=d3
  209.             d3=dd
  210.          endif
  211.          db1=d1
  212.          db2=d2
  213.          db3=d3
  214.          write(6,*) 'Val. propres de A scalee : new old  abs(old-new)'
  215.          write(6,*) ' vp1= ',d(1),db1,abs(db1-d(1))
  216.          write(6,*) ' vp2= ',d(2),db2,abs(db2-d(2))
  217.          write(6,*) ' vp3= ',d(3),db3,abs(db3-d(3))
  218.       endif
  219. *
  220. *     Faire le calcul des vecteurs propres avec S plutôt qu'avec A
  221. *     (moins d'erreur de cancellation)
  222. *     On utilise la méthode de deflation de Scherzinger et Dohrmann
  223. *     CMAME'08 car sinon, on n'a pas reussi a avoir une bonne
  224. *     orthogonalite des vecteurs propres lorsque les valeurs propres
  225. *     sont proches
  226. *
  227.       d(2)=ds2
  228. * Si sdet est positif, d(1) est la valeur propre la plus isolee
  229. * Si sdet est negatif, d(3) est la valeur propre la plus isolee
  230. * On recherche d'abord le vecteur propre associe a cette valeur
  231. *     propre
  232.       if (sdet.gt.0) then
  233.          igr=1
  234.       else
  235.          igr=3
  236.       endif
  237. * igr vaut 1 ou 3 quand ipe vaut 3 ou 1
  238.       ipe=4-igr
  239.       d(igr)=ds1
  240.       d(ipe)=ds3
  241.       if (ldbg) then
  242.          write(6,*) 'Val. propres de S'
  243.          write(6,*) ' vp1= ',d(1)
  244.          write(6,*) ' vp2= ',d(2)
  245.          write(6,*) ' vp3= ',d(3)
  246.       endif
  247. *
  248. *     CALCUL DES VECTEURS PROPRES
  249. *
  250. *     Raccourci lorsque les trois valeurs propres sont egales
  251. *     Dans ce cas S=0 et J2(S)=0
  252. *
  253.       xsca=max(xpre,xpet)
  254.       if (ldbg) then
  255.          write(6,*) 'xpre,xpet,xsca,xsj2=',xpre,xpet,xsca,xsj2
  256.       endif
  257. *
  258.       if (xsj2.lt.xsca) then
  259.          if (ldbg) write(6,*) '3 valeurs propres egales detectees'
  260.          do i=1,idima
  261.             do j=1,idima
  262.                if (i.eq.j) then
  263.                   x(j,i)=1.d0
  264.                else
  265.                   x(j,i)=0.d0
  266.                endif
  267.             enddo
  268.          enddo
  269.          goto 1000
  270.       endif
  271. *
  272. *     Faire le calcul des vecteurs propres avec S plutôt qu'avec A
  273. *     (moins d'erreur de cancellation)
  274. *     On utilise la méthode de deflation de Scherzinger et Dohrmann
  275. *     CMAME'08 car sinon, on n'a pas reussi a avoir une bonne
  276. *     orthogonalite des vecteurs propres lorsque les valeurs propres
  277. *     sont proches
  278. *
  279. *     RECHERCHE DU 1er VECTEUR PROPRE associe a la valeur propre isolee
  280. *     d(igr)
  281. *
  282.       call vectp2(s,d,igr,x)
  283.       if (lverif) then
  284. *     Verif
  285.          call vmorth(x,idima,3,xpre,r)
  286.          if (ierr.ne.0) then
  287.             write(6,*) '!!! X non orthonormale'
  288.             if (ldbg) then
  289.                write(6,*) 'Matrice des VPs X='
  290.                do ii=1,idima
  291.                   write (6,*) (X(ii,jj),jj=1,idima)
  292.                enddo
  293.                write(6,*) 'Matrice des Pscals R='
  294.                do ii=1,idima
  295.                   write (6,*) (R(ii,jj),jj=1,idima)
  296.                enddo
  297.             endif
  298.             goto 9999
  299.          endif
  300. *     vp(igr) est-il bien le vp associe a la valeur propre d(igr)
  301.          call vvectp(s,idima,3,d(igr),x(1,igr),xprev,v)
  302.          if (ierr.ne.0) then
  303.             write(6,*) '!!!   vp(igr) nest pas vecteur propre '
  304.             if (ldbg) then
  305.                write(6,*) ' v= (S - D(igr)) x(j,igr)  ='
  306.                write(6,*) (v(jj),jj=1,idima)
  307.             endif
  308.             goto 9999
  309.          endif
  310.       endif
  311. *
  312. *     FIN DE RECHERCHE DU 1er VECTEUR PROPRE
  313. *
  314. *
  315. *     Raccourci lorsque les deux autres valeurs propres sont egales
  316. *     Dans ce cas le discriminant est nul J4(S)=0
  317. *     Le developpement limite de lambda_2 et lambda_3 en J4 tendant vers
  318. *     0 est :
  319. *     lambda_2,3 \propto sqrt(J2) ( 1 + sqrt(J4) / (J2^3/2)  )
  320. *
  321.  
  322.       xsca=max((xsj2**3)*xpre,xpet)
  323.       if (ldbg) then
  324.          write(6,*) 'xsj2,xsj23,xpre,xsca,xsj4=',xsj2,xsj2**3,xpre,xsca
  325.      $        ,xsj4
  326.       endif
  327.       if (xsj4.lt.xsca) then
  328.          if (ldbg) write(6,*) '2 valeurs propres egales detectees'
  329. *     La matrice X contient une base donc on peut quitter prematurement
  330.          goto 1000
  331.       endif
  332. *
  333. *    RECHERCHE des 2emes et 3emes VECTEURS PROPRES
  334. *    les valeurs propres sont normalement numeriquement distinctes
  335. *
  336. *     Recherche un vecteur propre x(j,2) associe à la valeur propre
  337. *     simple d(2)  d'une matrice 3x3 S
  338. *
  339. *     L'image de S-d(2)I est de dimension 2 car d(2) est valeur
  340. *     propre simple. On cherche uniquement la partie de l'image
  341. *     orthogonale à vp(igr).
  342. *     Pour ce faire, on prend la plus grande des images des deux
  343. *     vecteurs de la base X différents de vp(igr). Ces deux images sont
  344. *     colinéaires car à la fois orthogonales à vp(igr) et vp(2) : w1
  345. *
  346. *     On choisit ensuite :  vp(2) = w1 \pvec vp(1)
  347. *                  puis  :  vp(3) = vp(1) \pvec vp(2)
  348. *
  349. *
  350.       call vectp3(s,d,igr,x)
  351. *
  352. * Verif
  353.       if (lverif) then
  354. * Orthonormalite des vps
  355.          call vmorth(x,idima,3,xpre,r)
  356.          if (ierr.ne.0) then
  357.             write(6,*) '!!! X non orthonormale'
  358.             if (ldbg) then
  359.                write(6,*) 'Matrice des VPs X='
  360.                do ii=1,idima
  361.                   write (6,*) (X(ii,jj),jj=1,idima)
  362.                enddo
  363.                write(6,*) 'Matrice des Pscals R='
  364.                do ii=1,idima
  365.                   write (6,*) (R(ii,jj),jj=1,idima)
  366.                enddo
  367.             endif
  368.             goto 9999
  369.          endif
  370.          do i=1,idima
  371. *     vp(i) est-il bien le vp associe a la valeur propre d(i)
  372.             call vvectp(s,idima,3,d(i),x(1,i),xprev,v)
  373.             if (ierr.ne.0) then
  374.                write(6,*) '!!!   vp(i) nest pas vecteur propre i=',i
  375.                if (ldbg) then
  376.                   write(6,*) ' v= (S - D(i)) x(j,i)  ='
  377.                   write(6,*) (v(jj),jj=1,idima)
  378.                endif
  379.                goto 9999
  380.             endif
  381.          enddo
  382.       endif
  383. *
  384.       if (.false.) then
  385. * Methode ancienne
  386.          dmaxa=max(abs(d(1)),abs(d(2)),abs(d(3)))
  387.          deps=max(dmaxa*xpre,xpetit*10.d0)
  388.          if (abs(d(1)-d(2)).le.deps) then
  389. *     valeur propre double
  390.             if (abs(d(2)-d(3)).le.deps) then
  391. *     valeur propre triple
  392.                if (ldbg) write(6,*) 'JACOB3 Valeur propre triple'
  393.                call vectp(s,d(1),x(1,1),3)
  394.             else
  395.                if (ldbg) write(6,*) 'JACOB3 Valeur propre double 1'
  396.                call vectp(a,d(1),x(1,1),2)
  397.                x(1,3)=x(2,1)*x(3,2)-x(3,1)*x(2,2)
  398.                x(2,3)=x(3,1)*x(1,2)-x(1,1)*x(3,2)
  399.                x(3,3)=x(1,1)*x(2,2)-x(2,1)*x(1,2)
  400. *old     call vectp(a,d(3),x(1,3),1)
  401.             endif
  402.          else
  403.             if (abs(d(2)-d(3)).le.deps) then
  404.                if (ldbg) write(6,*) 'JACOB3 Valeur propre double 2'
  405. *       valeur propre double
  406.                call vectp(a,d(2),x(1,2),2)
  407.                x(1,1)=x(2,2)*x(3,3)-x(3,2)*x(2,3)
  408.                x(2,1)=x(3,2)*x(1,3)-x(1,2)*x(3,3)
  409.                x(3,1)=x(1,2)*x(2,3)-x(2,2)*x(1,3)
  410. *old            call vectp(a,d(1),x(1,1),1)
  411.             else
  412. *       cas normal
  413.                if (ldbg) write(6,*) 'JACOB3 Cas normal'
  414.                call vectp(s,d(1),x(1,1),1)
  415.                call vectp(s,d(2),x(1,2),1)
  416.                call vectp(s,d(3),x(1,3),1)
  417.             endif
  418.          endif
  419.       endif
  420. *
  421. * On remet dans D les valeurs propres de A a la place de celles de S
  422. *
  423.  1000 continue
  424.       d(2)=das2*amax
  425.       if (sdet.gt.0) then
  426.          d(1)=das1*amax
  427.          d(3)=das3*amax
  428.       else
  429.          d(3)=das1*amax
  430.          d(1)=das3*amax
  431.       endif
  432.       return
  433. *
  434. * Error handling
  435. *
  436.  
  437.  9999 CONTINUE
  438.       MOTERR(1:8)='JACOB3'
  439.       call erreur(1127)
  440.       return
  441.       end
  442.  
  443.  

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