* F08WBF 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, LDVR, LWORK PARAMETER (LDA=NMAX,LDB=NMAX,LDVR=NMAX, + LWORK=NMAX*NB+2*NMAX*NMAX) * .. Local Scalars .. COMPLEX *16 EIG DOUBLE PRECISION ABNORM, ABNRM, BBNRM, EPS, ERBND, RCND, SMALL, + TOL INTEGER I, IHI, ILO, INFO, J, LWKOPT, N * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX), + B(LDB,NMAX), BETA(NMAX), DUMMY(1,1), + LSCALE(NMAX), RCONDE(NMAX), RCONDV(NMAX), + RSCALE(NMAX), VR(LDVR,NMAX), WORK(LWORK) INTEGER IWORK(NMAX+6) LOGICAL BWORK(NMAX) * .. External Functions .. DOUBLE PRECISION F06BNF, X02AJF, X02AMF EXTERNAL F06BNF, X02AJF, X02AMF * .. External Subroutines .. EXTERNAL DGGEVX * .. Intrinsic Functions .. INTRINSIC ABS, CMPLX, INT * .. Executable Statements .. WRITE (NOUT,*) 'F08WBF Example Program Results' * Skip heading in data file READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * Read in the matrices A and B * READ (NIN,*) ((A(I,J),J=1,N),I=1,N) READ (NIN,*) ((B(I,J),J=1,N),I=1,N) * * Solve the generalized eigenvalue problem * CALL DGGEVX('Balance','No vectors (left)','Vectors (right)', + 'Both reciprocal condition numbers',N,A,LDA,B,LDB, + ALPHAR,ALPHAI,BETA,DUMMY,1,VR,LDVR,ILO,IHI,LSCALE, + RSCALE,ABNRM,BBNRM,RCONDE,RCONDV,WORK,LWORK,IWORK, + BWORK,INFO) * IF (INFO.GT.0) THEN WRITE (NOUT,*) WRITE (NOUT,99999) 'Failure in DGGEVX. INFO =', INFO ELSE * * Compute the machine precision, the safe range parameter * SMALL and sqrt(ABNRM**2+BBNRM**2) * EPS = X02AJF() SMALL = X02AMF() ABNORM = F06BNF(ABNRM,BBNRM) TOL = EPS*ABNORM * * Print out eigenvalues and vectors and associated condition * number and bounds * DO 20 J = 1, N * * Print out information on the jth eigenvalue * WRITE (NOUT,*) IF ((ABS(ALPHAR(J))+ABS(ALPHAI(J)))*SMALL.GE.ABS(BETA(J)) + ) THEN WRITE (NOUT,99998) 'Eigenvalue(', J, ')', + ' is numerically infinite or undetermined', + 'ALPHAR(', J, ') = ', ALPHAR(J), ', ALPHAI(', J, + ') = ', ALPHAI(J), ', BETA(', J, ') = ', BETA(J) ELSE IF (ALPHAI(J).EQ.0.0D0) THEN WRITE (NOUT,99997) 'Eigenvalue(', J, ') = ', + ALPHAR(J)/BETA(J) ELSE EIG = CMPLX(ALPHAR(J),ALPHAI(J),KIND=KIND(A)) + /BETA(J) WRITE (NOUT,99996) 'Eigenvalue(', J, ') = ', EIG END IF END IF RCND = RCONDE(J) WRITE (NOUT,*) WRITE (NOUT,99995) 'Reciprocal condition number = ', RCND IF (RCND.GT.0.0D0) THEN ERBND = TOL/RCND WRITE (NOUT,99995) 'Error bound = ', + ERBND ELSE WRITE (NOUT,*) 'Error bound is infinite' END IF * * Print out information on the jth eigenvector * WRITE (NOUT,*) WRITE (NOUT,99994) 'Eigenvector(', J, ')' IF (ALPHAI(J).EQ.0.0D0) THEN WRITE (NOUT,99993) (VR(I,J),I=1,N) ELSE IF (ALPHAI(J).GT.0.0D0) THEN WRITE (NOUT,99992) (VR(I,J),VR(I,J+1),I=1,N) ELSE WRITE (NOUT,99992) (VR(I,J-1),-VR(I,J),I=1,N) END IF RCND = RCONDV(J) WRITE (NOUT,*) WRITE (NOUT,99995) 'Reciprocal condition number = ', RCND IF (RCND.GT.0.0D0) THEN ERBND = TOL/RCND WRITE (NOUT,99995) 'Error bound = ', + ERBND ELSE WRITE (NOUT,*) 'Error bound is infinite' END IF 20 CONTINUE * LWKOPT = INT(WORK(1)) IF (LWORK.LT.LWKOPT) THEN WRITE (NOUT,*) WRITE (NOUT,99991) 'Optimum workspace required = ', + LWKOPT, 'Workspace provided = ', LWORK END IF END IF ELSE WRITE (NOUT,*) WRITE (NOUT,*) 'NMAX too small' END IF * 99999 FORMAT (1X,A,I4) 99998 FORMAT (1X,A,I2,2A,/1X,2(A,I2,A,1P,E11.4),A,I2,A,1P,E11.4) 99997 FORMAT (1X,A,I2,A,1P,E11.4) 99996 FORMAT (1X,A,I2,A,'(',1P,E11.4,',',1P,E11.4,')') 99995 FORMAT (1X,A,1P,E8.1) 99994 FORMAT (1X,A,I2,A) 99993 FORMAT (1X,1P,E11.4) 99992 FORMAT (1X,'(',1P,E11.4,',',1P,E11.4,')') 99991 FORMAT (1X,A,I5,/1X,A,I5) END