c program DRSEVVUN c>> 1996-05-28 DRSEVVUN Krogh Added external statement. c>> 1994-10-19 DRSEVVUN Krogh Changes to use M77CON c>> 1994-09-22 DRSEVVUN CLL c>> 1992-04-23 CLL c>> 1992-03-04 DRSEVVUN Krogh Initial version. c Demonstrate unsymmetric eigenvalue/eigenvector subroutine SEVVUN. c ------------------------------------------------------------------ c--S replaces "?": DR?EVVUN, ?EVVUN, ?VECP, ?MATP, ?DOT c ------------------------------------------------------------------ integer I, J, LDA, N parameter (LDA = 3) integer IFLAG(LDA) real A(LDA,LDA), VR(LDA), VI(LDA), VEC(LDA, LDA) real WORK(LDA), ASAV(LDA, LDA), D(LDA, LDA), ANORM external SDOT real SDOT data (A(1,I), I=1,3) / 8.0e0, -1.0e0, -5.0e0 / data (A(2,I), I=1,3) / -4.0e0, 4.0e0, -2.0e0 / data (A(3,I), I=1,3) / 18.0e0, -5.0e0, -7.0e0 / data ANORM / 30.0e0 / data N / LDA / c ------------------------------------------------------------------ print*,'DRSEVVUN.. Demo driver for SEVVUN.' c c First copy A() to ASAV() for later residual check. c do 20 J = 1, N do 10 I = 1, N ASAV(I, J) = A(I, J) 10 continue 20 continue call SEVVUN(A(1, 1), LDA, N, VR, VI, VEC, IFLAG, WORK) print '(a, I2)',' IFLAG(1) =',IFLAG(1) if (IFLAG(1) .le. 2) then call SVECP(VR, N, ' Real part of the eigenvalues') call SVECP(VI, N, ' Imaginary part of the eigenvalues') call SMATP(VEC, LDA, N, N, * '0 Eigenvectors as columns or pairs of columns') c c As a check compute D = (ASAV*VEC - VEC*EVAL) / ANORM. c Expect D to be close to the machine precision. c do 50 J = 1, N if (VI(J) .eq. 0.0e0) then c Compute residual for a real eigenvalue and eigenvector do 30 I = 1, N D(I, J) = (SDOT(N, ASAV(I,1), LDA, VEC(1,J), 1) - * VEC(I,J) * VR(J)) / ANORM 30 continue else if (VI(J) .gt. 0.0e0) then do 40 I = 1, N D(I, J) = (SDOT(N,ASAV(I,1),LDA,VEC(1,J),1) - * VEC(I,J)*VR(J)+VEC(I,J+1)*VI(J))/ANORM D(I, J+1) = (SDOT(N,ASAV(I,1),LDA,VEC(1,J+1),1) - * VEC(I,J)*VI(J)-VEC(I,J+1)*VR(J))/ANORM 40 continue end if 50 continue call SMATP(D, LDA, N, N, * '0 Packed residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM') else print'(/a)', ' Failure in SEVVUN. ' end if stop end