*     F08SAF Example Program Text
*     Mark 21 Release.  NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NB, NMAX
      PARAMETER        (NB=64,NMAX=10)
      INTEGER          LDA, LDB, LWORK
      PARAMETER        (LDA=NMAX,LDB=NMAX,LWORK=(NB+2)*NMAX)
*     .. Local Scalars ..
      DOUBLE PRECISION ANORM, BNORM, EPS, RCOND, RCONDB, T1, T2, T3
      INTEGER          I, IFAIL, INFO, J, LWKOPT, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,NMAX), EERBND(NMAX),
     +                 RCONDZ(NMAX), W(NMAX), WORK(LWORK), ZERBND(NMAX)
      INTEGER          IWORK(NMAX)
*     .. External Functions ..
      DOUBLE PRECISION F06RCF, X02AJF
      EXTERNAL         F06RCF, X02AJF
*     .. External Subroutines ..
      EXTERNAL         DDISNA, DSYGV, DTRCON, X04CAF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS, INT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08SAF 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 = F06RCF('One norm','Upper',N,A,LDA,WORK)
         BNORM = F06RCF('One norm','Upper',N,B,LDB,WORK)
*
*        Solve the generalized symmetric eigenvalue problem
*        A*x = lambda*B*x (ITYPE = 1)
*
         CALL DSYGV(1,'Vectors','Upper',N,A,LDA,B,LDB,W,WORK,LWORK,INFO)
         LWKOPT = INT(WORK(1))
*
         IF (INFO.EQ.0) THEN
*
*           Print solution
*
            WRITE (NOUT,*) 'Eigenvalues'
            WRITE (NOUT,99999) (W(J),J=1,N)
*
            WRITE (NOUT,*)
            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 = X02AJF()
            IF (RCOND.GE.EPS) THEN
*
*              Call DDISNA (F08FLF) to estimate reciprocal condition
*              numbers for the eigenvectors of (A - lambda*B)
*
               CALL DDISNA('Eigenvectors',N,N,W,RCONDZ,INFO)
*
*              Compute the error estimates for the eigenvalues and
*              eigenvectors
*
               T1 = EPS/RCONDB
               T2 = ANORM/BNORM
               T3 = T2/RCOND
               DO 20 I = 1, N
                  EERBND(I) = T1*(T2+ABS(W(I)))
                  ZERBND(I) = T1*(T3+ABS(W(I)))/RCONDZ(I)
   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 DSYGV. INFO =', INFO
         END IF
*
*        Print workspace information
*
         IF (LWORK.LT.LWKOPT) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99995) 'Optimum workspace required = ', LWKOPT,
     +        'Workspace provided         = ', LWORK
         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)
99995 FORMAT (1X,A,I5,/1X,A,I5)
      END