Télécharger la plaquette Cast3M


fr en

Télécharger jacob2.eso

Retour à la liste

Numérotation des lignes : 

jacob2

  1. C JACOB2    SOURCE    GOUNAND   26/01/16    21:15:05     12450          
  2.       SUBROUTINE JACOB2(A,D,S)
  3. C=====================================================================
  4. C
  5. C       DIAGONALISATION DE A      SYMETRIQUE
  6. C   ENTREEES
  7. C        A(3,3) = MATRICE A DIAGONALISER  A(1,3)=A(2,3)=0.
  8. C                                         A(3,1)=A(3,2)=0.
  9. C   SORTIES
  10. C        D(3)   =  LES 3 VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3)
  11. C        S(3,3) =  VECTEURS PROPRES    S(I,3)=0. 0. 1.
  12. C
  13. C        RECUPERATION INCA FEVRIER 85    EBERSOLT
  14. C====================================================================
  15.       IMPLICIT INTEGER(I-N)
  16.       IMPLICIT REAL*8(A-H,O-Z)
  17. C
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC CCREEL
  21. C
  22.       DIMENSION A(3,*),D(*),S(3,*)
  23.       LOGICAL lverif,ldbg
  24.       DIMENSION r(3,3),xnr(3)
  25. C
  26. C Valeur particuliere : SQRT(2.)*0.5 = 1/SQRT(2.)
  27.       PARAMETER ( X1sr2 = 0.70710678118654752440084436210484904D0 )
  28. * Verifications
  29.       lverif=.false.
  30. * Impressions
  31.       ldbg=.false.
  32.       xpre=xzprec*3.D0
  33.       xpet=xpetit*10.D0
  34.       xprev=xzprec*10.D0
  35. *
  36.       D(3)=A(3,3)
  37.       S(3,1)=0.D0
  38.       S(3,2)=0.D0
  39.       S(1,3)=0.D0
  40.       S(2,3)=0.D0
  41.       S(3,3)=1.D0
  42. * Scaling de A
  43.       amax=0.d0
  44.       do i=1,2
  45.          do j=1,2
  46.             amax=max(amax,abs(a(j,i)))
  47.          enddo
  48.       enddo
  49. * Matrice nulle
  50.       if (amax.lt.xpet) then
  51.          D(1)=XZERO
  52.          D(2)=XZERO
  53.          S(1,1)=1.D0
  54.          S(2,1)=XZERO
  55.          S(1,2)=XZERO
  56.          S(2,2)=1.D0
  57.          return
  58.       endif
  59. C
  60.       amax1=1.D0/amax
  61.       do i=1,2
  62.          do j=1,2
  63.             R(j,i)=a(j,i)*amax1
  64.          enddo
  65.       enddo
  66. * Trace J1(A) et determinant J2(A) de A
  67.       xj1=R(1,1)+R(2,2)
  68.       xj2=R(1,1)*R(2,2)-R(1,2)**2
  69. *     Le discriminant reduit J4(A) du polynome caracteristique en l :
  70. *       x^2 -  J1(A) x + J2(A) = 0
  71. *     est : J4(A) = (J1(A) / 2) ^ 2 - J2(A)
  72. *
  73. *     Construisons S le deviateur de A : S = A - J1(A)/2 I
  74. *     On a J1(S)=0 et lambda_A = lambda_S + J1(A)/2
  75. *
  76. *     A et S ont les memes vecteurs propres
  77. *     Le polynome caracteristique de S est :
  78. *       x^2 + J2(S) = 0
  79. *     avec J2(S) = - [ ((a11-a22)/2)^2 + a12^2 ] = -J4(A)
  80. *     Cette expression de J4 est interessante car somme de carres donc positive
  81. *
  82.       xj4 = ((R(1,1)-R(2,2))/2)**2+R(1,2)**2
  83.       xsj4=sqrt(xj4)
  84.       stra=sign(1.d0,xj1)
  85. *  Formules de Fagnano modifiees pour les deux racines
  86. * cf. The Ins and Outs of Solving Quadratic Equations with Floating-Point
  87. *     Arithmetic, Frederic Goualard, HAL preprint
  88.       x1=(xj1/2)+stra*xsj4
  89.       if (x1.eq.0.d0) then
  90.          x2=0.d0
  91.       else
  92.          x2=xj2/x1
  93.       endif
  94.       if (stra.gt.0) then
  95.          d(1)=x1
  96.          d(2)=x2
  97.       else
  98.          d(1)=x2
  99.          d(2)=x1
  100.       endif
  101.       x3=xsj4
  102. *
  103. * Ancienne methode
  104. *
  105.       if (.false.) then
  106.          X1  =2.D0*R(1,2)
  107.          X2  =R(1,1)-R(2,2)
  108.          X3  =SQRT(X2*X2+X1*X1)
  109.          D(1)=0.5D0*(R(1,1)+R(2,2)+X3)
  110.          D(2)=D(1)-X3
  111.       endif
  112. *     Raccourci lorsque les deux valeurs propres sont egales
  113. *
  114.       xsca=max(xpre,xpet)
  115.       if (ldbg) then
  116.          write(6,*) 'xpre,xpet,xsca,x3=',xpre,xpet,xsca,x3
  117.       endif
  118. *
  119.       if (x3.lt.xsca) then
  120.          if (ldbg) write(6,*) '2 valeurs propres egales detectees'
  121.          do i=1,2
  122.             do j=1,2
  123.                if (i.eq.j) then
  124.                   s(j,i)=1.d0
  125.                else
  126.                   s(j,i)=0.d0
  127.                endif
  128.             enddo
  129.          enddo
  130.          goto 1000
  131.       endif
  132. C
  133. C Appliquer la methode de Scherzinger et Dohrmann CMAME'08
  134. * Recherche du plus grand vecteur de l'image de A - d(1) I :
  135. C
  136.       do i=1,2
  137.          R(i,i)=R(i,i)-d(1)
  138.       enddo
  139.       imax=0
  140.       xmax=0.d0
  141.       do i=1,2
  142.          xnr(i)=0.d0
  143.          do j=1,2
  144.             xnr(i)=xnr(i)+r(j,i)**2
  145.          enddo
  146.          if (xnr(i).gt.xmax) then
  147.             xmax=xnr(i)
  148.             imax=i
  149.          endif
  150.       enddo
  151.       igr=imax
  152. * Normalement pas possible si valeurs propres non egales
  153.       if (igr.eq.0) then
  154.          igr=1
  155.          r(1,1)=1.d0
  156.          r(2,1)=0.d0
  157.       else
  158.          xnorm=sqrt(xnr(igr))
  159.          do i=1,2
  160.             r(i,igr)=r(i,igr)/xnorm
  161.          enddo
  162.       endif
  163.       S(1,2)= r(1,igr)
  164.       S(2,2)= r(2,igr)
  165.       S(1,1)= S(2,2)
  166.       S(2,1)=-S(1,2)
  167.  
  168.       if (.false.) then
  169. *       Methode ancienne
  170.          R11=A(1,1)-D(1)
  171.          R21=A(1,2)
  172.          XNR1=SQRT(R11*R11+R21*R21)
  173.          R12=A(2,1)
  174.          R22=A(2,2)-D(1)
  175.          XNR2=SQRT(R12*R12+R22*R22)
  176.          IF (XNR1.GE.XNR2) THEN
  177.             I=1
  178.             T1=R11
  179.             T2=R21
  180.             XT=XNR1
  181.          ELSE
  182.             I=2
  183.             T1=R12
  184.             T2=R22
  185.             XT=XNR2
  186.          ENDIF
  187. C
  188.          XSCA=MAX(MAX(ABS(D(1)),ABS(D(2)))*XZPREC,XPETIT)*2
  189.          IF (XT.LE.XSCA) THEN
  190. *         write(6,*) 'Cas 0'
  191.             X5    = X1sr2
  192.             S(1,1)=X5
  193.             S(2,1)=SIGN(X5,X1)
  194.             S(1,2)=-S(2,1)
  195.             S(2,2)=X5
  196.          ELSE
  197. *
  198.             S1=T1/XT
  199.             S2=T2/XT
  200. *
  201. *       IF (I.EQ.1) THEN
  202. *            write(6,*) 'Cas 1'
  203.             S(1,1)= S2
  204.             S(2,1)=-S1
  205.             S(1,2)= S1
  206.             S(2,2)= S2
  207. *         ELSE
  208. *            write(6,*) 'Cas 2'
  209. *            S(1,1)= S1
  210. *            S(2,1)= S2
  211. *            S(1,2)=-S2
  212. *            S(2,2)= S1
  213. *         ENDIF
  214.          ENDIF
  215.       endif
  216.  1000 continue
  217. * Scaling des valeurs propres
  218.       d(1)=d(1)*amax
  219.       d(2)=d(2)*amax
  220.       if (lverif) then
  221. * Orthonormalite
  222.          call vmorth(s,2,3,xpre,r)
  223.          if (ierr.ne.0) then
  224.             write(6,*) '!!! X non orthonormale'
  225.             if (ldbg) then
  226.                write(6,*) 'Matrice des VPs s='
  227.                do ii=1,2
  228.                   write (6,*) (s(ii,jj),jj=1,2)
  229.                enddo
  230.                write(6,*) 'Matrice des Pscals R='
  231.                do ii=1,2
  232.                   write (6,*) (R(ii,jj),jj=1,2)
  233.                enddo
  234.             endif
  235.             goto 9999
  236.          endif
  237. * Vecteur propre
  238.          xprevp=amax*xprev
  239.          do i=1,2
  240. *     vp(i) est-il bien le vp associe a la valeur propre d(i)
  241. *     xnr espace de stockage ici
  242.             call vvectp(a,2,3,d(i),s(1,i),xprevp,xnr)
  243.             if (ierr.ne.0) then
  244.                write(6,*) '!!!   vp(i) nest pas vecteur propre i=',i
  245.                if (ldbg) then
  246.                   write(6,*) ' v= (S - D(i)) s(j,i)  ='
  247.                   write(6,*) (xnr(jj),jj=1,2)
  248.                endif
  249.                goto 9999
  250.             endif
  251.          enddo
  252.       endif
  253.       RETURN
  254. *
  255. * Error handling
  256. *
  257.  
  258.  9999 CONTINUE
  259.       MOTERR(1:8)='JACOB2'
  260.       call erreur(1127)
  261.       RETURN
  262.       END
  263.  
  264.  

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