* F08TAF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=10) CHARACTER UPLO PARAMETER (UPLO='U') * .. Local Scalars .. DOUBLE PRECISION ANORM, BNORM, EPS, RCOND, RCONDB, T1, T2 INTEGER I, INFO, J, N * .. Local Arrays .. DOUBLE PRECISION AP((NMAX*(NMAX+1))/2), BP((NMAX*(NMAX+1))/2), + DUMMY(1,1), EERBND(NMAX), W(NMAX), WORK(3*NMAX) INTEGER IWORK(NMAX) * .. External Functions .. DOUBLE PRECISION F06RDF, X02AJF EXTERNAL F06RDF, X02AJF * .. External Subroutines .. EXTERNAL DSPGV, DTPCON * .. Intrinsic Functions .. INTRINSIC ABS * .. Executable Statements .. WRITE (NOUT,*) 'F08TAF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read the upper or lower triangular parts of the matrices A and * B from data file * IF (UPLO.EQ.'U') THEN READ (NIN,*) ((AP(I+(J*(J-1))/2),J=I,N),I=1,N) READ (NIN,*) ((BP(I+(J*(J-1))/2),J=I,N),I=1,N) ELSE IF (UPLO.EQ.'L') THEN READ (NIN,*) ((AP(I+((2*N-J)*(J-1))/2),J=1,I),I=1,N) READ (NIN,*) ((BP(I+((2*N-J)*(J-1))/2),J=1,I),I=1,N) END IF * * Compute the one-norms of the symmetric matrices A and B * ANORM = F06RDF('One norm',UPLO,N,AP,WORK) BNORM = F06RDF('One norm',UPLO,N,BP,WORK) * * Solve the generalized symmetric eigenvalue problem * A*x = lambda*B*x (ITYPE = 1) * CALL DSPGV(1,'No vectors',UPLO,N,AP,BP,W,DUMMY,1,WORK,INFO) * IF (INFO.EQ.0) THEN * * Print solution * WRITE (NOUT,*) 'Eigenvalues' WRITE (NOUT,99999) (W(J),J=1,N) * * Call DTPCON (F07UGF) to estimate the reciprocal condition * number of the Cholesky factor of B. Note that: * cond(B) = 1/RCOND**2 * CALL DTPCON('One norm',UPLO,'Non-unit',N,BP,RCOND,WORK, + IWORK,INFO) * * Print the reciprocal condition number of B * RCONDB = RCOND**2 WRITE (NOUT,*) WRITE (NOUT,*) + 'Estimate of reciprocal condition number for B' WRITE (NOUT,99998) RCONDB * * Get the machine precision, EPS, and if RCONDB is not less * than EPS**2, compute error estimates for the eigenvalues * EPS = X02AJF() IF (RCOND.GE.EPS) THEN T1 = EPS/RCONDB T2 = ANORM/BNORM DO 20 I = 1, N EERBND(I) = T1*(T2+ABS(W(I))) 20 CONTINUE * * Print the approximate error bounds for the eigenvalues * WRITE (NOUT,*) WRITE (NOUT,*) 'Error estimates for the eigenvalues' WRITE (NOUT,99998) (EERBND(I),I=1,N) ELSE WRITE (NOUT,*) WRITE (NOUT,*) 'B is very ill-conditioned, error ', + 'estimates have not been computed' END IF ELSE IF (INFO.GT.N .AND. INFO.LE.2*N) THEN I = INFO - N WRITE (NOUT,99997) 'The leading minor of order ', I, + ' of B is not positive definite' ELSE WRITE (NOUT,99996) 'Failure in DSPGV. INFO =', INFO END IF ELSE WRITE (NOUT,*) 'NMAX too small' END IF * 99999 FORMAT (3X,(6F11.4)) 99998 FORMAT (4X,1P,6E11.1) 99997 FORMAT (1X,A,I4,A) 99996 FORMAT (1X,A,I4) END