* DSYGVD Example Program Text * NAG Copyright 2005. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NB, NMAX PARAMETER (NB=64,NMAX=10) INTEGER LDA, LDB, LIWORK, LWORK PARAMETER (LDA=NMAX,LDB=NMAX,LIWORK=3+5*NMAX, + LWORK=1+(NB+6+2*NMAX)*NMAX) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D+0) * .. Local Scalars .. DOUBLE PRECISION ANORM, BNORM, EPS, RCOND, RCONDB, T1, T2, T3 INTEGER I, IFAIL, INFO, J, LIWOPT, LWOPT, N * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), B(LDB,NMAX), EERBND(NMAX), + RCONDZ(NMAX), W(NMAX), WORK(LWORK), ZERBND(NMAX) INTEGER IWORK(LIWORK) * .. External Functions .. DOUBLE PRECISION DLAMCH, DLANSY EXTERNAL DLAMCH, DLANSY * .. External Subroutines .. EXTERNAL DDISNA, DSYGVD, DTRCON, X04CAF * .. Intrinsic Functions .. INTRINSIC ABS * .. Executable Statements .. WRITE (NOUT,*) 'DSYGVD Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read the upper triangular parts of the matrices A and B * READ (NIN,*) ((A(I,J),J=I,N),I=1,N) READ (NIN,*) ((B(I,J),J=I,N),I=1,N) * * Compute the one-norms of the symmetric matrices A and B * ANORM = DLANSY('One norm','Upper',N,A,LDA,WORK) BNORM = DLANSY('One norm','Upper',N,B,LDB,WORK) * * Solve the generalized symmetric eigenvalue problem * A*B*x = lambda*x (ITYPE = 2) * CALL DSYGVD(2,'Vectors','Upper',N,A,LDA,B,LDB,W,WORK,LWORK, + IWORK,LIWORK,INFO) LWOPT = WORK(1) LIWOPT = IWORK(1) * IF (INFO.EQ.0) THEN * * Print solution * WRITE (NOUT,*) 'Eigenvalues' WRITE (NOUT,99999) (W(J),J=1,N) * IFAIL = 0 CALL X04CAF('General',' ',N,N,A,LDA,'Eigenvectors',IFAIL) * * Call DTRCON (F07TGF) to estimate the reciprocal condition * number of the Cholesky factor of B. Note that: * cond(B) = 1/RCOND**2 * CALL DTRCON('One norm','Upper','Non-unit',N,B,LDB,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 and * eigenvectors * EPS = DLAMCH('Eps') IF (RCOND.GE.EPS) THEN * * Call DDISNA (F08FLF) to estimate reciprocal condition * numbers for the eigenvectors of (A*B - lambda*I) * CALL DDISNA('Eigenvectors',N,N,W,RCONDZ,INFO) * * Compute the error estimates for the eigenvalues and * eigenvectors * T1 = ONE/RCOND T2 = EPS*T1 T3 = ANORM*BNORM DO 20 I = 1, N EERBND(I) = EPS*(T3+ABS(W(I))/RCONDB) ZERBND(I) = T2*(T3/RCONDZ(I)+T1) 20 CONTINUE * * Print the approximate error bounds for the eigenvalues * and vectors * WRITE (NOUT,*) WRITE (NOUT,*) 'Error estimates for the eigenvalues' WRITE (NOUT,99998) (EERBND(I),I=1,N) WRITE (NOUT,*) WRITE (NOUT,*) 'Error estimates for the eigenvectors' WRITE (NOUT,99998) (ZERBND(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 DSYGVD. INFO =', INFO END IF * * Print workspace information * IF (LWORK.LT.LWOPT) THEN WRITE (NOUT,*) WRITE (NOUT,99995) 'Optimum workspace required = ', LWOPT, + 'Workspace provided = ', LWORK END IF IF (LIWORK.LT.LIWOPT) THEN WRITE (NOUT,*) WRITE (NOUT,99995) 'Integer workspace required = ', LIWOPT, + 'Integer workspace provided = ', LIWORK END IF ELSE WRITE (NOUT,*) 'NMAX too small' END IF STOP * 99999 FORMAT (3X,(6F11.4)) 99998 FORMAT (4X,1P,6E11.1) 99997 FORMAT (1X,A,I4,A) 99996 FORMAT (1X,A,I4) 99995 FORMAT (1X,A,I5,/1X,A,I5) END