Télécharger la plaquette Cast3M


fr en

Télécharger jacod2.eso

Retour à la liste

Numérotation des lignes : 

jacod2

  1. C JACOD2    SOURCE    GOUNAND   26/01/14    21:15:06     12449          
  2.       SUBROUTINE JACOD2(A,D)
  3. C=====================================================================
  4. C
  5. C     DIAGONALISATION DE A      SYMETRIQUE
  6. C
  7. C     sans sortir les vecteurs propres
  8. C
  9. C     2025/10 Gounand : on reprend et simplifie jacob2.eso
  10. C
  11. C   ENTREEES
  12. C        A(3,3) = MATRICE A DIAGONALISER  A(1,3)=A(2,3)=0.
  13. C                                         A(3,1)=A(3,2)=0.
  14. C   SORTIES
  15. C        D(3)   =  LES 3 VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3)
  16. C
  17. C        RECUPERATION INCA FEVRIER 85    EBERSOLT
  18. C====================================================================
  19.       IMPLICIT INTEGER(I-N)
  20.       IMPLICIT REAL*8(A-H,O-Z)
  21. C
  22. -INC PPARAM
  23. -INC CCOPTIO
  24. -INC CCREEL
  25. C
  26.       DIMENSION A(3,*),D(*)
  27.       DIMENSION r(3,3)
  28. *
  29.       xpet=xpetit*10.D0
  30. *
  31.       D(3)=A(3,3)
  32. * Scaling de A
  33.       amax=0.d0
  34.       do i=1,2
  35.          do j=1,2
  36.             amax=max(amax,abs(a(j,i)))
  37.          enddo
  38.       enddo
  39. * Matrice nulle
  40.       if (amax.lt.xpet) then
  41.          D(1)=XZERO
  42.          D(2)=XZERO
  43.          return
  44.       endif
  45. C
  46.       amax1=1.D0/amax
  47.       do i=1,2
  48.          do j=1,2
  49.             R(j,i)=a(j,i)*amax1
  50.          enddo
  51.       enddo
  52. * Trace J1(A) et determinant J2(A) de A
  53.       xj1=R(1,1)+R(2,2)
  54.       xj2=R(1,1)*R(2,2)-R(1,2)**2
  55. *     Le discriminant reduit J4(A) du polynome caracteristique en l :
  56. *       x^2 -  J1(A) x + J2(A) = 0
  57. *     est : J4(A) = (J1(A) / 2) ^ 2 - J2(A)
  58. *
  59. *     Construisons S le deviateur de A : S = A - J1(A)/2 I
  60. *     On a J1(S)=0 et lambda_A = lambda_S + J1(A)/2
  61. *
  62. *     A et S ont les memes vecteurs propres
  63. *     Le polynome caracteristique de S est :
  64. *       x^2 + J2(S) = 0
  65. *     avec J2(S) = - [ ((a11-a22)/2)^2 + a12^2 ] = -J4(A)
  66. *     Cette expression de J4 est interessante car somme de carres donc positive
  67. *
  68.       xj4 = ((R(1,1)-R(2,2))/2)**2+R(1,2)**2
  69.       xsj4=sqrt(xj4)
  70.       stra=sign(1.d0,xj1)
  71. * Formules de Fagnano modifiees pour les deux racines
  72. * cf. The Ins and Outs of Solving Quadratic Equations with Floating-Point
  73. *     Arithmetic, Frederic Goualard, HAL preprint
  74.       x1=(xj1/2)+stra*xsj4
  75.       if (x1.eq.0.d0) then
  76.          x2=0.d0
  77.       else
  78.          x2=xj2/x1
  79.       endif
  80.       if (stra.gt.0) then
  81.          d(1)=x1
  82.          d(2)=x2
  83.       else
  84.          d(1)=x2
  85.          d(2)=x1
  86.       endif
  87. * Scaling des valeurs propres
  88.       d(1)=d(1)*amax
  89.       d(2)=d(2)*amax
  90.       RETURN
  91.       END
  92.  
  93.  

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