C JACOB2 SOURCE GOUNAND 26/01/16 21:15:05 12450 SUBROUTINE JACOB2(A,D,S) C===================================================================== C C DIAGONALISATION DE A SYMETRIQUE C ENTREEES C A(3,3) = MATRICE A DIAGONALISER A(1,3)=A(2,3)=0. C A(3,1)=A(3,2)=0. C SORTIES C D(3) = LES 3 VALEURS PROPRES D(1) > D(2) ET D(3)=A(3,3) C S(3,3) = VECTEURS PROPRES S(I,3)=0. 0. 1. C C RECUPERATION INCA FEVRIER 85 EBERSOLT C==================================================================== IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C -INC PPARAM -INC CCOPTIO -INC CCREEL C DIMENSION A(3,*),D(*),S(3,*) LOGICAL lverif,ldbg DIMENSION r(3,3),xnr(3) C C Valeur particuliere : SQRT(2.)*0.5 = 1/SQRT(2.) PARAMETER ( X1sr2 = 0.70710678118654752440084436210484904D0 ) * Verifications lverif=.false. * Impressions ldbg=.false. xpre=xzprec*3.D0 xpet=xpetit*10.D0 xprev=xzprec*10.D0 * D(3)=A(3,3) S(3,1)=0.D0 S(3,2)=0.D0 S(1,3)=0.D0 S(2,3)=0.D0 S(3,3)=1.D0 * Scaling de A amax=0.d0 do i=1,2 do j=1,2 amax=max(amax,abs(a(j,i))) enddo enddo * Matrice nulle if (amax.lt.xpet) then D(1)=XZERO D(2)=XZERO S(1,1)=1.D0 S(2,1)=XZERO S(1,2)=XZERO S(2,2)=1.D0 return endif C amax1=1.D0/amax do i=1,2 do j=1,2 R(j,i)=a(j,i)*amax1 enddo enddo * Trace J1(A) et determinant J2(A) de A xj1=R(1,1)+R(2,2) xj2=R(1,1)*R(2,2)-R(1,2)**2 * Le discriminant reduit J4(A) du polynome caracteristique en l : * x^2 - J1(A) x + J2(A) = 0 * est : J4(A) = (J1(A) / 2) ^ 2 - J2(A) * * Construisons S le deviateur de A : S = A - J1(A)/2 I * On a J1(S)=0 et lambda_A = lambda_S + J1(A)/2 * * A et S ont les memes vecteurs propres * Le polynome caracteristique de S est : * x^2 + J2(S) = 0 * avec J2(S) = - [ ((a11-a22)/2)^2 + a12^2 ] = -J4(A) * Cette expression de J4 est interessante car somme de carres donc positive * xj4 = ((R(1,1)-R(2,2))/2)**2+R(1,2)**2 xsj4=sqrt(xj4) stra=sign(1.d0,xj1) * Formules de Fagnano modifiees pour les deux racines * cf. The Ins and Outs of Solving Quadratic Equations with Floating-Point * Arithmetic, Frederic Goualard, HAL preprint x1=(xj1/2)+stra*xsj4 if (x1.eq.0.d0) then x2=0.d0 else x2=xj2/x1 endif if (stra.gt.0) then d(1)=x1 d(2)=x2 else d(1)=x2 d(2)=x1 endif x3=xsj4 * * Ancienne methode * if (.false.) then X1 =2.D0*R(1,2) X2 =R(1,1)-R(2,2) X3 =SQRT(X2*X2+X1*X1) D(1)=0.5D0*(R(1,1)+R(2,2)+X3) D(2)=D(1)-X3 endif * Raccourci lorsque les deux valeurs propres sont egales * xsca=max(xpre,xpet) if (ldbg) then write(6,*) 'xpre,xpet,xsca,x3=',xpre,xpet,xsca,x3 endif * if (x3.lt.xsca) then if (ldbg) write(6,*) '2 valeurs propres egales detectees' do i=1,2 do j=1,2 if (i.eq.j) then s(j,i)=1.d0 else s(j,i)=0.d0 endif enddo enddo goto 1000 endif C C Appliquer la methode de Scherzinger et Dohrmann CMAME'08 * Recherche du plus grand vecteur de l'image de A - d(1) I : C do i=1,2 R(i,i)=R(i,i)-d(1) enddo imax=0 xmax=0.d0 do i=1,2 xnr(i)=0.d0 do j=1,2 xnr(i)=xnr(i)+r(j,i)**2 enddo if (xnr(i).gt.xmax) then xmax=xnr(i) imax=i endif enddo igr=imax * Normalement pas possible si valeurs propres non egales if (igr.eq.0) then igr=1 r(1,1)=1.d0 r(2,1)=0.d0 else xnorm=sqrt(xnr(igr)) do i=1,2 r(i,igr)=r(i,igr)/xnorm enddo endif S(1,2)= r(1,igr) S(2,2)= r(2,igr) S(1,1)= S(2,2) S(2,1)=-S(1,2) if (.false.) then * Methode ancienne R11=A(1,1)-D(1) R21=A(1,2) XNR1=SQRT(R11*R11+R21*R21) R12=A(2,1) R22=A(2,2)-D(1) XNR2=SQRT(R12*R12+R22*R22) IF (XNR1.GE.XNR2) THEN I=1 T1=R11 T2=R21 XT=XNR1 ELSE I=2 T1=R12 T2=R22 XT=XNR2 ENDIF C XSCA=MAX(MAX(ABS(D(1)),ABS(D(2)))*XZPREC,XPETIT)*2 IF (XT.LE.XSCA) THEN * write(6,*) 'Cas 0' X5 = X1sr2 S(1,1)=X5 S(2,1)=SIGN(X5,X1) S(1,2)=-S(2,1) S(2,2)=X5 ELSE * S1=T1/XT S2=T2/XT * * IF (I.EQ.1) THEN * write(6,*) 'Cas 1' S(1,1)= S2 S(2,1)=-S1 S(1,2)= S1 S(2,2)= S2 * ELSE * write(6,*) 'Cas 2' * S(1,1)= S1 * S(2,1)= S2 * S(1,2)=-S2 * S(2,2)= S1 * ENDIF ENDIF endif 1000 continue * Scaling des valeurs propres d(1)=d(1)*amax d(2)=d(2)*amax if (lverif) then * Orthonormalite call vmorth(s,2,3,xpre,r) if (ierr.ne.0) then write(6,*) '!!! X non orthonormale' if (ldbg) then write(6,*) 'Matrice des VPs s=' do ii=1,2 write (6,*) (s(ii,jj),jj=1,2) enddo write(6,*) 'Matrice des Pscals R=' do ii=1,2 write (6,*) (R(ii,jj),jj=1,2) enddo endif goto 9999 endif * Vecteur propre xprevp=amax*xprev do i=1,2 * vp(i) est-il bien le vp associe a la valeur propre d(i) * xnr espace de stockage ici call vvectp(a,2,3,d(i),s(1,i),xprevp,xnr) if (ierr.ne.0) then write(6,*) '!!! vp(i) nest pas vecteur propre i=',i if (ldbg) then write(6,*) ' v= (S - D(i)) s(j,i) =' write(6,*) (xnr(jj),jj=1,2) endif goto 9999 endif enddo endif RETURN * * Error handling * 9999 CONTINUE MOTERR(1:8)='JACOB2' call erreur(1127) RETURN END