* F08YYF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=8) INTEGER LDA, LDB, LDVL, LDVR, LWORK PARAMETER (LDA=NMAX,LDB=NMAX,LDVL=NMAX,LDVR=NMAX, + LWORK=2*NMAX*NMAX) * .. Local Scalars .. DOUBLE PRECISION EPS, SNORM, STNRM, TNORM INTEGER I, INFO, J, M, N * .. Local Arrays .. COMPLEX *16 A(LDA,NMAX), B(LDB,NMAX), VL(LDVL,NMAX), + VR(LDVR,NMAX), WORK(LWORK) DOUBLE PRECISION DIF(NMAX), RWORK(2*NMAX), S(NMAX) INTEGER IWORK(NMAX+2) LOGICAL SELECT(1) * .. External Functions .. DOUBLE PRECISION F06BNF, F06UAF, X02AJF EXTERNAL F06BNF, F06UAF, X02AJF * .. External Subroutines .. EXTERNAL ZTGEVC, ZTGSNA * .. Executable Statements .. WRITE (NOUT,*) 'F08YYF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read A and B from data file * READ (NIN,*) ((A(I,J),J=1,N),I=1,N) READ (NIN,*) ((B(I,J),J=1,N),I=1,N) * * Calculate the left and right generalized eigenvectors of the * pair (A,B). * CALL ZTGEVC('Both','All',SELECT,N,A,LDA,B,LDB,VL,LDVL,VR,LDVR, + N,M,WORK,RWORK,INFO) * * Estimate condition numbers for all the generalized eigenvalues * and right eigenvectors of the pair (A,B) * CALL ZTGSNA('Both','All',SELECT,N,A,LDA,B,LDB,VL,LDVL,VR,LDVR, + S,DIF,N,M,WORK,LWORK,IWORK,INFO) * * Print condition numbers of eigenvalues and right eigenvectors * WRITE (NOUT,*) 'S' WRITE (NOUT,99999) (S(I),I=1,M) WRITE (NOUT,*) WRITE (NOUT,*) 'DIF' WRITE (NOUT,99999) (DIF(I),I=1,M) * * Calculate approximate error estimates * * Compute the 1-norms of A and B and then compute * SQRT(SNORM**2 + TNORM**2) * EPS = X02AJF() SNORM = F06UAF('1-norm',N,N,A,LDA,RWORK) TNORM = F06UAF('1-norm',N,N,B,LDB,RWORK) STNRM = F06BNF(SNORM,TNORM) WRITE (NOUT,*) WRITE (NOUT,*) + 'Approximate error estimates for eigenvalues of (A,B)' WRITE (NOUT,99999) (EPS*STNRM/S(I),I=1,M) WRITE (NOUT,*) WRITE (NOUT,*) + 'Approximate error estimates for right eigenvectors of (A,B)' WRITE (NOUT,99999) (EPS*STNRM/DIF(I),I=1,M) ELSE WRITE (NOUT,*) 'NMAX too small' END IF * 99999 FORMAT ((3X,1P,7E11.1)) END