* F08YXF Example Program Text * Mark 20 Release. NAG Copyright 2001. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX, LDA, LDB, LDVL, LDVR, LWORK PARAMETER (NMAX=10,LDA=NMAX,LDB=NMAX,LDVL=NMAX,LDVR=NMAX, + LWORK=6*NMAX) COMPLEX *16 CONE, CZERO PARAMETER (CONE=1.0D0,CZERO=0.0D0) * .. Local Scalars .. COMPLEX *16 E INTEGER I, ICOLS, IFAIL, IHI, ILO, INFO, IROWS, J, JWORK, + M, N LOGICAL ILEFT, IRIGHT CHARACTER COMPQ, COMPZ, HOWMNY, JOB, SIDE * .. Local Arrays .. COMPLEX *16 A(LDA,NMAX), ALPHA(NMAX), B(LDB,NMAX), + BETA(NMAX), TAU(NMAX), VL(LDVL,LDVL), + VR(LDVR,LDVR), WORK(LWORK) DOUBLE PRECISION LSCALE(NMAX), RSCALE(NMAX), RWORK(6*NMAX) LOGICAL SELECT(NMAX) CHARACTER CLABS(1), RLABS(1) * .. External Subroutines .. EXTERNAL F06TFF, F06THF, X04DBF, ZGEQRF, ZGGBAK, ZGGBAL, + ZGGHRD, ZHGEQZ, ZTGEVC, ZUNGQR, ZUNMQR * .. Intrinsic Functions .. INTRINSIC AIMAG, DBLE, NINT * .. Executable Statements .. WRITE (NOUT,*) 'F08YXF Example Program Results' * * ILEFT is TRUE if left eigenvectors are required * IRIGHT is TRUE if right eigenvectors are required * ILEFT = .TRUE. IRIGHT = .TRUE. * * Skip heading in data file * READ (NIN,*) READ (NIN,*) N IF (N.LE.NMAX) THEN * * READ matrix A from data file * READ (NIN,*) ((A(I,J),J=1,N),I=1,N) * * READ matrix B from data file * READ (NIN,*) ((B(I,J),J=1,N),I=1,N) * * Balance matrix pair (A,B) * JOB = 'B' CALL ZGGBAL(JOB,N,A,LDA,B,LDB,ILO,IHI,LSCALE,RSCALE,RWORK,INFO) * * Matrix A after balancing * IFAIL = 0 CALL X04DBF('General',' ',N,N,A,LDA,'Bracketed','F7.4', + 'Matrix A after balancing','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) WRITE (NOUT,*) * * Matrix B after balancing * IFAIL = 0 CALL X04DBF('General',' ',N,N,B,LDB,'Bracketed','F7.4', + 'Matrix B after balancing','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) WRITE (NOUT,*) * * Reduce B to triangular form using QR * IROWS = IHI + 1 - ILO ICOLS = N + 1 - ILO CALL ZGEQRF(IROWS,ICOLS,B(ILO,ILO),LDB,TAU,WORK,LWORK,INFO) * * Apply the orthogonal transformation to A * CALL ZUNMQR('L','C',IROWS,ICOLS,IROWS,B(ILO,ILO),LDB,TAU, + A(ILO,ILO),LDA,WORK,LWORK,INFO) * * Initialize VL (for left eigenvectors) * IF (ILEFT) THEN * CALL F06THF('General',N,N,CZERO,CONE,VL,LDVL) CALL F06TFF('Lower',IROWS-1,IROWS-1,B(ILO+1,ILO),LDB, + VL(ILO+1,ILO),LDVL) CALL ZUNGQR(IROWS,IROWS,IROWS,VL(ILO,ILO),LDVL,TAU,WORK, + LWORK,INFO) END IF * * Initialize VR for right eigenvectors * IF (IRIGHT) CALL F06THF('General',N,N,CZERO,CONE,VR,LDVR) * * Compute the generalized Hessenberg form of (A,B) * COMPQ = 'V' COMPZ = 'V' CALL ZGGHRD(COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,VL,LDVL,VR,LDVR, + INFO) * * Matrix A in generalized Hessenberg form * IFAIL = 0 CALL X04DBF('General',' ',N,N,A,LDA,'Bracketed','F7.3', + 'Matrix A in Hessenberg form','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) WRITE (NOUT,*) * * Matrix B in generalized Hessenberg form * IFAIL = 0 CALL X04DBF('General',' ',N,N,B,LDB,'Bracketed','F7.3', + 'Matrix B in Hessenberg form','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) * * Routine ZHGEQZ * Workspace query: JWORK = -1 * JWORK = -1 JOB = 'S' CALL ZHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,ALPHA,BETA, + VL,LDVL,VR,LDVR,WORK,JWORK,RWORK,INFO) WRITE (NOUT,*) WRITE (NOUT,99999) NINT(DBLE(WORK(1))) WRITE (NOUT,99998) LWORK WRITE (NOUT,*) * * Compute the generalized Schur form * if the workspace LWORK is adequate * IF (NINT(DBLE(WORK(1))).LE.LWORK) THEN CALL ZHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,ALPHA, + BETA,VL,LDVL,VR,LDVR,WORK,LWORK,RWORK,INFO) * * Print the generalized eigenvalues * Note: the actual values of beta are real and non-negative * WRITE (NOUT,99997) DO 20 I = 1, N IF (DBLE(BETA(I)).NE.0.0D0) THEN E = ALPHA(I)/BETA(I) WRITE (NOUT,99995) I, '(', DBLE(E), ',', AIMAG(E), ')' ELSE WRITE (NOUT,99996) I END IF 20 CONTINUE WRITE (NOUT,*) * * Compute left and right generalized eigenvectors * of the balanced matrix * HOWMNY = 'B' IF (ILEFT .AND. IRIGHT) THEN SIDE = 'B' ELSE IF (ILEFT) THEN SIDE = 'L' ELSE IF (IRIGHT) THEN SIDE = 'R' END IF * CALL ZTGEVC(SIDE,HOWMNY,SELECT,N,A,LDA,B,LDB,VL,LDVL,VR, + LDVR,N,M,WORK,RWORK,INFO) * * Compute right eigenvectors of the original matrix * IF (IRIGHT) THEN JOB = 'B' SIDE = 'R' * CALL ZGGBAK(JOB,SIDE,N,ILO,IHI,LSCALE,RSCALE,N,VR,LDVR, + INFO) * * Print the right eigenvectors * IFAIL = 0 CALL X04DBF('General',' ',N,N,VR,LDVR,'Bracketed','F7.4', + 'Right eigenvectors','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) WRITE (NOUT,*) END IF * * Compute left eigenvectors of the original matrix * IF (IRIGHT) THEN JOB = 'B' SIDE = 'L' * CALL ZGGBAK(JOB,SIDE,N,ILO,IHI,LSCALE,RSCALE,N,VL,LDVL, + INFO) * * Print the left eigenvectors * IFAIL = 0 CALL X04DBF('General',' ',N,N,VL,LDVL,'Bracketed','F7.4', + 'Left eigenvectors','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) END IF ELSE WRITE (NOUT,99994) END IF END IF * 99999 FORMAT (1X,'Minimal required LWORK = ',I6) 99998 FORMAT (1X,'Actual value of LWORK = ',I6) 99997 FORMAT (1X,'Generalized eigenvalues') 99996 FORMAT (1X,I4,' Infinite eigenvalue') 99995 FORMAT (1X,I4,5X,A,F7.3,A,F7.3,A) 99994 FORMAT (1X,'Insufficient workspace for array WORK',/' in F08XSF/', + 'ZHGEQZ') END