* F08NBF 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, LDVL, LDVR, LWORK PARAMETER (LDA=NMAX,LDVL=NMAX,LDVR=NMAX,LWORK=(2+NB)*NMAX) * .. Local Scalars .. COMPLEX *16 EIG DOUBLE PRECISION ABNRM, EPS, ERBND, RCND, TOL INTEGER I, IHI, ILO, INFO, J, LWKOPT, N * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), RCONDE(NMAX), RCONDV(NMAX), + SCALE(NMAX), VL(LDVL,NMAX), VR(LDVR,NMAX), + WI(NMAX), WORK(LWORK), WR(NMAX) INTEGER IWORK(2*NMAX-2) * .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF * .. External Subroutines .. EXTERNAL DGEEVX * .. Intrinsic Functions .. INTRINSIC CMPLX * .. Executable Statements .. WRITE (NOUT,*) 'F08NBF Example Program Results' * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read the matrix A from data file * READ (NIN,*) ((A(I,J),J=1,N),I=1,N) * * Solve the eigenvalue problem * CALL DGEEVX('Balance','Vectors (left)','Vectors (right)', + 'Both reciprocal condition numbers',N,A,LDA,WR,WI, + VL,LDVL,VR,LDVR,ILO,IHI,SCALE,ABNRM,RCONDE,RCONDV, + WORK,LWORK,IWORK,INFO) * IF (INFO.EQ.0) THEN * * Compute the machine precision * EPS = X02AJF() TOL = EPS*ABNRM * * Print the eigenvalues and vectors, and associated condition * number and bounds * DO 20 J = 1, N * * Print information on jth eigenvalue * WRITE (NOUT,*) IF (WI(J).EQ.0.0D0) THEN WRITE (NOUT,99999) 'Eigenvalue(', J, ') = ', WR(J) ELSE EIG = CMPLX(WR(J),WI(J),KIND=KIND(EPS)) WRITE (NOUT,99998) 'Eigenvalue(', J, ') = ', EIG END IF RCND = RCONDE(J) WRITE (NOUT,*) WRITE (NOUT,99997) 'Reciprocal condition number = ', RCND IF (RCND.GT.0.0D0) THEN ERBND = TOL/RCND WRITE (NOUT,99997) 'Error bound = ', + ERBND ELSE WRITE (NOUT,*) 'Error bound is infinite' END IF * * Print information on jth eigenvector * WRITE (NOUT,*) WRITE (NOUT,99996) 'Eigenvector(', J, ')' IF (WI(J).EQ.0.0D0) THEN WRITE (NOUT,99995) (VR(I,J),I=1,N) ELSE IF (WI(J).GT.0.0D0) THEN WRITE (NOUT,99994) (VR(I,J),VR(I,J+1),I=1,N) ELSE WRITE (NOUT,99994) (VR(I,J-1),-VR(I,J),I=1,N) END IF RCND = RCONDV(J) WRITE (NOUT,*) WRITE (NOUT,99997) 'Reciprocal condition number = ', RCND IF (RCND.GT.0.0D0) THEN ERBND = TOL/RCND WRITE (NOUT,99997) 'Error bound = ', + ERBND ELSE WRITE (NOUT,*) 'Error bound is infinite' END IF 20 CONTINUE ELSE WRITE (NOUT,*) WRITE (NOUT,99993) 'Failure in DGEEVX. INFO = ', INFO END IF * * Print workspace information * LWKOPT = WORK(1) IF (LWORK.LT.LWKOPT) THEN WRITE (NOUT,*) WRITE (NOUT,99992) 'Optimum workspace required = ', LWKOPT, + 'Workspace provided = ', LWORK END IF ELSE WRITE (NOUT,*) WRITE (NOUT,*) 'NMAX too small' END IF STOP * 99999 FORMAT (1X,A,I2,A,1P,E11.4) 99998 FORMAT (1X,A,I2,A,'(',1P,E11.4,',',1P,E11.4,')') 99997 FORMAT (1X,A,1P,E8.1) 99996 FORMAT (1X,A,I2,A) 99995 FORMAT (1X,1P,E11.4) 99994 FORMAT (1X,'(',1P,E11.4,',',1P,E11.4,')') 99993 FORMAT (1X,A,I4) 99992 FORMAT (1X,A,I5,/1X,A,I5) END