* 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 LDS, LDT, LDVL, LDVR, LWORK PARAMETER (LDS=NMAX,LDT=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 S(LDS,NMAX), T(LDT,NMAX), VL(LDVL,NMAX), + VR(LDVR,NMAX), WORK(LWORK) DOUBLE PRECISION DIF(NMAX), RWORK(2*NMAX), SCON(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 S and T from data file * READ (NIN,*) ((S(I,J),J=1,N),I=1,N) READ (NIN,*) ((T(I,J),J=1,N),I=1,N) * * Calculate the left and right generalized eigenvectors of the * pair (S,T). * CALL ZTGEVC('Both','All',SELECT,N,S,LDS,T,LDT,VL,LDVL,VR,LDVR, + N,M,WORK,RWORK,INFO) * * Estimate condition numbers for all the generalized eigenvalues * and right eigenvectors of the pair (S,T) * CALL ZTGSNA('Both','All',SELECT,N,S,LDS,T,LDT,VL,LDVL,VR,LDVR, + SCON,DIF,N,M,WORK,LWORK,IWORK,INFO) * * Print condition numbers of eigenvalues and right eigenvectors * WRITE (NOUT,*) 'SCON' WRITE (NOUT,99999) (SCON(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 S and T and then compute * SQRT(SNORM**2 + TNORM**2) * EPS = X02AJF() SNORM = F06UAF('1-norm',N,N,S,LDS,RWORK) TNORM = F06UAF('1-norm',N,N,T,LDT,RWORK) STNRM = F06BNF(SNORM,TNORM) WRITE (NOUT,*) WRITE (NOUT,*) + 'Approximate error estimates for eigenvalues of (S,T)' WRITE (NOUT,99999) (EPS*TNORM/SCON(I),I=1,M) WRITE (NOUT,*) WRITE (NOUT,*) + 'Approximate error estimates for right eigenvectors of (S,T)' WRITE (NOUT,99999) (EPS*TNORM/DIF(I),I=1,M) ELSE WRITE (NOUT,*) 'NMAX too small' END IF STOP * 99999 FORMAT ((3X,1P,7E11.1)) END