* F08JAF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=20) INTEGER LDZ PARAMETER (LDZ=NMAX) * .. Local Scalars .. DOUBLE PRECISION EERRBD, EPS INTEGER I, IFAIL, INFO, N * .. Local Arrays .. DOUBLE PRECISION D(NMAX), E(NMAX-1), RCONDZ(NMAX), WORK(2*NMAX-2), + Z(LDZ,NMAX), ZERRBD(NMAX) * .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF * .. External Subroutines .. EXTERNAL DDISNA, DSTEV, X04CAF * .. Intrinsic Functions .. INTRINSIC ABS, MAX * .. Executable Statements .. WRITE (NOUT,*) 'F08JAF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read the diagonal and off-diagonal elements of the matrix A * from data file * READ (NIN,*) (D(I),I=1,N) READ (NIN,*) (E(I),I=1,N-1) * * Solve the symmetric tridiagonal eigenvalue problem * CALL DSTEV('Vectors',N,D,E,Z,LDZ,WORK,INFO) * IF (INFO.EQ.0) THEN * * Print solution * WRITE (NOUT,*) 'Eigenvalues' WRITE (NOUT,99999) (D(I),I=1,N) * IFAIL = 0 CALL X04CAF('General',' ',N,N,Z,LDZ,'Eigenvectors',IFAIL) * * Get the machine precision, EPS and compute the approximate * error bound for the computed eigenvalues. Note that for * the 2-norm, max( abs(D(i)) ) = norm(A), and since the * eigenvalues are returned in ascending order * max( abs(D(i)) ) = max( abs(D(1)), abs(D(n))) * EPS = X02AJF() EERRBD = EPS*MAX(ABS(D(1)),ABS(D(N))) * * Call DDISNA (F08FLF) to estimate reciprocal condition * numbers for the eigenvectors * CALL DDISNA('Eigenvectors',N,N,D,RCONDZ,INFO) * * Compute the error estimates for the eigenvectors * DO 20 I = 1, N ZERRBD(I) = EERRBD/RCONDZ(I) 20 CONTINUE * * Print the approximate error bounds for the eigenvalues * and vectors * WRITE (NOUT,*) WRITE (NOUT,*) 'Error estimate for the eigenvalues' WRITE (NOUT,99998) EERRBD WRITE (NOUT,*) WRITE (NOUT,*) 'Error estimates for the eigenvectors' WRITE (NOUT,99998) (ZERRBD(I),I=1,N) ELSE WRITE (NOUT,99997) 'Failure in DSTEV. INFO =', INFO END IF ELSE WRITE (NOUT,*) 'NMAX too small' END IF * 99999 FORMAT (3X,(8F8.4)) 99998 FORMAT (4X,1P,6E11.1) 99997 FORMAT (1X,A,I4) END