Mercurial > octave-dspies
changeset 8408:15c23c1c3c18
add missing blas & lapack sources
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Mon, 15 Dec 2008 13:49:16 +0100 |
| parents | 096c22ce2b0b |
| children | 7f6efc2b75d9 |
| files | libcruft/ChangeLog libcruft/blas/zsyrk.f libcruft/lapack/Makefile.in libcruft/lapack/cggbak.f libcruft/lapack/cggev.f libcruft/lapack/cgghrd.f libcruft/lapack/chgeqz.f libcruft/lapack/ctgevc.f libcruft/lapack/dggev.f libcruft/lapack/sggev.f libcruft/lapack/zggbak.f libcruft/lapack/zggev.f libcruft/lapack/zgghrd.f libcruft/lapack/zhgeqz.f libcruft/lapack/ztgevc.f |
| diffstat | 15 files changed, 6063 insertions(+), 85 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Mon Dec 15 13:49:11 2008 +0100 +++ b/libcruft/ChangeLog Mon Dec 15 13:49:16 2008 +0100 @@ -1,3 +1,12 @@ +2008-12-15 Jaroslav Hajek <highegg@gmail.com> + + * blas/zsyrk.f: New source. + * lapack/cggbak.f, lapack/cggev.f, lapack/cgghrd.f, lapack/chgeqz.f, + lapack/ctgevc.f, lapack/dggev.f, lapack/sggev.f, lapack/zggbak.f, + lapack/zggev.f, lapack/zgghrd.f, lapack/zhgeqz.f, lapack/ztgevc.f: + New sources. + * lapack/Makefile.in: Include them. + 2008-08-12 Thomas Treichl <Thomas.Treichl@gmx.net> * blas/icamax.f, blas/isamax.f: New files.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/zsyrk.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,294 @@ + SUBROUTINE ZSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC) +* .. Scalar Arguments .. + DOUBLE COMPLEX ALPHA,BETA + INTEGER K,LDA,LDC,N + CHARACTER TRANS,UPLO +* .. +* .. Array Arguments .. + DOUBLE COMPLEX A(LDA,*),C(LDC,*) +* .. +* +* Purpose +* ======= +* +* ZSYRK performs one of the symmetric rank k operations +* +* C := alpha*A*A' + beta*C, +* +* or +* +* C := alpha*A'*A + beta*C, +* +* where alpha and beta are scalars, C is an n by n symmetric matrix +* and A is an n by k matrix in the first case and a k by n matrix +* in the second case. +* +* Arguments +* ========== +* +* UPLO - CHARACTER*1. +* On entry, UPLO specifies whether the upper or lower +* triangular part of the array C is to be referenced as +* follows: +* +* UPLO = 'U' or 'u' Only the upper triangular part of C +* is to be referenced. +* +* UPLO = 'L' or 'l' Only the lower triangular part of C +* is to be referenced. +* +* Unchanged on exit. +* +* TRANS - CHARACTER*1. +* On entry, TRANS specifies the operation to be performed as +* follows: +* +* TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. +* +* TRANS = 'T' or 't' C := alpha*A'*A + beta*C. +* +* Unchanged on exit. +* +* N - INTEGER. +* On entry, N specifies the order of the matrix C. N must be +* at least zero. +* Unchanged on exit. +* +* K - INTEGER. +* On entry with TRANS = 'N' or 'n', K specifies the number +* of columns of the matrix A, and on entry with +* TRANS = 'T' or 't', K specifies the number of rows of the +* matrix A. K must be at least zero. +* Unchanged on exit. +* +* ALPHA - COMPLEX*16 . +* On entry, ALPHA specifies the scalar alpha. +* Unchanged on exit. +* +* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is +* k when TRANS = 'N' or 'n', and is n otherwise. +* Before entry with TRANS = 'N' or 'n', the leading n by k +* part of the array A must contain the matrix A, otherwise +* the leading k by n part of the array A must contain the +* matrix A. +* Unchanged on exit. +* +* LDA - INTEGER. +* On entry, LDA specifies the first dimension of A as declared +* in the calling (sub) program. When TRANS = 'N' or 'n' +* then LDA must be at least max( 1, n ), otherwise LDA must +* be at least max( 1, k ). +* Unchanged on exit. +* +* BETA - COMPLEX*16 . +* On entry, BETA specifies the scalar beta. +* Unchanged on exit. +* +* C - COMPLEX*16 array of DIMENSION ( LDC, n ). +* Before entry with UPLO = 'U' or 'u', the leading n by n +* upper triangular part of the array C must contain the upper +* triangular part of the symmetric matrix and the strictly +* lower triangular part of C is not referenced. On exit, the +* upper triangular part of the array C is overwritten by the +* upper triangular part of the updated matrix. +* Before entry with UPLO = 'L' or 'l', the leading n by n +* lower triangular part of the array C must contain the lower +* triangular part of the symmetric matrix and the strictly +* upper triangular part of C is not referenced. On exit, the +* lower triangular part of the array C is overwritten by the +* lower triangular part of the updated matrix. +* +* LDC - INTEGER. +* On entry, LDC specifies the first dimension of C as declared +* in the calling (sub) program. LDC must be at least +* max( 1, n ). +* Unchanged on exit. +* +* +* Level 3 Blas routine. +* +* -- Written on 8-February-1989. +* Jack Dongarra, Argonne National Laboratory. +* Iain Duff, AERE Harwell. +* Jeremy Du Croz, Numerical Algorithms Group Ltd. +* Sven Hammarling, Numerical Algorithms Group Ltd. +* +* +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Local Scalars .. + DOUBLE COMPLEX TEMP + INTEGER I,INFO,J,L,NROWA + LOGICAL UPPER +* .. +* .. Parameters .. + DOUBLE COMPLEX ONE + PARAMETER (ONE= (1.0D+0,0.0D+0)) + DOUBLE COMPLEX ZERO + PARAMETER (ZERO= (0.0D+0,0.0D+0)) +* .. +* +* Test the input parameters. +* + IF (LSAME(TRANS,'N')) THEN + NROWA = N + ELSE + NROWA = K + END IF + UPPER = LSAME(UPLO,'U') +* + INFO = 0 + IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN + INFO = 1 + ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. + + (.NOT.LSAME(TRANS,'T'))) THEN + INFO = 2 + ELSE IF (N.LT.0) THEN + INFO = 3 + ELSE IF (K.LT.0) THEN + INFO = 4 + ELSE IF (LDA.LT.MAX(1,NROWA)) THEN + INFO = 7 + ELSE IF (LDC.LT.MAX(1,N)) THEN + INFO = 10 + END IF + IF (INFO.NE.0) THEN + CALL XERBLA('ZSYRK ',INFO) + RETURN + END IF +* +* Quick return if possible. +* + IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. + + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN +* +* And when alpha.eq.zero. +* + IF (ALPHA.EQ.ZERO) THEN + IF (UPPER) THEN + IF (BETA.EQ.ZERO) THEN + DO 20 J = 1,N + DO 10 I = 1,J + C(I,J) = ZERO + 10 CONTINUE + 20 CONTINUE + ELSE + DO 40 J = 1,N + DO 30 I = 1,J + C(I,J) = BETA*C(I,J) + 30 CONTINUE + 40 CONTINUE + END IF + ELSE + IF (BETA.EQ.ZERO) THEN + DO 60 J = 1,N + DO 50 I = J,N + C(I,J) = ZERO + 50 CONTINUE + 60 CONTINUE + ELSE + DO 80 J = 1,N + DO 70 I = J,N + C(I,J) = BETA*C(I,J) + 70 CONTINUE + 80 CONTINUE + END IF + END IF + RETURN + END IF +* +* Start the operations. +* + IF (LSAME(TRANS,'N')) THEN +* +* Form C := alpha*A*A' + beta*C. +* + IF (UPPER) THEN + DO 130 J = 1,N + IF (BETA.EQ.ZERO) THEN + DO 90 I = 1,J + C(I,J) = ZERO + 90 CONTINUE + ELSE IF (BETA.NE.ONE) THEN + DO 100 I = 1,J + C(I,J) = BETA*C(I,J) + 100 CONTINUE + END IF + DO 120 L = 1,K + IF (A(J,L).NE.ZERO) THEN + TEMP = ALPHA*A(J,L) + DO 110 I = 1,J + C(I,J) = C(I,J) + TEMP*A(I,L) + 110 CONTINUE + END IF + 120 CONTINUE + 130 CONTINUE + ELSE + DO 180 J = 1,N + IF (BETA.EQ.ZERO) THEN + DO 140 I = J,N + C(I,J) = ZERO + 140 CONTINUE + ELSE IF (BETA.NE.ONE) THEN + DO 150 I = J,N + C(I,J) = BETA*C(I,J) + 150 CONTINUE + END IF + DO 170 L = 1,K + IF (A(J,L).NE.ZERO) THEN + TEMP = ALPHA*A(J,L) + DO 160 I = J,N + C(I,J) = C(I,J) + TEMP*A(I,L) + 160 CONTINUE + END IF + 170 CONTINUE + 180 CONTINUE + END IF + ELSE +* +* Form C := alpha*A'*A + beta*C. +* + IF (UPPER) THEN + DO 210 J = 1,N + DO 200 I = 1,J + TEMP = ZERO + DO 190 L = 1,K + TEMP = TEMP + A(L,I)*A(L,J) + 190 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = ALPHA*TEMP + ELSE + C(I,J) = ALPHA*TEMP + BETA*C(I,J) + END IF + 200 CONTINUE + 210 CONTINUE + ELSE + DO 240 J = 1,N + DO 230 I = J,N + TEMP = ZERO + DO 220 L = 1,K + TEMP = TEMP + A(L,I)*A(L,J) + 220 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = ALPHA*TEMP + ELSE + C(I,J) = ALPHA*TEMP + BETA*C(I,J) + END IF + 230 CONTINUE + 240 CONTINUE + END IF + END IF +* + RETURN +* +* End of ZSYRK . +* + END
--- a/libcruft/lapack/Makefile.in Mon Dec 15 13:49:11 2008 +0100 +++ b/libcruft/lapack/Makefile.in Mon Dec 15 13:49:16 2008 +0100 @@ -26,91 +26,129 @@ EXTERNAL_DISTFILES = $(DISTFILES) -FSRC = cbdsqr.f csrscl.f cgbcon.f cgbtf2.f cgbtrf.f cgbtrs.f \ - cgebak.f cgebal.f cgebd2.f cgebrd.f cgecon.f cgeesx.f cgeev.f \ - cgehd2.f cgehrd.f cgelq2.f cgelqf.f cgelsd.f cgelss.f cgelsy.f \ - cgeqp3.f cgeqpf.f cgeqr2.f cgeqrf.f cgesv.f cgesvd.f cgetf2.f \ - cgetrf.f cgetri.f cgetrs.f cggbal.f cgtsv.f cgttrf.f cgttrs.f \ - cgtts2.f cheev.f chetd2.f chetrd.f chseqr.f clabrd.f clacgv.f \ - clacn2.f clacon.f clacpy.f cladiv.f clahqr.f clahr2.f clahrd.f \ - claic1.f clals0.f clalsa.f clalsd.f clange.f clanhe.f clanhs.f \ - clantr.f claqp2.f claqps.f claqr0.f claqr1.f claqr2.f claqr3.f \ - claqr4.f claqr5.f clarf.f clarfb.f clarfg.f clarft.f clarfx.f \ - clartg.f clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f \ - classq.f claswp.f clatbs.f clatrd.f clatrs.f clatrz.f clauu2.f \ - clauum.f cpbcon.f cpbtf2.f cpbtrf.f cpbtrs.f cpocon.f cpotf2.f \ - cpotrf.f cpotri.f cpotrs.f cptsv.f cpttrf.f cpttrs.f cptts2.f crot.f \ - csteqr.f ctrcon.f ctrevc.f ctrexc.f ctrsen.f ctrsyl.f ctrti2.f \ - ctrtri.f ctrtrs.f ctzrzf.f cung2l.f cung2r.f cungbr.f cunghr.f \ - cungl2.f cunglq.f cungql.f cungqr.f cungtr.f cunm2r.f cunmbr.f \ - cunml2.f cunmlq.f cunmqr.f cunmr3.f cunmrz.f \ - dbdsqr.f dgbcon.f dgbtf2.f dgbtrf.f dgbtrs.f dgebak.f dgebal.f \ - dgebd2.f dgebrd.f dgecon.f dgeesx.f dgeev.f dgehd2.f dgehrd.f \ - dgelq2.f dgelqf.f dgelsd.f dgelss.f dgelsy.f dgeqp3.f dgeqpf.f \ - dgeqr2.f dgeqrf.f dgesv.f dgesvd.f dgetf2.f dgetrf.f dgetri.f \ - dgetrs.f dggbak.f dggbal.f dgghrd.f dgtsv.f dgttrf.f dgttrs.f \ - dgtts2.f dhgeqz.f dhseqr.f dlabad.f dlabrd.f dlacn2.f dlacon.f \ - dlacpy.f dladiv.f dlae2.f dlaed6.f dlaev2.f dlaexc.f dlag2.f \ - dlahqr.f dlahr2.f dlahrd.f dlaic1.f dlaln2.f dlals0.f dlalsa.f \ - dlalsd.f dlamc1.f dlamc2.f dlamc3.f dlamc4.f dlamc5.f dlamch.f \ - dlamrg.f dlange.f dlanhs.f dlanst.f dlansy.f dlantr.f dlanv2.f \ - dlapy2.f dlapy3.f dlaqp2.f dlaqps.f dlaqr0.f dlaqr1.f dlaqr2.f \ - dlaqr3.f dlaqr4.f dlaqr5.f dlarf.f dlarfb.f dlarfg.f dlarft.f \ - dlarfx.f dlartg.f dlarz.f dlarzb.f dlarzt.f dlas2.f dlascl.f \ - dlasd0.f dlasd1.f dlasd2.f dlasd3.f dlasd4.f dlasd5.f dlasd6.f \ - dlasd7.f dlasd8.f dlasda.f dlasdq.f dlasdt.f dlaset.f dlasq1.f \ - dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f dlasr.f dlasrt.f \ - dlassq.f dlasv2.f dlaswp.f dlasy2.f dlatbs.f dlatrd.f dlatrs.f \ - dlatrz.f dlauu2.f dlauum.f dlazq3.f dlazq4.f dorg2l.f dorg2r.f \ - dorgbr.f dorghr.f dorgl2.f dorglq.f dorgql.f dorgqr.f dorgtr.f \ - dorm2r.f dormbr.f dorml2.f dormlq.f dormqr.f dormr3.f dormrz.f \ - dpbcon.f dpbtf2.f dpbtrf.f dpbtrs.f dpocon.f dpotf2.f dpotrf.f \ - dpotri.f dpotrs.f dptsv.f dpttrf.f dpttrs.f dptts2.f drscl.f \ - dsteqr.f dsterf.f dsyev.f dsytd2.f dsytrd.f dtgevc.f dtrcon.f \ - dtrevc.f dtrexc.f dtrsen.f dtrsyl.f dtrti2.f dtrtri.f dtrtrs.f \ - dtzrzf.f dzsum1.f icmax1.f ieeeck.f ilaenv.f iparmq.f izmax1.f \ - sbdsqr.f sgbcon.f sgbtf2.f sgbtrf.f sgbtrs.f sgebak.f sgebal.f \ - sgebd2.f sgebrd.f sgecon.f sgeesx.f sgeev.f sgehd2.f sgehrd.f \ - sgelq2.f sgelqf.f sgelsd.f sgelss.f sgelsy.f sgeqp3.f sgeqpf.f \ - sgeqr2.f sgeqrf.f sgesv.f sgesvd.f sgetf2.f sgetrf.f sgetri.f \ - sgetrs.f sggbak.f sggbal.f sgghrd.f sgtsv.f sgttrf.f sgttrs.f \ - sgtts2.f shgeqz.f shseqr.f slabad.f slabrd.f slacn2.f slacon.f \ - slacpy.f sladiv.f slae2.f slaed6.f slaev2.f slaexc.f slag2.f \ - slahqr.f slahr2.f slahrd.f slaic1.f slaln2.f slals0.f slalsa.f \ - slalsd.f slamc1.f slamc2.f slamc3.f slamc4.f slamc5.f slamch.f \ - slamrg.f slange.f slanhs.f slanst.f slansy.f slantr.f slanv2.f \ - slapy2.f slapy3.f slaqp2.f slaqps.f slaqr0.f slaqr1.f slaqr2.f \ - slaqr3.f slaqr4.f slaqr5.f slarf.f slarfb.f slarfg.f slarft.f \ - slarfx.f slartg.f slarz.f slarzb.f slarzt.f slas2.f slascl.f \ - slasd0.f slasd1.f slasd2.f slasd3.f slasd4.f slasd5.f slasd6.f \ - slasd7.f slasd8.f slasda.f slasdq.f slasdt.f slaset.f slasq1.f \ - slasq2.f slasq3.f slasq4.f slasq5.f slasq6.f slasr.f slasrt.f \ - slassq.f slasv2.f slaswp.f slasy2.f slatbs.f slatrd.f slatrs.f \ - slatrz.f slauu2.f slauum.f slazq3.f slazq4.f sorg2l.f sorg2r.f \ - sorgbr.f sorghr.f sorgl2.f sorglq.f sorgql.f sorgqr.f sorgtr.f \ - sorm2r.f sormbr.f sorml2.f sormlq.f sormqr.f sormr3.f sormrz.f \ - spbcon.f spbtf2.f spbtrf.f spbtrs.f spocon.f spotf2.f spotrf.f \ - spotri.f spotrs.f sptsv.f spttrf.f spttrs.f sptts2.f srscl.f \ - ssteqr.f ssterf.f ssyev.f ssytd2.f ssytrd.f stgevc.f strcon.f \ - strevc.f strexc.f strsen.f strsyl.f strti2.f strtri.f strtrs.f \ - stzrzf.f scsum1.f zbdsqr.f zdrscl.f zgbcon.f zgbtf2.f zgbtrf.f zgbtrs.f \ - zgebak.f zgebal.f zgebd2.f zgebrd.f zgecon.f zgeesx.f zgeev.f \ - zgehd2.f zgehrd.f zgelq2.f zgelqf.f zgelsd.f zgelss.f zgelsy.f \ - zgeqp3.f zgeqpf.f zgeqr2.f zgeqrf.f zgesv.f zgesvd.f zgetf2.f \ - zgetrf.f zgetri.f zgetrs.f zggbal.f zgtsv.f zgttrf.f zgttrs.f \ - zgtts2.f zheev.f zhetd2.f zhetrd.f zhseqr.f zlabrd.f zlacgv.f \ - zlacn2.f zlacon.f zlacpy.f zladiv.f zlahqr.f zlahr2.f zlahrd.f \ - zlaic1.f zlals0.f zlalsa.f zlalsd.f zlange.f zlanhe.f zlanhs.f \ - zlantr.f zlaqp2.f zlaqps.f zlaqr0.f zlaqr1.f zlaqr2.f zlaqr3.f \ - zlaqr4.f zlaqr5.f zlarf.f zlarfb.f zlarfg.f zlarft.f zlarfx.f \ - zlartg.f zlarz.f zlarzb.f zlarzt.f zlascl.f zlaset.f zlasr.f \ - zlassq.f zlaswp.f zlatbs.f zlatrd.f zlatrs.f zlatrz.f zlauu2.f \ - zlauum.f zpbcon.f zpbtf2.f zpbtrf.f zpbtrs.f zpocon.f zpotf2.f \ - zpotrf.f zpotri.f zpotrs.f zptsv.f zpttrf.f zpttrs.f zptts2.f zrot.f \ - zsteqr.f ztrcon.f ztrevc.f ztrexc.f ztrsen.f ztrsyl.f ztrti2.f \ - ztrtri.f ztrtrs.f ztzrzf.f zung2l.f zung2r.f zungbr.f zunghr.f \ - zungl2.f zunglq.f zungql.f zungqr.f zungtr.f zunm2r.f zunmbr.f \ - zunml2.f zunmlq.f zunmqr.f zunmr3.f zunmrz.f +FSRC = + cbdsqr.f cgbcon.f cgbtf2.f cgbtrf.f cgbtrs.f cgbtrs.f \ + cgebak.f cgebal.f cgebd2.f cgebrd.f cgecon.f cgecon.f \ + cgeesx.f cgeev.f cgehd2.f cgehrd.f cgelq2.f cgelq2.f \ + cgelqf.f cgelsd.f cgelss.f cgelsy.f cgeqp3.f cgeqp3.f \ + cgeqpf.f cgeqr2.f cgeqrf.f cgesvd.f cgesv.f cgesv.f \ + cgetf2.f cgetrf.f cgetri.f cgetrs.f cggbak.f cggbak.f \ + cggbal.f cggev.f cgghrd.f cgtsv.f cgttrf.f cgttrf.f \ + cgttrs.f cgtts2.f cheev.f chetd2.f chetrd.f chetrd.f \ + chgeqz.f chseqr.f clabrd.f clacgv.f clacn2.f clacn2.f \ + clacon.f clacpy.f cladiv.f clahqr.f clahr2.f clahr2.f \ + clahrd.f claic1.f clals0.f clalsa.f clalsd.f clalsd.f \ + clange.f clanhe.f clanhs.f clantr.f claqp2.f claqp2.f \ + claqps.f claqr0.f claqr1.f claqr2.f claqr3.f claqr3.f \ + claqr4.f claqr5.f clarfb.f clarf.f clarfg.f clarfg.f \ + clarft.f clarfx.f clartg.f clarzb.f clarz.f clarz.f \ + clarzt.f clascl.f claset.f clasr.f classq.f classq.f \ + claswp.f clatbs.f clatrd.f clatrs.f clatrz.f clatrz.f \ + clauu2.f clauum.f cpbcon.f cpbtf2.f cpbtrf.f cpbtrf.f \ + cpbtrs.f cpocon.f cpotf2.f cpotrf.f cpotri.f cpotri.f \ + cpotrs.f cptsv.f cpttrf.f cpttrs.f cptts2.f cptts2.f \ + crot.f csrscl.f csteqr.f ctgevc.f ctrcon.f ctrcon.f \ + ctrevc.f ctrexc.f ctrsen.f ctrsyl.f ctrti2.f ctrti2.f \ + ctrtri.f ctrtrs.f ctzrzf.f cung2l.f cung2r.f cung2r.f \ + cungbr.f cunghr.f cungl2.f cunglq.f cungql.f cungql.f \ + cungqr.f cungtr.f cunm2r.f cunmbr.f cunml2.f cunml2.f \ + cunmlq.f cunmqr.f cunmr3.f cunmrz.f dbdsqr.f dbdsqr.f \ + dgbcon.f dgbtf2.f dgbtrf.f dgbtrs.f dgebak.f dgebak.f \ + dgebal.f dgebd2.f dgebrd.f dgecon.f dgeesx.f dgeesx.f \ + dgeev.f dgehd2.f dgehrd.f dgelq2.f dgelqf.f dgelqf.f \ + dgelsd.f dgelss.f dgelsy.f dgeqp3.f dgeqpf.f dgeqpf.f \ + dgeqr2.f dgeqrf.f dgesvd.f dgesv.f dgetf2.f dgetf2.f \ + dgetrf.f dgetri.f dgetrs.f dggbak.f dggbal.f dggbal.f \ + dggev.f dgghrd.f dgtsv.f dgttrf.f dgttrs.f dgttrs.f \ + dgtts2.f dhgeqz.f dhseqr.f dlabad.f dlabrd.f dlabrd.f \ + dlacn2.f dlacon.f dlacpy.f dladiv.f dlae2.f dlae2.f \ + dlaed6.f dlaev2.f dlaexc.f dlag2.f dlahqr.f dlahqr.f \ + dlahr2.f dlahrd.f dlaic1.f dlaln2.f dlals0.f dlals0.f \ + dlalsa.f dlalsd.f dlamc1.f dlamc2.f dlamc3.f dlamc3.f \ + dlamc4.f dlamc5.f dlamch.f dlamrg.f dlange.f dlange.f \ + dlanhs.f dlanst.f dlansy.f dlantr.f dlanv2.f dlanv2.f \ + dlapy2.f dlapy3.f dlaqp2.f dlaqps.f dlaqr0.f dlaqr0.f \ + dlaqr1.f dlaqr2.f dlaqr3.f dlaqr4.f dlaqr5.f dlaqr5.f \ + dlarfb.f dlarf.f dlarfg.f dlarft.f dlarfx.f dlarfx.f \ + dlartg.f dlarzb.f dlarz.f dlarzt.f dlas2.f dlas2.f \ + dlascl.f dlasd0.f dlasd1.f dlasd2.f dlasd3.f dlasd3.f \ + dlasd4.f dlasd5.f dlasd6.f dlasd7.f dlasd8.f dlasd8.f \ + dlasda.f dlasdq.f dlasdt.f dlaset.f dlasq1.f dlasq1.f \ + dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f dlasq6.f \ + dlasr.f dlasrt.f dlassq.f dlasv2.f dlaswp.f dlaswp.f \ + dlasy2.f dlatbs.f dlatrd.f dlatrs.f dlatrz.f dlatrz.f \ + dlauu2.f dlauum.f dlazq3.f dlazq4.f dorg2l.f dorg2l.f \ + dorg2r.f dorgbr.f dorghr.f dorgl2.f dorglq.f dorglq.f \ + dorgql.f dorgqr.f dorgtr.f dorm2r.f dormbr.f dormbr.f \ + dorml2.f dormlq.f dormqr.f dormr3.f dormrz.f dormrz.f \ + dpbcon.f dpbtf2.f dpbtrf.f dpbtrs.f dpocon.f dpocon.f \ + dpotf2.f dpotrf.f dpotri.f dpotrs.f dptsv.f dptsv.f \ + dpttrf.f dpttrs.f dptts2.f drscl.f dsteqr.f dsteqr.f \ + dsterf.f dsyev.f dsytd2.f dsytrd.f dtgevc.f dtgevc.f \ + dtrcon.f dtrevc.f dtrexc.f dtrsen.f dtrsyl.f dtrsyl.f \ + dtrti2.f dtrtri.f dtrtrs.f dtzrzf.f dzsum1.f dzsum1.f \ + icmax1.f ieeeck.f ilaenv.f iparmq.f izmax1.f izmax1.f \ + sbdsqr.f scsum1.f sgbcon.f sgbtf2.f sgbtrf.f sgbtrf.f \ + sgbtrs.f sgebak.f sgebal.f sgebd2.f sgebrd.f sgebrd.f \ + sgecon.f sgeesx.f sgeev.f sgehd2.f sgehrd.f sgehrd.f \ + sgelq2.f sgelqf.f sgelsd.f sgelss.f sgelsy.f sgelsy.f \ + sgeqp3.f sgeqpf.f sgeqr2.f sgeqrf.f sgesvd.f sgesvd.f \ + sgesv.f sgetf2.f sgetrf.f sgetri.f sgetrs.f sgetrs.f \ + sggbak.f sggbal.f sggev.f sgghrd.f sgtsv.f sgtsv.f \ + sgttrf.f sgttrs.f sgtts2.f shgeqz.f shseqr.f shseqr.f \ + slabad.f slabrd.f slacn2.f slacon.f slacpy.f slacpy.f \ + sladiv.f slae2.f slaed6.f slaev2.f slaexc.f slaexc.f \ + slag2.f slahqr.f slahr2.f slahrd.f slaic1.f slaic1.f \ + slaln2.f slals0.f slalsa.f slalsd.f slamc1.f slamc1.f \ + slamc2.f slamc3.f slamc4.f slamc5.f slamch.f slamch.f \ + slamrg.f slange.f slanhs.f slanst.f slansy.f slansy.f \ + slantr.f slanv2.f slapy2.f slapy3.f slaqp2.f slaqp2.f \ + slaqps.f slaqr0.f slaqr1.f slaqr2.f slaqr3.f slaqr3.f \ + slaqr4.f slaqr5.f slarfb.f slarf.f slarfg.f slarfg.f \ + slarft.f slarfx.f slartg.f slarzb.f slarz.f slarz.f \ + slarzt.f slas2.f slascl.f slasd0.f slasd1.f slasd1.f \ + slasd2.f slasd3.f slasd4.f slasd5.f slasd6.f slasd6.f \ + slasd7.f slasd8.f slasda.f slasdq.f slasdt.f slasdt.f \ + slaset.f slasq1.f slasq2.f slasq3.f slasq4.f slasq4.f \ + slasq5.f slasq6.f slasr.f slasrt.f slassq.f slassq.f \ + slasv2.f slaswp.f slasy2.f slatbs.f slatrd.f slatrd.f \ + slatrs.f slatrz.f slauu2.f slauum.f slazq3.f slazq3.f \ + slazq4.f sorg2l.f sorg2r.f sorgbr.f sorghr.f sorghr.f \ + sorgl2.f sorglq.f sorgql.f sorgqr.f sorgtr.f sorgtr.f \ + sorm2r.f sormbr.f sorml2.f sormlq.f sormqr.f sormqr.f \ + sormr3.f sormrz.f spbcon.f spbtf2.f spbtrf.f spbtrf.f \ + spbtrs.f spocon.f spotf2.f spotrf.f spotri.f spotri.f \ + spotrs.f sptsv.f spttrf.f spttrs.f sptts2.f sptts2.f \ + srscl.f ssteqr.f ssterf.f ssyev.f ssytd2.f ssytd2.f \ + ssytrd.f stgevc.f strcon.f strevc.f strexc.f strexc.f \ + strsen.f strsyl.f strti2.f strtri.f strtrs.f strtrs.f \ + stzrzf.f zbdsqr.f zdrscl.f zgbcon.f zgbtf2.f zgbtf2.f \ + zgbtrf.f zgbtrs.f zgebak.f zgebal.f zgebd2.f zgebd2.f \ + zgebrd.f zgecon.f zgeesx.f zgeev.f zgehd2.f zgehd2.f \ + zgehrd.f zgelq2.f zgelqf.f zgelsd.f zgelss.f zgelss.f \ + zgelsy.f zgeqp3.f zgeqpf.f zgeqr2.f zgeqrf.f zgeqrf.f \ + zgesvd.f zgesv.f zgetf2.f zgetrf.f zgetri.f zgetri.f \ + zgetrs.f zggbak.f zggbal.f zggev.f zgghrd.f zgghrd.f \ + zgtsv.f zgttrf.f zgttrs.f zgtts2.f zheev.f zheev.f \ + zhetd2.f zhetrd.f zhgeqz.f zhseqr.f zlabrd.f zlabrd.f \ + zlacgv.f zlacn2.f zlacon.f zlacpy.f zladiv.f zladiv.f \ + zlahqr.f zlahr2.f zlahrd.f zlaic1.f zlals0.f zlals0.f \ + zlalsa.f zlalsd.f zlange.f zlanhe.f zlanhs.f zlanhs.f \ + zlantr.f zlaqp2.f zlaqps.f zlaqr0.f zlaqr1.f zlaqr1.f \ + zlaqr2.f zlaqr3.f zlaqr4.f zlaqr5.f zlarfb.f zlarfb.f \ + zlarf.f zlarfg.f zlarft.f zlarfx.f zlartg.f zlartg.f \ + zlarzb.f zlarz.f zlarzt.f zlascl.f zlaset.f zlaset.f \ + zlasr.f zlassq.f zlaswp.f zlatbs.f zlatrd.f zlatrd.f \ + zlatrs.f zlatrz.f zlauu2.f zlauum.f zpbcon.f zpbcon.f \ + zpbtf2.f zpbtrf.f zpbtrs.f zpocon.f zpotf2.f zpotf2.f \ + zpotrf.f zpotri.f zpotrs.f zptsv.f zpttrf.f zpttrf.f \ + zpttrs.f zptts2.f zrot.f zsteqr.f ztgevc.f ztgevc.f \ + ztrcon.f ztrevc.f ztrexc.f ztrsen.f ztrsyl.f ztrsyl.f \ + ztrti2.f ztrtri.f ztrtrs.f ztzrzf.f zung2l.f zung2l.f \ + zung2r.f zungbr.f zunghr.f zungl2.f zunglq.f zunglq.f \ + zungql.f zungqr.f zungtr.f zunm2r.f zunmbr.f zunmbr.f \ + zunml2.f zunmlq.f zunmqr.f zunmr3.f zunmrz.f zunmrz.f + include $(TOPDIR)/Makeconf
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/cggbak.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,220 @@ + SUBROUTINE CGGBAK( JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, + $ LDV, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOB, SIDE + INTEGER IHI, ILO, INFO, LDV, M, N +* .. +* .. Array Arguments .. + REAL LSCALE( * ), RSCALE( * ) + COMPLEX V( LDV, * ) +* .. +* +* Purpose +* ======= +* +* CGGBAK forms the right or left eigenvectors of a complex generalized +* eigenvalue problem A*x = lambda*B*x, by backward transformation on +* the computed eigenvectors of the balanced pair of matrices output by +* CGGBAL. +* +* Arguments +* ========= +* +* JOB (input) CHARACTER*1 +* Specifies the type of backward transformation required: +* = 'N': do nothing, return immediately; +* = 'P': do backward transformation for permutation only; +* = 'S': do backward transformation for scaling only; +* = 'B': do backward transformations for both permutation and +* scaling. +* JOB must be the same as the argument JOB supplied to CGGBAL. +* +* SIDE (input) CHARACTER*1 +* = 'R': V contains right eigenvectors; +* = 'L': V contains left eigenvectors. +* +* N (input) INTEGER +* The number of rows of the matrix V. N >= 0. +* +* ILO (input) INTEGER +* IHI (input) INTEGER +* The integers ILO and IHI determined by CGGBAL. +* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. +* +* LSCALE (input) REAL array, dimension (N) +* Details of the permutations and/or scaling factors applied +* to the left side of A and B, as returned by CGGBAL. +* +* RSCALE (input) REAL array, dimension (N) +* Details of the permutations and/or scaling factors applied +* to the right side of A and B, as returned by CGGBAL. +* +* M (input) INTEGER +* The number of columns of the matrix V. M >= 0. +* +* V (input/output) COMPLEX array, dimension (LDV,M) +* On entry, the matrix of right or left eigenvectors to be +* transformed, as returned by CTGEVC. +* On exit, V is overwritten by the transformed eigenvectors. +* +* LDV (input) INTEGER +* The leading dimension of the matrix V. LDV >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* Further Details +* =============== +* +* See R.C. Ward, Balancing the generalized eigenvalue problem, +* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LEFTV, RIGHTV + INTEGER I, K +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CSSCAL, CSWAP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + RIGHTV = LSAME( SIDE, 'R' ) + LEFTV = LSAME( SIDE, 'L' ) +* + INFO = 0 + IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. + $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN + INFO = -1 + ELSE IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( ILO.LT.1 ) THEN + INFO = -4 + ELSE IF( N.EQ.0 .AND. IHI.EQ.0 .AND. ILO.NE.1 ) THEN + INFO = -4 + ELSE IF( N.GT.0 .AND. ( IHI.LT.ILO .OR. IHI.GT.MAX( 1, N ) ) ) + $ THEN + INFO = -5 + ELSE IF( N.EQ.0 .AND. ILO.EQ.1 .AND. IHI.NE.0 ) THEN + INFO = -5 + ELSE IF( M.LT.0 ) THEN + INFO = -8 + ELSE IF( LDV.LT.MAX( 1, N ) ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGGBAK', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN + IF( M.EQ.0 ) + $ RETURN + IF( LSAME( JOB, 'N' ) ) + $ RETURN +* + IF( ILO.EQ.IHI ) + $ GO TO 30 +* +* Backward balance +* + IF( LSAME( JOB, 'S' ) .OR. LSAME( JOB, 'B' ) ) THEN +* +* Backward transformation on right eigenvectors +* + IF( RIGHTV ) THEN + DO 10 I = ILO, IHI + CALL CSSCAL( M, RSCALE( I ), V( I, 1 ), LDV ) + 10 CONTINUE + END IF +* +* Backward transformation on left eigenvectors +* + IF( LEFTV ) THEN + DO 20 I = ILO, IHI + CALL CSSCAL( M, LSCALE( I ), V( I, 1 ), LDV ) + 20 CONTINUE + END IF + END IF +* +* Backward permutation +* + 30 CONTINUE + IF( LSAME( JOB, 'P' ) .OR. LSAME( JOB, 'B' ) ) THEN +* +* Backward permutation on right eigenvectors +* + IF( RIGHTV ) THEN + IF( ILO.EQ.1 ) + $ GO TO 50 + DO 40 I = ILO - 1, 1, -1 + K = RSCALE( I ) + IF( K.EQ.I ) + $ GO TO 40 + CALL CSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 40 CONTINUE +* + 50 CONTINUE + IF( IHI.EQ.N ) + $ GO TO 70 + DO 60 I = IHI + 1, N + K = RSCALE( I ) + IF( K.EQ.I ) + $ GO TO 60 + CALL CSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 60 CONTINUE + END IF +* +* Backward permutation on left eigenvectors +* + 70 CONTINUE + IF( LEFTV ) THEN + IF( ILO.EQ.1 ) + $ GO TO 90 + DO 80 I = ILO - 1, 1, -1 + K = LSCALE( I ) + IF( K.EQ.I ) + $ GO TO 80 + CALL CSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 80 CONTINUE +* + 90 CONTINUE + IF( IHI.EQ.N ) + $ GO TO 110 + DO 100 I = IHI + 1, N + K = LSCALE( I ) + IF( K.EQ.I ) + $ GO TO 100 + CALL CSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 100 CONTINUE + END IF + END IF +* + 110 CONTINUE +* + RETURN +* +* End of CGGBAK +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/cggev.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,454 @@ + SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, + $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO ) +* +* -- LAPACK driver routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOBVL, JOBVR + INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N +* .. +* .. Array Arguments .. + REAL RWORK( * ) + COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), + $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), + $ WORK( * ) +* .. +* +* Purpose +* ======= +* +* CGGEV computes for a pair of N-by-N complex nonsymmetric matrices +* (A,B), the generalized eigenvalues, and optionally, the left and/or +* right generalized eigenvectors. +* +* A generalized eigenvalue for a pair of matrices (A,B) is a scalar +* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is +* singular. It is usually represented as the pair (alpha,beta), as +* there is a reasonable interpretation for beta=0, and even for both +* being zero. +* +* The right generalized eigenvector v(j) corresponding to the +* generalized eigenvalue lambda(j) of (A,B) satisfies +* +* A * v(j) = lambda(j) * B * v(j). +* +* The left generalized eigenvector u(j) corresponding to the +* generalized eigenvalues lambda(j) of (A,B) satisfies +* +* u(j)**H * A = lambda(j) * u(j)**H * B +* +* where u(j)**H is the conjugate-transpose of u(j). +* +* Arguments +* ========= +* +* JOBVL (input) CHARACTER*1 +* = 'N': do not compute the left generalized eigenvectors; +* = 'V': compute the left generalized eigenvectors. +* +* JOBVR (input) CHARACTER*1 +* = 'N': do not compute the right generalized eigenvectors; +* = 'V': compute the right generalized eigenvectors. +* +* N (input) INTEGER +* The order of the matrices A, B, VL, and VR. N >= 0. +* +* A (input/output) COMPLEX array, dimension (LDA, N) +* On entry, the matrix A in the pair (A,B). +* On exit, A has been overwritten. +* +* LDA (input) INTEGER +* The leading dimension of A. LDA >= max(1,N). +* +* B (input/output) COMPLEX array, dimension (LDB, N) +* On entry, the matrix B in the pair (A,B). +* On exit, B has been overwritten. +* +* LDB (input) INTEGER +* The leading dimension of B. LDB >= max(1,N). +* +* ALPHA (output) COMPLEX array, dimension (N) +* BETA (output) COMPLEX array, dimension (N) +* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the +* generalized eigenvalues. +* +* Note: the quotients ALPHA(j)/BETA(j) may easily over- or +* underflow, and BETA(j) may even be zero. Thus, the user +* should avoid naively computing the ratio alpha/beta. +* However, ALPHA will be always less than and usually +* comparable with norm(A) in magnitude, and BETA always less +* than and usually comparable with norm(B). +* +* VL (output) COMPLEX array, dimension (LDVL,N) +* If JOBVL = 'V', the left generalized eigenvectors u(j) are +* stored one after another in the columns of VL, in the same +* order as their eigenvalues. +* Each eigenvector is scaled so the largest component has +* abs(real part) + abs(imag. part) = 1. +* Not referenced if JOBVL = 'N'. +* +* LDVL (input) INTEGER +* The leading dimension of the matrix VL. LDVL >= 1, and +* if JOBVL = 'V', LDVL >= N. +* +* VR (output) COMPLEX array, dimension (LDVR,N) +* If JOBVR = 'V', the right generalized eigenvectors v(j) are +* stored one after another in the columns of VR, in the same +* order as their eigenvalues. +* Each eigenvector is scaled so the largest component has +* abs(real part) + abs(imag. part) = 1. +* Not referenced if JOBVR = 'N'. +* +* LDVR (input) INTEGER +* The leading dimension of the matrix VR. LDVR >= 1, and +* if JOBVR = 'V', LDVR >= N. +* +* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,2*N). +* For good performance, LWORK must generally be larger. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* RWORK (workspace/output) REAL array, dimension (8*N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value. +* =1,...,N: +* The QZ iteration failed. No eigenvectors have been +* calculated, but ALPHA(j) and BETA(j) should be +* correct for j=INFO+1,...,N. +* > N: =N+1: other then QZ iteration failed in SHGEQZ, +* =N+2: error return from STGEVC. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), + $ CONE = ( 1.0E0, 0.0E0 ) ) +* .. +* .. Local Scalars .. + LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY + CHARACTER CHTEMP + INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO, + $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR, + $ LWKMIN, LWKOPT + REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, + $ SMLNUM, TEMP + COMPLEX X +* .. +* .. Local Arrays .. + LOGICAL LDUMMA( 1 ) +* .. +* .. External Subroutines .. + EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY, + $ CLASCL, CLASET, CTGEVC, CUNGQR, CUNMQR, SLABAD, + $ XERBLA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + REAL CLANGE, SLAMCH + EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, MAX, REAL, SQRT +* .. +* .. Statement Functions .. + REAL ABS1 +* .. +* .. Statement Function definitions .. + ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) +* .. +* .. Executable Statements .. +* +* Decode the input arguments +* + IF( LSAME( JOBVL, 'N' ) ) THEN + IJOBVL = 1 + ILVL = .FALSE. + ELSE IF( LSAME( JOBVL, 'V' ) ) THEN + IJOBVL = 2 + ILVL = .TRUE. + ELSE + IJOBVL = -1 + ILVL = .FALSE. + END IF +* + IF( LSAME( JOBVR, 'N' ) ) THEN + IJOBVR = 1 + ILVR = .FALSE. + ELSE IF( LSAME( JOBVR, 'V' ) ) THEN + IJOBVR = 2 + ILVR = .TRUE. + ELSE + IJOBVR = -1 + ILVR = .FALSE. + END IF + ILV = ILVL .OR. ILVR +* +* Test the input arguments +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( IJOBVL.LE.0 ) THEN + INFO = -1 + ELSE IF( IJOBVR.LE.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN + INFO = -11 + ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN + INFO = -13 + END IF +* +* Compute workspace +* (Note: Comments in the code beginning "Workspace:" describe the +* minimal amount of workspace needed at that point in the code, +* as well as the preferred amount for good performance. +* NB refers to the optimal block size for the immediately +* following subroutine, as returned by ILAENV. The workspace is +* computed assuming ILO = 1 and IHI = N, the worst case.) +* + IF( INFO.EQ.0 ) THEN + LWKMIN = MAX( 1, 2*N ) + LWKOPT = MAX( 1, N + N*ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) ) + LWKOPT = MAX( LWKOPT, N + + $ N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, 0 ) ) + IF( ILVL ) THEN + LWKOPT = MAX( LWKOPT, N + + $ N*ILAENV( 1, 'CUNGQR', ' ', N, 1, N, -1 ) ) + END IF + WORK( 1 ) = LWKOPT +* + IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) + $ INFO = -15 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGGEV ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Get machine constants +* + EPS = SLAMCH( 'E' )*SLAMCH( 'B' ) + SMLNUM = SLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM + CALL SLABAD( SMLNUM, BIGNUM ) + SMLNUM = SQRT( SMLNUM ) / EPS + BIGNUM = ONE / SMLNUM +* +* Scale A if max element outside range [SMLNUM,BIGNUM] +* + ANRM = CLANGE( 'M', N, N, A, LDA, RWORK ) + ILASCL = .FALSE. + IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN + ANRMTO = SMLNUM + ILASCL = .TRUE. + ELSE IF( ANRM.GT.BIGNUM ) THEN + ANRMTO = BIGNUM + ILASCL = .TRUE. + END IF + IF( ILASCL ) + $ CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) +* +* Scale B if max element outside range [SMLNUM,BIGNUM] +* + BNRM = CLANGE( 'M', N, N, B, LDB, RWORK ) + ILBSCL = .FALSE. + IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN + BNRMTO = SMLNUM + ILBSCL = .TRUE. + ELSE IF( BNRM.GT.BIGNUM ) THEN + BNRMTO = BIGNUM + ILBSCL = .TRUE. + END IF + IF( ILBSCL ) + $ CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) +* +* Permute the matrices A, B to isolate eigenvalues if possible +* (Real Workspace: need 6*N) +* + ILEFT = 1 + IRIGHT = N + 1 + IRWRK = IRIGHT + N + CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR ) +* +* Reduce B to triangular form (QR decomposition of B) +* (Complex Workspace: need N, prefer N*NB) +* + IROWS = IHI + 1 - ILO + IF( ILV ) THEN + ICOLS = N + 1 - ILO + ELSE + ICOLS = IROWS + END IF + ITAU = 1 + IWRK = ITAU + IROWS + CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) +* +* Apply the orthogonal transformation to matrix A +* (Complex Workspace: need N, prefer N*NB) +* + CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, + $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), + $ LWORK+1-IWRK, IERR ) +* +* Initialize VL +* (Complex Workspace: need N, prefer N*NB) +* + IF( ILVL ) THEN + CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL ) + IF( IROWS.GT.1 ) THEN + CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, + $ VL( ILO+1, ILO ), LDVL ) + END IF + CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, + $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) + END IF +* +* Initialize VR +* + IF( ILVR ) + $ CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR ) +* +* Reduce to generalized Hessenberg form +* + IF( ILV ) THEN +* +* Eigenvectors requested -- work on whole matrix. +* + CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, + $ LDVL, VR, LDVR, IERR ) + ELSE + CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, + $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR ) + END IF +* +* Perform QZ algorithm (Compute eigenvalues, and optionally, the +* Schur form and Schur vectors) +* (Complex Workspace: need N) +* (Real Workspace: need N) +* + IWRK = ITAU + IF( ILV ) THEN + CHTEMP = 'S' + ELSE + CHTEMP = 'E' + END IF + CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, + $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ), + $ LWORK+1-IWRK, RWORK( IRWRK ), IERR ) + IF( IERR.NE.0 ) THEN + IF( IERR.GT.0 .AND. IERR.LE.N ) THEN + INFO = IERR + ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN + INFO = IERR - N + ELSE + INFO = N + 1 + END IF + GO TO 70 + END IF +* +* Compute Eigenvectors +* (Real Workspace: need 2*N) +* (Complex Workspace: need 2*N) +* + IF( ILV ) THEN + IF( ILVL ) THEN + IF( ILVR ) THEN + CHTEMP = 'B' + ELSE + CHTEMP = 'L' + END IF + ELSE + CHTEMP = 'R' + END IF +* + CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, + $ VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ), + $ IERR ) + IF( IERR.NE.0 ) THEN + INFO = N + 2 + GO TO 70 + END IF +* +* Undo balancing on VL and VR and normalization +* (Workspace: none needed) +* + IF( ILVL ) THEN + CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), N, VL, LDVL, IERR ) + DO 30 JC = 1, N + TEMP = ZERO + DO 10 JR = 1, N + TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) ) + 10 CONTINUE + IF( TEMP.LT.SMLNUM ) + $ GO TO 30 + TEMP = ONE / TEMP + DO 20 JR = 1, N + VL( JR, JC ) = VL( JR, JC )*TEMP + 20 CONTINUE + 30 CONTINUE + END IF + IF( ILVR ) THEN + CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), N, VR, LDVR, IERR ) + DO 60 JC = 1, N + TEMP = ZERO + DO 40 JR = 1, N + TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) ) + 40 CONTINUE + IF( TEMP.LT.SMLNUM ) + $ GO TO 60 + TEMP = ONE / TEMP + DO 50 JR = 1, N + VR( JR, JC ) = VR( JR, JC )*TEMP + 50 CONTINUE + 60 CONTINUE + END IF + END IF +* +* Undo scaling if necessary +* + IF( ILASCL ) + $ CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) +* + IF( ILBSCL ) + $ CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) +* + 70 CONTINUE + WORK( 1 ) = LWKOPT +* + RETURN +* +* End of CGGEV +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/cgghrd.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,264 @@ + SUBROUTINE CGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, + $ LDQ, Z, LDZ, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER COMPQ, COMPZ + INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ Z( LDZ, * ) +* .. +* +* Purpose +* ======= +* +* CGGHRD reduces a pair of complex matrices (A,B) to generalized upper +* Hessenberg form using unitary transformations, where A is a +* general matrix and B is upper triangular. The form of the generalized +* eigenvalue problem is +* A*x = lambda*B*x, +* and B is typically made upper triangular by computing its QR +* factorization and moving the unitary matrix Q to the left side +* of the equation. +* +* This subroutine simultaneously reduces A to a Hessenberg matrix H: +* Q**H*A*Z = H +* and transforms B to another upper triangular matrix T: +* Q**H*B*Z = T +* in order to reduce the problem to its standard form +* H*y = lambda*T*y +* where y = Z**H*x. +* +* The unitary matrices Q and Z are determined as products of Givens +* rotations. They may either be formed explicitly, or they may be +* postmultiplied into input matrices Q1 and Z1, so that +* Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H +* Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H +* If Q1 is the unitary matrix from the QR factorization of B in the +* original equation A*x = lambda*B*x, then CGGHRD reduces the original +* problem to generalized Hessenberg form. +* +* Arguments +* ========= +* +* COMPQ (input) CHARACTER*1 +* = 'N': do not compute Q; +* = 'I': Q is initialized to the unit matrix, and the +* unitary matrix Q is returned; +* = 'V': Q must contain a unitary matrix Q1 on entry, +* and the product Q1*Q is returned. +* +* COMPZ (input) CHARACTER*1 +* = 'N': do not compute Q; +* = 'I': Q is initialized to the unit matrix, and the +* unitary matrix Q is returned; +* = 'V': Q must contain a unitary matrix Q1 on entry, +* and the product Q1*Q is returned. +* +* N (input) INTEGER +* The order of the matrices A and B. N >= 0. +* +* ILO (input) INTEGER +* IHI (input) INTEGER +* ILO and IHI mark the rows and columns of A which are to be +* reduced. It is assumed that A is already upper triangular +* in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are +* normally set by a previous call to CGGBAL; otherwise they +* should be set to 1 and N respectively. +* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. +* +* A (input/output) COMPLEX array, dimension (LDA, N) +* On entry, the N-by-N general matrix to be reduced. +* On exit, the upper triangle and the first subdiagonal of A +* are overwritten with the upper Hessenberg matrix H, and the +* rest is set to zero. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* B (input/output) COMPLEX array, dimension (LDB, N) +* On entry, the N-by-N upper triangular matrix B. +* On exit, the upper triangular matrix T = Q**H B Z. The +* elements below the diagonal are set to zero. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,N). +* +* Q (input/output) COMPLEX array, dimension (LDQ, N) +* On entry, if COMPQ = 'V', the unitary matrix Q1, typically +* from the QR factorization of B. +* On exit, if COMPQ='I', the unitary matrix Q, and if +* COMPQ = 'V', the product Q1*Q. +* Not referenced if COMPQ='N'. +* +* LDQ (input) INTEGER +* The leading dimension of the array Q. +* LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. +* +* Z (input/output) COMPLEX array, dimension (LDZ, N) +* On entry, if COMPZ = 'V', the unitary matrix Z1. +* On exit, if COMPZ='I', the unitary matrix Z, and if +* COMPZ = 'V', the product Z1*Z. +* Not referenced if COMPZ='N'. +* +* LDZ (input) INTEGER +* The leading dimension of the array Z. +* LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* Further Details +* =============== +* +* This routine reduces A to Hessenberg and B to triangular form by +* an unblocked reduction, as described in _Matrix_Computations_, +* by Golub and van Loan (Johns Hopkins Press). +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX CONE, CZERO + PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), + $ CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL ILQ, ILZ + INTEGER ICOMPQ, ICOMPZ, JCOL, JROW + REAL C + COMPLEX CTEMP, S +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CLARTG, CLASET, CROT, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC CONJG, MAX +* .. +* .. Executable Statements .. +* +* Decode COMPQ +* + IF( LSAME( COMPQ, 'N' ) ) THEN + ILQ = .FALSE. + ICOMPQ = 1 + ELSE IF( LSAME( COMPQ, 'V' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 2 + ELSE IF( LSAME( COMPQ, 'I' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 3 + ELSE + ICOMPQ = 0 + END IF +* +* Decode COMPZ +* + IF( LSAME( COMPZ, 'N' ) ) THEN + ILZ = .FALSE. + ICOMPZ = 1 + ELSE IF( LSAME( COMPZ, 'V' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 2 + ELSE IF( LSAME( COMPZ, 'I' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 3 + ELSE + ICOMPZ = 0 + END IF +* +* Test the input parameters. +* + INFO = 0 + IF( ICOMPQ.LE.0 ) THEN + INFO = -1 + ELSE IF( ICOMPZ.LE.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( ILO.LT.1 ) THEN + INFO = -4 + ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -9 + ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN + INFO = -11 + ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN + INFO = -13 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGGHRD', -INFO ) + RETURN + END IF +* +* Initialize Q and Z if desired. +* + IF( ICOMPQ.EQ.3 ) + $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) + IF( ICOMPZ.EQ.3 ) + $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) +* +* Quick return if possible +* + IF( N.LE.1 ) + $ RETURN +* +* Zero out lower triangle of B +* + DO 20 JCOL = 1, N - 1 + DO 10 JROW = JCOL + 1, N + B( JROW, JCOL ) = CZERO + 10 CONTINUE + 20 CONTINUE +* +* Reduce A and B +* + DO 40 JCOL = ILO, IHI - 2 +* + DO 30 JROW = IHI, JCOL + 2, -1 +* +* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) +* + CTEMP = A( JROW-1, JCOL ) + CALL CLARTG( CTEMP, A( JROW, JCOL ), C, S, + $ A( JROW-1, JCOL ) ) + A( JROW, JCOL ) = CZERO + CALL CROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA, + $ A( JROW, JCOL+1 ), LDA, C, S ) + CALL CROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB, + $ B( JROW, JROW-1 ), LDB, C, S ) + IF( ILQ ) + $ CALL CROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, + $ CONJG( S ) ) +* +* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) +* + CTEMP = B( JROW, JROW ) + CALL CLARTG( CTEMP, B( JROW, JROW-1 ), C, S, + $ B( JROW, JROW ) ) + B( JROW, JROW-1 ) = CZERO + CALL CROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S ) + CALL CROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C, + $ S ) + IF( ILZ ) + $ CALL CROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S ) + 30 CONTINUE + 40 CONTINUE +* + RETURN +* +* End of CGGHRD +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/chgeqz.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,758 @@ + SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, + $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, + $ RWORK, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER COMPQ, COMPZ, JOB + INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N +* .. +* .. Array Arguments .. + REAL RWORK( * ) + COMPLEX ALPHA( * ), BETA( * ), H( LDH, * ), + $ Q( LDQ, * ), T( LDT, * ), WORK( * ), + $ Z( LDZ, * ) +* .. +* +* Purpose +* ======= +* +* CHGEQZ computes the eigenvalues of a complex matrix pair (H,T), +* where H is an upper Hessenberg matrix and T is upper triangular, +* using the single-shift QZ method. +* Matrix pairs of this type are produced by the reduction to +* generalized upper Hessenberg form of a complex matrix pair (A,B): +* +* A = Q1*H*Z1**H, B = Q1*T*Z1**H, +* +* as computed by CGGHRD. +* +* If JOB='S', then the Hessenberg-triangular pair (H,T) is +* also reduced to generalized Schur form, +* +* H = Q*S*Z**H, T = Q*P*Z**H, +* +* where Q and Z are unitary matrices and S and P are upper triangular. +* +* Optionally, the unitary matrix Q from the generalized Schur +* factorization may be postmultiplied into an input matrix Q1, and the +* unitary matrix Z may be postmultiplied into an input matrix Z1. +* If Q1 and Z1 are the unitary matrices from CGGHRD that reduced +* the matrix pair (A,B) to generalized Hessenberg form, then the output +* matrices Q1*Q and Z1*Z are the unitary factors from the generalized +* Schur factorization of (A,B): +* +* A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. +* +* To avoid overflow, eigenvalues of the matrix pair (H,T) +* (equivalently, of (A,B)) are computed as a pair of complex values +* (alpha,beta). If beta is nonzero, lambda = alpha / beta is an +* eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP) +* A*x = lambda*B*x +* and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the +* alternate form of the GNEP +* mu*A*y = B*y. +* The values of alpha and beta for the i-th eigenvalue can be read +* directly from the generalized Schur form: alpha = S(i,i), +* beta = P(i,i). +* +* Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix +* Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), +* pp. 241--256. +* +* Arguments +* ========= +* +* JOB (input) CHARACTER*1 +* = 'E': Compute eigenvalues only; +* = 'S': Computer eigenvalues and the Schur form. +* +* COMPQ (input) CHARACTER*1 +* = 'N': Left Schur vectors (Q) are not computed; +* = 'I': Q is initialized to the unit matrix and the matrix Q +* of left Schur vectors of (H,T) is returned; +* = 'V': Q must contain a unitary matrix Q1 on entry and +* the product Q1*Q is returned. +* +* COMPZ (input) CHARACTER*1 +* = 'N': Right Schur vectors (Z) are not computed; +* = 'I': Q is initialized to the unit matrix and the matrix Z +* of right Schur vectors of (H,T) is returned; +* = 'V': Z must contain a unitary matrix Z1 on entry and +* the product Z1*Z is returned. +* +* N (input) INTEGER +* The order of the matrices H, T, Q, and Z. N >= 0. +* +* ILO (input) INTEGER +* IHI (input) INTEGER +* ILO and IHI mark the rows and columns of H which are in +* Hessenberg form. It is assumed that A is already upper +* triangular in rows and columns 1:ILO-1 and IHI+1:N. +* If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. +* +* H (input/output) COMPLEX array, dimension (LDH, N) +* On entry, the N-by-N upper Hessenberg matrix H. +* On exit, if JOB = 'S', H contains the upper triangular +* matrix S from the generalized Schur factorization. +* If JOB = 'E', the diagonal of H matches that of S, but +* the rest of H is unspecified. +* +* LDH (input) INTEGER +* The leading dimension of the array H. LDH >= max( 1, N ). +* +* T (input/output) COMPLEX array, dimension (LDT, N) +* On entry, the N-by-N upper triangular matrix T. +* On exit, if JOB = 'S', T contains the upper triangular +* matrix P from the generalized Schur factorization. +* If JOB = 'E', the diagonal of T matches that of P, but +* the rest of T is unspecified. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= max( 1, N ). +* +* ALPHA (output) COMPLEX array, dimension (N) +* The complex scalars alpha that define the eigenvalues of +* GNEP. ALPHA(i) = S(i,i) in the generalized Schur +* factorization. +* +* BETA (output) COMPLEX array, dimension (N) +* The real non-negative scalars beta that define the +* eigenvalues of GNEP. BETA(i) = P(i,i) in the generalized +* Schur factorization. +* +* Together, the quantities alpha = ALPHA(j) and beta = BETA(j) +* represent the j-th eigenvalue of the matrix pair (A,B), in +* one of the forms lambda = alpha/beta or mu = beta/alpha. +* Since either lambda or mu may overflow, they should not, +* in general, be computed. +* +* Q (input/output) COMPLEX array, dimension (LDQ, N) +* On entry, if COMPZ = 'V', the unitary matrix Q1 used in the +* reduction of (A,B) to generalized Hessenberg form. +* On exit, if COMPZ = 'I', the unitary matrix of left Schur +* vectors of (H,T), and if COMPZ = 'V', the unitary matrix of +* left Schur vectors of (A,B). +* Not referenced if COMPZ = 'N'. +* +* LDQ (input) INTEGER +* The leading dimension of the array Q. LDQ >= 1. +* If COMPQ='V' or 'I', then LDQ >= N. +* +* Z (input/output) COMPLEX array, dimension (LDZ, N) +* On entry, if COMPZ = 'V', the unitary matrix Z1 used in the +* reduction of (A,B) to generalized Hessenberg form. +* On exit, if COMPZ = 'I', the unitary matrix of right Schur +* vectors of (H,T), and if COMPZ = 'V', the unitary matrix of +* right Schur vectors of (A,B). +* Not referenced if COMPZ = 'N'. +* +* LDZ (input) INTEGER +* The leading dimension of the array Z. LDZ >= 1. +* If COMPZ='V' or 'I', then LDZ >= N. +* +* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) +* On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* RWORK (workspace) REAL array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* = 1,...,N: the QZ iteration did not converge. (H,T) is not +* in Schur form, but ALPHA(i) and BETA(i), +* i=INFO+1,...,N should be correct. +* = N+1,...,2*N: the shift calculation failed. (H,T) is not +* in Schur form, but ALPHA(i) and BETA(i), +* i=INFO-N+1,...,N should be correct. +* +* Further Details +* =============== +* +* We assume that complex ABS works as long as its value is less than +* overflow. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), + $ CONE = ( 1.0E+0, 0.0E+0 ) ) + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + REAL HALF + PARAMETER ( HALF = 0.5E+0 ) +* .. +* .. Local Scalars .. + LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY + INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, + $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, + $ JR, MAXIT + REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, + $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP + COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, + $ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1, + $ U12, X +* .. +* .. External Functions .. + LOGICAL LSAME + REAL CLANHS, SLAMCH + EXTERNAL LSAME, CLANHS, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT +* .. +* .. Statement Functions .. + REAL ABS1 +* .. +* .. Statement Function definitions .. + ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) +* .. +* .. Executable Statements .. +* +* Decode JOB, COMPQ, COMPZ +* + IF( LSAME( JOB, 'E' ) ) THEN + ILSCHR = .FALSE. + ISCHUR = 1 + ELSE IF( LSAME( JOB, 'S' ) ) THEN + ILSCHR = .TRUE. + ISCHUR = 2 + ELSE + ISCHUR = 0 + END IF +* + IF( LSAME( COMPQ, 'N' ) ) THEN + ILQ = .FALSE. + ICOMPQ = 1 + ELSE IF( LSAME( COMPQ, 'V' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 2 + ELSE IF( LSAME( COMPQ, 'I' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 3 + ELSE + ICOMPQ = 0 + END IF +* + IF( LSAME( COMPZ, 'N' ) ) THEN + ILZ = .FALSE. + ICOMPZ = 1 + ELSE IF( LSAME( COMPZ, 'V' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 2 + ELSE IF( LSAME( COMPZ, 'I' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 3 + ELSE + ICOMPZ = 0 + END IF +* +* Check Argument Values +* + INFO = 0 + WORK( 1 ) = MAX( 1, N ) + LQUERY = ( LWORK.EQ.-1 ) + IF( ISCHUR.EQ.0 ) THEN + INFO = -1 + ELSE IF( ICOMPQ.EQ.0 ) THEN + INFO = -2 + ELSE IF( ICOMPZ.EQ.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( ILO.LT.1 ) THEN + INFO = -5 + ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN + INFO = -6 + ELSE IF( LDH.LT.N ) THEN + INFO = -8 + ELSE IF( LDT.LT.N ) THEN + INFO = -10 + ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN + INFO = -14 + ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN + INFO = -16 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -18 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CHGEQZ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* +c WORK( 1 ) = CMPLX( 1 ) + IF( N.LE.0 ) THEN + WORK( 1 ) = CMPLX( 1 ) + RETURN + END IF +* +* Initialize Q and Z +* + IF( ICOMPQ.EQ.3 ) + $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) + IF( ICOMPZ.EQ.3 ) + $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) +* +* Machine Constants +* + IN = IHI + 1 - ILO + SAFMIN = SLAMCH( 'S' ) + ULP = SLAMCH( 'E' )*SLAMCH( 'B' ) + ANORM = CLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK ) + BNORM = CLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK ) + ATOL = MAX( SAFMIN, ULP*ANORM ) + BTOL = MAX( SAFMIN, ULP*BNORM ) + ASCALE = ONE / MAX( SAFMIN, ANORM ) + BSCALE = ONE / MAX( SAFMIN, BNORM ) +* +* +* Set Eigenvalues IHI+1:N +* + DO 10 J = IHI + 1, N + ABSB = ABS( T( J, J ) ) + IF( ABSB.GT.SAFMIN ) THEN + SIGNBC = CONJG( T( J, J ) / ABSB ) + T( J, J ) = ABSB + IF( ILSCHR ) THEN + CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 ) + CALL CSCAL( J, SIGNBC, H( 1, J ), 1 ) + ELSE + H( J, J ) = H( J, J )*SIGNBC + END IF + IF( ILZ ) + $ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 ) + ELSE + T( J, J ) = CZERO + END IF + ALPHA( J ) = H( J, J ) + BETA( J ) = T( J, J ) + 10 CONTINUE +* +* If IHI < ILO, skip QZ steps +* + IF( IHI.LT.ILO ) + $ GO TO 190 +* +* MAIN QZ ITERATION LOOP +* +* Initialize dynamic indices +* +* Eigenvalues ILAST+1:N have been found. +* Column operations modify rows IFRSTM:whatever +* Row operations modify columns whatever:ILASTM +* +* If only eigenvalues are being computed, then +* IFRSTM is the row of the last splitting row above row ILAST; +* this is always at least ILO. +* IITER counts iterations since the last eigenvalue was found, +* to tell when to use an extraordinary shift. +* MAXIT is the maximum number of QZ sweeps allowed. +* + ILAST = IHI + IF( ILSCHR ) THEN + IFRSTM = 1 + ILASTM = N + ELSE + IFRSTM = ILO + ILASTM = IHI + END IF + IITER = 0 + ESHIFT = CZERO + MAXIT = 30*( IHI-ILO+1 ) +* + DO 170 JITER = 1, MAXIT +* +* Check for too many iterations. +* + IF( JITER.GT.MAXIT ) + $ GO TO 180 +* +* Split the matrix if possible. +* +* Two tests: +* 1: H(j,j-1)=0 or j=ILO +* 2: T(j,j)=0 +* +* Special case: j=ILAST +* + IF( ILAST.EQ.ILO ) THEN + GO TO 60 + ELSE + IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN + H( ILAST, ILAST-1 ) = CZERO + GO TO 60 + END IF + END IF +* + IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN + T( ILAST, ILAST ) = CZERO + GO TO 50 + END IF +* +* General case: j<ILAST +* + DO 40 J = ILAST - 1, ILO, -1 +* +* Test 1: for H(j,j-1)=0 or j=ILO +* + IF( J.EQ.ILO ) THEN + ILAZRO = .TRUE. + ELSE + IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN + H( J, J-1 ) = CZERO + ILAZRO = .TRUE. + ELSE + ILAZRO = .FALSE. + END IF + END IF +* +* Test 2: for T(j,j)=0 +* + IF( ABS( T( J, J ) ).LT.BTOL ) THEN + T( J, J ) = CZERO +* +* Test 1a: Check for 2 consecutive small subdiagonals in A +* + ILAZR2 = .FALSE. + IF( .NOT.ILAZRO ) THEN + IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1, + $ J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) ) + $ ILAZR2 = .TRUE. + END IF +* +* If both tests pass (1 & 2), i.e., the leading diagonal +* element of B in the block is zero, split a 1x1 block off +* at the top. (I.e., at the J-th row/column) The leading +* diagonal element of the remainder can also be zero, so +* this may have to be done repeatedly. +* + IF( ILAZRO .OR. ILAZR2 ) THEN + DO 20 JCH = J, ILAST - 1 + CTEMP = H( JCH, JCH ) + CALL CLARTG( CTEMP, H( JCH+1, JCH ), C, S, + $ H( JCH, JCH ) ) + H( JCH+1, JCH ) = CZERO + CALL CROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, + $ H( JCH+1, JCH+1 ), LDH, C, S ) + CALL CROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, + $ T( JCH+1, JCH+1 ), LDT, C, S ) + IF( ILQ ) + $ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, + $ C, CONJG( S ) ) + IF( ILAZR2 ) + $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C + ILAZR2 = .FALSE. + IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN + IF( JCH+1.GE.ILAST ) THEN + GO TO 60 + ELSE + IFIRST = JCH + 1 + GO TO 70 + END IF + END IF + T( JCH+1, JCH+1 ) = CZERO + 20 CONTINUE + GO TO 50 + ELSE +* +* Only test 2 passed -- chase the zero to T(ILAST,ILAST) +* Then process as in the case T(ILAST,ILAST)=0 +* + DO 30 JCH = J, ILAST - 1 + CTEMP = T( JCH, JCH+1 ) + CALL CLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S, + $ T( JCH, JCH+1 ) ) + T( JCH+1, JCH+1 ) = CZERO + IF( JCH.LT.ILASTM-1 ) + $ CALL CROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, + $ T( JCH+1, JCH+2 ), LDT, C, S ) + CALL CROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, + $ H( JCH+1, JCH-1 ), LDH, C, S ) + IF( ILQ ) + $ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, + $ C, CONJG( S ) ) + CTEMP = H( JCH+1, JCH ) + CALL CLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S, + $ H( JCH+1, JCH ) ) + H( JCH+1, JCH-1 ) = CZERO + CALL CROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, + $ H( IFRSTM, JCH-1 ), 1, C, S ) + CALL CROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, + $ T( IFRSTM, JCH-1 ), 1, C, S ) + IF( ILZ ) + $ CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, + $ C, S ) + 30 CONTINUE + GO TO 50 + END IF + ELSE IF( ILAZRO ) THEN +* +* Only test 1 passed -- work on J:ILAST +* + IFIRST = J + GO TO 70 + END IF +* +* Neither test passed -- try next J +* + 40 CONTINUE +* +* (Drop-through is "impossible") +* + INFO = 2*N + 1 + GO TO 210 +* +* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a +* 1x1 block. +* + 50 CONTINUE + CTEMP = H( ILAST, ILAST ) + CALL CLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S, + $ H( ILAST, ILAST ) ) + H( ILAST, ILAST-1 ) = CZERO + CALL CROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, + $ H( IFRSTM, ILAST-1 ), 1, C, S ) + CALL CROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, + $ T( IFRSTM, ILAST-1 ), 1, C, S ) + IF( ILZ ) + $ CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) +* +* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA +* + 60 CONTINUE + ABSB = ABS( T( ILAST, ILAST ) ) + IF( ABSB.GT.SAFMIN ) THEN + SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB ) + T( ILAST, ILAST ) = ABSB + IF( ILSCHR ) THEN + CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 ) + CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ), + $ 1 ) + ELSE + H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC + END IF + IF( ILZ ) + $ CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 ) + ELSE + T( ILAST, ILAST ) = CZERO + END IF + ALPHA( ILAST ) = H( ILAST, ILAST ) + BETA( ILAST ) = T( ILAST, ILAST ) +* +* Go to next block -- exit if finished. +* + ILAST = ILAST - 1 + IF( ILAST.LT.ILO ) + $ GO TO 190 +* +* Reset counters +* + IITER = 0 + ESHIFT = CZERO + IF( .NOT.ILSCHR ) THEN + ILASTM = ILAST + IF( IFRSTM.GT.ILAST ) + $ IFRSTM = ILO + END IF + GO TO 160 +* +* QZ step +* +* This iteration only involves rows/columns IFIRST:ILAST. We +* assume IFIRST < ILAST, and that the diagonal of B is non-zero. +* + 70 CONTINUE + IITER = IITER + 1 + IF( .NOT.ILSCHR ) THEN + IFRSTM = IFIRST + END IF +* +* Compute the Shift. +* +* At this point, IFIRST < ILAST, and the diagonal elements of +* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in +* magnitude) +* + IF( ( IITER / 10 )*10.NE.IITER ) THEN +* +* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of +* the bottom-right 2x2 block of A inv(B) which is nearest to +* the bottom-right element. +* +* We factor B as U*D, where U has unit diagonals, and +* compute (A*inv(D))*inv(U). +* + U12 = ( BSCALE*T( ILAST-1, ILAST ) ) / + $ ( BSCALE*T( ILAST, ILAST ) ) + AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / + $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) + AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / + $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) + AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / + $ ( BSCALE*T( ILAST, ILAST ) ) + AD22 = ( ASCALE*H( ILAST, ILAST ) ) / + $ ( BSCALE*T( ILAST, ILAST ) ) + ABI22 = AD22 - U12*AD21 +* + T1 = HALF*( AD11+ABI22 ) + RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 ) + TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) + + $ AIMAG( T1-ABI22 )*AIMAG( RTDISC ) + IF( TEMP.LE.ZERO ) THEN + SHIFT = T1 + RTDISC + ELSE + SHIFT = T1 - RTDISC + END IF + ELSE +* +* Exceptional shift. Chosen for no particularly good reason. +* + ESHIFT = ESHIFT + CONJG( ( ASCALE*H( ILAST-1, ILAST ) ) / + $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) ) + SHIFT = ESHIFT + END IF +* +* Now check for two consecutive small subdiagonals. +* + DO 80 J = ILAST - 1, IFIRST + 1, -1 + ISTART = J + CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) ) + TEMP = ABS1( CTEMP ) + TEMP2 = ASCALE*ABS1( H( J+1, J ) ) + TEMPR = MAX( TEMP, TEMP2 ) + IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN + TEMP = TEMP / TEMPR + TEMP2 = TEMP2 / TEMPR + END IF + IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL ) + $ GO TO 90 + 80 CONTINUE +* + ISTART = IFIRST + CTEMP = ASCALE*H( IFIRST, IFIRST ) - + $ SHIFT*( BSCALE*T( IFIRST, IFIRST ) ) + 90 CONTINUE +* +* Do an implicit-shift QZ sweep. +* +* Initial Q +* + CTEMP2 = ASCALE*H( ISTART+1, ISTART ) + CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) +* +* Sweep +* + DO 150 J = ISTART, ILAST - 1 + IF( J.GT.ISTART ) THEN + CTEMP = H( J, J-1 ) + CALL CLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) + H( J+1, J-1 ) = CZERO + END IF +* + DO 100 JC = J, ILASTM + CTEMP = C*H( J, JC ) + S*H( J+1, JC ) + H( J+1, JC ) = -CONJG( S )*H( J, JC ) + C*H( J+1, JC ) + H( J, JC ) = CTEMP + CTEMP2 = C*T( J, JC ) + S*T( J+1, JC ) + T( J+1, JC ) = -CONJG( S )*T( J, JC ) + C*T( J+1, JC ) + T( J, JC ) = CTEMP2 + 100 CONTINUE + IF( ILQ ) THEN + DO 110 JR = 1, N + CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 ) + Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) + Q( JR, J ) = CTEMP + 110 CONTINUE + END IF +* + CTEMP = T( J+1, J+1 ) + CALL CLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) + T( J+1, J ) = CZERO +* + DO 120 JR = IFRSTM, MIN( J+2, ILAST ) + CTEMP = C*H( JR, J+1 ) + S*H( JR, J ) + H( JR, J ) = -CONJG( S )*H( JR, J+1 ) + C*H( JR, J ) + H( JR, J+1 ) = CTEMP + 120 CONTINUE + DO 130 JR = IFRSTM, J + CTEMP = C*T( JR, J+1 ) + S*T( JR, J ) + T( JR, J ) = -CONJG( S )*T( JR, J+1 ) + C*T( JR, J ) + T( JR, J+1 ) = CTEMP + 130 CONTINUE + IF( ILZ ) THEN + DO 140 JR = 1, N + CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) + Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J ) + Z( JR, J+1 ) = CTEMP + 140 CONTINUE + END IF + 150 CONTINUE +* + 160 CONTINUE +* + 170 CONTINUE +* +* Drop-through = non-convergence +* + 180 CONTINUE + INFO = ILAST + GO TO 210 +* +* Successful completion of all QZ steps +* + 190 CONTINUE +* +* Set Eigenvalues 1:ILO-1 +* + DO 200 J = 1, ILO - 1 + ABSB = ABS( T( J, J ) ) + IF( ABSB.GT.SAFMIN ) THEN + SIGNBC = CONJG( T( J, J ) / ABSB ) + T( J, J ) = ABSB + IF( ILSCHR ) THEN + CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 ) + CALL CSCAL( J, SIGNBC, H( 1, J ), 1 ) + ELSE + H( J, J ) = H( J, J )*SIGNBC + END IF + IF( ILZ ) + $ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 ) + ELSE + T( J, J ) = CZERO + END IF + ALPHA( J ) = H( J, J ) + BETA( J ) = T( J, J ) + 200 CONTINUE +* +* Normal Termination +* + INFO = 0 +* +* Exit (other than argument error) -- return optimal workspace size +* + 210 CONTINUE + WORK( 1 ) = CMPLX( N ) + RETURN +* +* End of CHGEQZ +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/ctgevc.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,633 @@ + SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, + $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER HOWMNY, SIDE + INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N +* .. +* .. Array Arguments .. + LOGICAL SELECT( * ) + REAL RWORK( * ) + COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ), + $ VR( LDVR, * ), WORK( * ) +* .. +* +* +* Purpose +* ======= +* +* CTGEVC computes some or all of the right and/or left eigenvectors of +* a pair of complex matrices (S,P), where S and P are upper triangular. +* Matrix pairs of this type are produced by the generalized Schur +* factorization of a complex matrix pair (A,B): +* +* A = Q*S*Z**H, B = Q*P*Z**H +* +* as computed by CGGHRD + CHGEQZ. +* +* The right eigenvector x and the left eigenvector y of (S,P) +* corresponding to an eigenvalue w are defined by: +* +* S*x = w*P*x, (y**H)*S = w*(y**H)*P, +* +* where y**H denotes the conjugate tranpose of y. +* The eigenvalues are not input to this routine, but are computed +* directly from the diagonal elements of S and P. +* +* This routine returns the matrices X and/or Y of right and left +* eigenvectors of (S,P), or the products Z*X and/or Q*Y, +* where Z and Q are input matrices. +* If Q and Z are the unitary factors from the generalized Schur +* factorization of a matrix pair (A,B), then Z*X and Q*Y +* are the matrices of right and left eigenvectors of (A,B). +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'R': compute right eigenvectors only; +* = 'L': compute left eigenvectors only; +* = 'B': compute both right and left eigenvectors. +* +* HOWMNY (input) CHARACTER*1 +* = 'A': compute all right and/or left eigenvectors; +* = 'B': compute all right and/or left eigenvectors, +* backtransformed by the matrices in VR and/or VL; +* = 'S': compute selected right and/or left eigenvectors, +* specified by the logical array SELECT. +* +* SELECT (input) LOGICAL array, dimension (N) +* If HOWMNY='S', SELECT specifies the eigenvectors to be +* computed. The eigenvector corresponding to the j-th +* eigenvalue is computed if SELECT(j) = .TRUE.. +* Not referenced if HOWMNY = 'A' or 'B'. +* +* N (input) INTEGER +* The order of the matrices S and P. N >= 0. +* +* S (input) COMPLEX array, dimension (LDS,N) +* The upper triangular matrix S from a generalized Schur +* factorization, as computed by CHGEQZ. +* +* LDS (input) INTEGER +* The leading dimension of array S. LDS >= max(1,N). +* +* P (input) COMPLEX array, dimension (LDP,N) +* The upper triangular matrix P from a generalized Schur +* factorization, as computed by CHGEQZ. P must have real +* diagonal elements. +* +* LDP (input) INTEGER +* The leading dimension of array P. LDP >= max(1,N). +* +* VL (input/output) COMPLEX array, dimension (LDVL,MM) +* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must +* contain an N-by-N matrix Q (usually the unitary matrix Q +* of left Schur vectors returned by CHGEQZ). +* On exit, if SIDE = 'L' or 'B', VL contains: +* if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); +* if HOWMNY = 'B', the matrix Q*Y; +* if HOWMNY = 'S', the left eigenvectors of (S,P) specified by +* SELECT, stored consecutively in the columns of +* VL, in the same order as their eigenvalues. +* Not referenced if SIDE = 'R'. +* +* LDVL (input) INTEGER +* The leading dimension of array VL. LDVL >= 1, and if +* SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N. +* +* VR (input/output) COMPLEX array, dimension (LDVR,MM) +* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must +* contain an N-by-N matrix Q (usually the unitary matrix Z +* of right Schur vectors returned by CHGEQZ). +* On exit, if SIDE = 'R' or 'B', VR contains: +* if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); +* if HOWMNY = 'B', the matrix Z*X; +* if HOWMNY = 'S', the right eigenvectors of (S,P) specified by +* SELECT, stored consecutively in the columns of +* VR, in the same order as their eigenvalues. +* Not referenced if SIDE = 'L'. +* +* LDVR (input) INTEGER +* The leading dimension of the array VR. LDVR >= 1, and if +* SIDE = 'R' or 'B', LDVR >= N. +* +* MM (input) INTEGER +* The number of columns in the arrays VL and/or VR. MM >= M. +* +* M (output) INTEGER +* The number of columns in the arrays VL and/or VR actually +* used to store the eigenvectors. If HOWMNY = 'A' or 'B', M +* is set to N. Each selected eigenvector occupies one column. +* +* WORK (workspace) COMPLEX array, dimension (2*N) +* +* RWORK (workspace) REAL array, dimension (2*N) +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), + $ CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP, + $ LSA, LSB + INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC, + $ J, JE, JR + REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG, + $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA, + $ SCALE, SMALL, TEMP, ULP, XMAX + COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X +* .. +* .. External Functions .. + LOGICAL LSAME + REAL SLAMCH + COMPLEX CLADIV + EXTERNAL LSAME, SLAMCH, CLADIV +* .. +* .. External Subroutines .. + EXTERNAL CGEMV, SLABAD, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL +* .. +* .. Statement Functions .. + REAL ABS1 +* .. +* .. Statement Function definitions .. + ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) +* .. +* .. Executable Statements .. +* +* Decode and Test the input parameters +* + IF( LSAME( HOWMNY, 'A' ) ) THEN + IHWMNY = 1 + ILALL = .TRUE. + ILBACK = .FALSE. + ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN + IHWMNY = 2 + ILALL = .FALSE. + ILBACK = .FALSE. + ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN + IHWMNY = 3 + ILALL = .TRUE. + ILBACK = .TRUE. + ELSE + IHWMNY = -1 + END IF +* + IF( LSAME( SIDE, 'R' ) ) THEN + ISIDE = 1 + COMPL = .FALSE. + COMPR = .TRUE. + ELSE IF( LSAME( SIDE, 'L' ) ) THEN + ISIDE = 2 + COMPL = .TRUE. + COMPR = .FALSE. + ELSE IF( LSAME( SIDE, 'B' ) ) THEN + ISIDE = 3 + COMPL = .TRUE. + COMPR = .TRUE. + ELSE + ISIDE = -1 + END IF +* + INFO = 0 + IF( ISIDE.LT.0 ) THEN + INFO = -1 + ELSE IF( IHWMNY.LT.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( LDS.LT.MAX( 1, N ) ) THEN + INFO = -6 + ELSE IF( LDP.LT.MAX( 1, N ) ) THEN + INFO = -8 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CTGEVC', -INFO ) + RETURN + END IF +* +* Count the number of eigenvectors +* + IF( .NOT.ILALL ) THEN + IM = 0 + DO 10 J = 1, N + IF( SELECT( J ) ) + $ IM = IM + 1 + 10 CONTINUE + ELSE + IM = N + END IF +* +* Check diagonal of B +* + ILBBAD = .FALSE. + DO 20 J = 1, N + IF( AIMAG( P( J, J ) ).NE.ZERO ) + $ ILBBAD = .TRUE. + 20 CONTINUE +* + IF( ILBBAD ) THEN + INFO = -7 + ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN + INFO = -10 + ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN + INFO = -12 + ELSE IF( MM.LT.IM ) THEN + INFO = -13 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CTGEVC', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + M = IM + IF( N.EQ.0 ) + $ RETURN +* +* Machine Constants +* + SAFMIN = SLAMCH( 'Safe minimum' ) + BIG = ONE / SAFMIN + CALL SLABAD( SAFMIN, BIG ) + ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) + SMALL = SAFMIN*N / ULP + BIG = ONE / SMALL + BIGNUM = ONE / ( SAFMIN*N ) +* +* Compute the 1-norm of each column of the strictly upper triangular +* part of A and B to check for possible overflow in the triangular +* solver. +* + ANORM = ABS1( S( 1, 1 ) ) + BNORM = ABS1( P( 1, 1 ) ) + RWORK( 1 ) = ZERO + RWORK( N+1 ) = ZERO + DO 40 J = 2, N + RWORK( J ) = ZERO + RWORK( N+J ) = ZERO + DO 30 I = 1, J - 1 + RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) ) + RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) ) + 30 CONTINUE + ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) ) + BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) ) + 40 CONTINUE +* + ASCALE = ONE / MAX( ANORM, SAFMIN ) + BSCALE = ONE / MAX( BNORM, SAFMIN ) +* +* Left eigenvectors +* + IF( COMPL ) THEN + IEIG = 0 +* +* Main loop over eigenvalues +* + DO 140 JE = 1, N + IF( ILALL ) THEN + ILCOMP = .TRUE. + ELSE + ILCOMP = SELECT( JE ) + END IF + IF( ILCOMP ) THEN + IEIG = IEIG + 1 +* + IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND. + $ ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN +* +* Singular matrix pencil -- return unit eigenvector +* + DO 50 JR = 1, N + VL( JR, IEIG ) = CZERO + 50 CONTINUE + VL( IEIG, IEIG ) = CONE + GO TO 140 + END IF +* +* Non-singular eigenvalue: +* Compute coefficients a and b in +* H +* y ( a A - b B ) = 0 +* + TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE, + $ ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN ) + SALPHA = ( TEMP*S( JE, JE ) )*ASCALE + SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE + ACOEFF = SBETA*ASCALE + BCOEFF = SALPHA*BSCALE +* +* Scale to avoid underflow +* + LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL + LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. + $ SMALL +* + SCALE = ONE + IF( LSA ) + $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) + IF( LSB ) + $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* + $ MIN( BNORM, BIG ) ) + IF( LSA .OR. LSB ) THEN + SCALE = MIN( SCALE, ONE / + $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), + $ ABS1( BCOEFF ) ) ) ) + IF( LSA ) THEN + ACOEFF = ASCALE*( SCALE*SBETA ) + ELSE + ACOEFF = SCALE*ACOEFF + END IF + IF( LSB ) THEN + BCOEFF = BSCALE*( SCALE*SALPHA ) + ELSE + BCOEFF = SCALE*BCOEFF + END IF + END IF +* + ACOEFA = ABS( ACOEFF ) + BCOEFA = ABS1( BCOEFF ) + XMAX = ONE + DO 60 JR = 1, N + WORK( JR ) = CZERO + 60 CONTINUE + WORK( JE ) = CONE + DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) +* +* H +* Triangular solve of (a A - b B) y = 0 +* +* H +* (rowwise in (a A - b B) , or columnwise in a A - b B) +* + DO 100 J = JE + 1, N +* +* Compute +* j-1 +* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) +* k=je +* (Scale if necessary) +* + TEMP = ONE / XMAX + IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM* + $ TEMP ) THEN + DO 70 JR = JE, J - 1 + WORK( JR ) = TEMP*WORK( JR ) + 70 CONTINUE + XMAX = ONE + END IF + SUMA = CZERO + SUMB = CZERO +* + DO 80 JR = JE, J - 1 + SUMA = SUMA + CONJG( S( JR, J ) )*WORK( JR ) + SUMB = SUMB + CONJG( P( JR, J ) )*WORK( JR ) + 80 CONTINUE + SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB +* +* Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) ) +* +* with scaling and perturbation of the denominator +* + D = CONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) ) + IF( ABS1( D ).LE.DMIN ) + $ D = CMPLX( DMIN ) +* + IF( ABS1( D ).LT.ONE ) THEN + IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN + TEMP = ONE / ABS1( SUM ) + DO 90 JR = JE, J - 1 + WORK( JR ) = TEMP*WORK( JR ) + 90 CONTINUE + XMAX = TEMP*XMAX + SUM = TEMP*SUM + END IF + END IF + WORK( J ) = CLADIV( -SUM, D ) + XMAX = MAX( XMAX, ABS1( WORK( J ) ) ) + 100 CONTINUE +* +* Back transform eigenvector if HOWMNY='B'. +* + IF( ILBACK ) THEN + CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL, + $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 ) + ISRC = 2 + IBEG = 1 + ELSE + ISRC = 1 + IBEG = JE + END IF +* +* Copy and scale eigenvector into column of VL +* + XMAX = ZERO + DO 110 JR = IBEG, N + XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) + 110 CONTINUE +* + IF( XMAX.GT.SAFMIN ) THEN + TEMP = ONE / XMAX + DO 120 JR = IBEG, N + VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) + 120 CONTINUE + ELSE + IBEG = N + 1 + END IF +* + DO 130 JR = 1, IBEG - 1 + VL( JR, IEIG ) = CZERO + 130 CONTINUE +* + END IF + 140 CONTINUE + END IF +* +* Right eigenvectors +* + IF( COMPR ) THEN + IEIG = IM + 1 +* +* Main loop over eigenvalues +* + DO 250 JE = N, 1, -1 + IF( ILALL ) THEN + ILCOMP = .TRUE. + ELSE + ILCOMP = SELECT( JE ) + END IF + IF( ILCOMP ) THEN + IEIG = IEIG - 1 +* + IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND. + $ ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN +* +* Singular matrix pencil -- return unit eigenvector +* + DO 150 JR = 1, N + VR( JR, IEIG ) = CZERO + 150 CONTINUE + VR( IEIG, IEIG ) = CONE + GO TO 250 + END IF +* +* Non-singular eigenvalue: +* Compute coefficients a and b in +* +* ( a A - b B ) x = 0 +* + TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE, + $ ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN ) + SALPHA = ( TEMP*S( JE, JE ) )*ASCALE + SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE + ACOEFF = SBETA*ASCALE + BCOEFF = SALPHA*BSCALE +* +* Scale to avoid underflow +* + LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL + LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. + $ SMALL +* + SCALE = ONE + IF( LSA ) + $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) + IF( LSB ) + $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* + $ MIN( BNORM, BIG ) ) + IF( LSA .OR. LSB ) THEN + SCALE = MIN( SCALE, ONE / + $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), + $ ABS1( BCOEFF ) ) ) ) + IF( LSA ) THEN + ACOEFF = ASCALE*( SCALE*SBETA ) + ELSE + ACOEFF = SCALE*ACOEFF + END IF + IF( LSB ) THEN + BCOEFF = BSCALE*( SCALE*SALPHA ) + ELSE + BCOEFF = SCALE*BCOEFF + END IF + END IF +* + ACOEFA = ABS( ACOEFF ) + BCOEFA = ABS1( BCOEFF ) + XMAX = ONE + DO 160 JR = 1, N + WORK( JR ) = CZERO + 160 CONTINUE + WORK( JE ) = CONE + DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) +* +* Triangular solve of (a A - b B) x = 0 (columnwise) +* +* WORK(1:j-1) contains sums w, +* WORK(j+1:JE) contains x +* + DO 170 JR = 1, JE - 1 + WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE ) + 170 CONTINUE + WORK( JE ) = CONE +* + DO 210 J = JE - 1, 1, -1 +* +* Form x(j) := - w(j) / d +* with scaling and perturbation of the denominator +* + D = ACOEFF*S( J, J ) - BCOEFF*P( J, J ) + IF( ABS1( D ).LE.DMIN ) + $ D = CMPLX( DMIN ) +* + IF( ABS1( D ).LT.ONE ) THEN + IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN + TEMP = ONE / ABS1( WORK( J ) ) + DO 180 JR = 1, JE + WORK( JR ) = TEMP*WORK( JR ) + 180 CONTINUE + END IF + END IF +* + WORK( J ) = CLADIV( -WORK( J ), D ) +* + IF( J.GT.1 ) THEN +* +* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling +* + IF( ABS1( WORK( J ) ).GT.ONE ) THEN + TEMP = ONE / ABS1( WORK( J ) ) + IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE. + $ BIGNUM*TEMP ) THEN + DO 190 JR = 1, JE + WORK( JR ) = TEMP*WORK( JR ) + 190 CONTINUE + END IF + END IF +* + CA = ACOEFF*WORK( J ) + CB = BCOEFF*WORK( J ) + DO 200 JR = 1, J - 1 + WORK( JR ) = WORK( JR ) + CA*S( JR, J ) - + $ CB*P( JR, J ) + 200 CONTINUE + END IF + 210 CONTINUE +* +* Back transform eigenvector if HOWMNY='B'. +* + IF( ILBACK ) THEN + CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1, + $ CZERO, WORK( N+1 ), 1 ) + ISRC = 2 + IEND = N + ELSE + ISRC = 1 + IEND = JE + END IF +* +* Copy and scale eigenvector into column of VR +* + XMAX = ZERO + DO 220 JR = 1, IEND + XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) + 220 CONTINUE +* + IF( XMAX.GT.SAFMIN ) THEN + TEMP = ONE / XMAX + DO 230 JR = 1, IEND + VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) + 230 CONTINUE + ELSE + IEND = 0 + END IF +* + DO 240 JR = IEND + 1, N + VR( JR, IEIG ) = CZERO + 240 CONTINUE +* + END IF + 250 CONTINUE + END IF +* + RETURN +* +* End of CTGEVC +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dggev.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,489 @@ + SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, + $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) +* +* -- LAPACK driver routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOBVL, JOBVR + INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), + $ B( LDB, * ), BETA( * ), VL( LDVL, * ), + $ VR( LDVR, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B) +* the generalized eigenvalues, and optionally, the left and/or right +* generalized eigenvectors. +* +* A generalized eigenvalue for a pair of matrices (A,B) is a scalar +* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is +* singular. It is usually represented as the pair (alpha,beta), as +* there is a reasonable interpretation for beta=0, and even for both +* being zero. +* +* The right eigenvector v(j) corresponding to the eigenvalue lambda(j) +* of (A,B) satisfies +* +* A * v(j) = lambda(j) * B * v(j). +* +* The left eigenvector u(j) corresponding to the eigenvalue lambda(j) +* of (A,B) satisfies +* +* u(j)**H * A = lambda(j) * u(j)**H * B . +* +* where u(j)**H is the conjugate-transpose of u(j). +* +* +* Arguments +* ========= +* +* JOBVL (input) CHARACTER*1 +* = 'N': do not compute the left generalized eigenvectors; +* = 'V': compute the left generalized eigenvectors. +* +* JOBVR (input) CHARACTER*1 +* = 'N': do not compute the right generalized eigenvectors; +* = 'V': compute the right generalized eigenvectors. +* +* N (input) INTEGER +* The order of the matrices A, B, VL, and VR. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) +* On entry, the matrix A in the pair (A,B). +* On exit, A has been overwritten. +* +* LDA (input) INTEGER +* The leading dimension of A. LDA >= max(1,N). +* +* B (input/output) DOUBLE PRECISION array, dimension (LDB, N) +* On entry, the matrix B in the pair (A,B). +* On exit, B has been overwritten. +* +* LDB (input) INTEGER +* The leading dimension of B. LDB >= max(1,N). +* +* ALPHAR (output) DOUBLE PRECISION array, dimension (N) +* ALPHAI (output) DOUBLE PRECISION array, dimension (N) +* BETA (output) DOUBLE PRECISION array, dimension (N) +* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will +* be the generalized eigenvalues. If ALPHAI(j) is zero, then +* the j-th eigenvalue is real; if positive, then the j-th and +* (j+1)-st eigenvalues are a complex conjugate pair, with +* ALPHAI(j+1) negative. +* +* Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) +* may easily over- or underflow, and BETA(j) may even be zero. +* Thus, the user should avoid naively computing the ratio +* alpha/beta. However, ALPHAR and ALPHAI will be always less +* than and usually comparable with norm(A) in magnitude, and +* BETA always less than and usually comparable with norm(B). +* +* VL (output) DOUBLE PRECISION array, dimension (LDVL,N) +* If JOBVL = 'V', the left eigenvectors u(j) are stored one +* after another in the columns of VL, in the same order as +* their eigenvalues. If the j-th eigenvalue is real, then +* u(j) = VL(:,j), the j-th column of VL. If the j-th and +* (j+1)-th eigenvalues form a complex conjugate pair, then +* u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1). +* Each eigenvector is scaled so the largest component has +* abs(real part)+abs(imag. part)=1. +* Not referenced if JOBVL = 'N'. +* +* LDVL (input) INTEGER +* The leading dimension of the matrix VL. LDVL >= 1, and +* if JOBVL = 'V', LDVL >= N. +* +* VR (output) DOUBLE PRECISION array, dimension (LDVR,N) +* If JOBVR = 'V', the right eigenvectors v(j) are stored one +* after another in the columns of VR, in the same order as +* their eigenvalues. If the j-th eigenvalue is real, then +* v(j) = VR(:,j), the j-th column of VR. If the j-th and +* (j+1)-th eigenvalues form a complex conjugate pair, then +* v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1). +* Each eigenvector is scaled so the largest component has +* abs(real part)+abs(imag. part)=1. +* Not referenced if JOBVR = 'N'. +* +* LDVR (input) INTEGER +* The leading dimension of the matrix VR. LDVR >= 1, and +* if JOBVR = 'V', LDVR >= N. +* +* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,8*N). +* For good performance, LWORK must generally be larger. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value. +* = 1,...,N: +* The QZ iteration failed. No eigenvectors have been +* calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) +* should be correct for j=INFO+1,...,N. +* > N: =N+1: other than QZ iteration failed in DHGEQZ. +* =N+2: error return from DTGEVC. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY + CHARACTER CHTEMP + INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO, + $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK, + $ MINWRK + DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, + $ SMLNUM, TEMP +* .. +* .. Local Arrays .. + LOGICAL LDUMMA( 1 ) +* .. +* .. External Subroutines .. + EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD, + $ DLACPY,DLASCL, DLASET, DORGQR, DORMQR, DTGEVC, + $ XERBLA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + DOUBLE PRECISION DLAMCH, DLANGE + EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SQRT +* .. +* .. Executable Statements .. +* +* Decode the input arguments +* + IF( LSAME( JOBVL, 'N' ) ) THEN + IJOBVL = 1 + ILVL = .FALSE. + ELSE IF( LSAME( JOBVL, 'V' ) ) THEN + IJOBVL = 2 + ILVL = .TRUE. + ELSE + IJOBVL = -1 + ILVL = .FALSE. + END IF +* + IF( LSAME( JOBVR, 'N' ) ) THEN + IJOBVR = 1 + ILVR = .FALSE. + ELSE IF( LSAME( JOBVR, 'V' ) ) THEN + IJOBVR = 2 + ILVR = .TRUE. + ELSE + IJOBVR = -1 + ILVR = .FALSE. + END IF + ILV = ILVL .OR. ILVR +* +* Test the input arguments +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( IJOBVL.LE.0 ) THEN + INFO = -1 + ELSE IF( IJOBVR.LE.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN + INFO = -12 + ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN + INFO = -14 + END IF +* +* Compute workspace +* (Note: Comments in the code beginning "Workspace:" describe the +* minimal amount of workspace needed at that point in the code, +* as well as the preferred amount for good performance. +* NB refers to the optimal block size for the immediately +* following subroutine, as returned by ILAENV. The workspace is +* computed assuming ILO = 1 and IHI = N, the worst case.) +* + IF( INFO.EQ.0 ) THEN + MINWRK = MAX( 1, 8*N ) + MAXWRK = MAX( 1, N*( 7 + + $ ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) ) ) + MAXWRK = MAX( MAXWRK, N*( 7 + + $ ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) ) ) + IF( ILVL ) THEN + MAXWRK = MAX( MAXWRK, N*( 7 + + $ ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) ) + END IF + WORK( 1 ) = MAXWRK +* + IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) + $ INFO = -16 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGGEV ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Get machine constants +* + EPS = DLAMCH( 'P' ) + SMLNUM = DLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) + SMLNUM = SQRT( SMLNUM ) / EPS + BIGNUM = ONE / SMLNUM +* +* Scale A if max element outside range [SMLNUM,BIGNUM] +* + ANRM = DLANGE( 'M', N, N, A, LDA, WORK ) + ILASCL = .FALSE. + IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN + ANRMTO = SMLNUM + ILASCL = .TRUE. + ELSE IF( ANRM.GT.BIGNUM ) THEN + ANRMTO = BIGNUM + ILASCL = .TRUE. + END IF + IF( ILASCL ) + $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) +* +* Scale B if max element outside range [SMLNUM,BIGNUM] +* + BNRM = DLANGE( 'M', N, N, B, LDB, WORK ) + ILBSCL = .FALSE. + IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN + BNRMTO = SMLNUM + ILBSCL = .TRUE. + ELSE IF( BNRM.GT.BIGNUM ) THEN + BNRMTO = BIGNUM + ILBSCL = .TRUE. + END IF + IF( ILBSCL ) + $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) +* +* Permute the matrices A, B to isolate eigenvalues if possible +* (Workspace: need 6*N) +* + ILEFT = 1 + IRIGHT = N + 1 + IWRK = IRIGHT + N + CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ), + $ WORK( IRIGHT ), WORK( IWRK ), IERR ) +* +* Reduce B to triangular form (QR decomposition of B) +* (Workspace: need N, prefer N*NB) +* + IROWS = IHI + 1 - ILO + IF( ILV ) THEN + ICOLS = N + 1 - ILO + ELSE + ICOLS = IROWS + END IF + ITAU = IWRK + IWRK = ITAU + IROWS + CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) +* +* Apply the orthogonal transformation to matrix A +* (Workspace: need N, prefer N*NB) +* + CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, + $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), + $ LWORK+1-IWRK, IERR ) +* +* Initialize VL +* (Workspace: need N, prefer N*NB) +* + IF( ILVL ) THEN + CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL ) + IF( IROWS.GT.1 ) THEN + CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, + $ VL( ILO+1, ILO ), LDVL ) + END IF + CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, + $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) + END IF +* +* Initialize VR +* + IF( ILVR ) + $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR ) +* +* Reduce to generalized Hessenberg form +* (Workspace: none needed) +* + IF( ILV ) THEN +* +* Eigenvectors requested -- work on whole matrix. +* + CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, + $ LDVL, VR, LDVR, IERR ) + ELSE + CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, + $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR ) + END IF +* +* Perform QZ algorithm (Compute eigenvalues, and optionally, the +* Schur forms and Schur vectors) +* (Workspace: need N) +* + IWRK = ITAU + IF( ILV ) THEN + CHTEMP = 'S' + ELSE + CHTEMP = 'E' + END IF + CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, + $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) + IF( IERR.NE.0 ) THEN + IF( IERR.GT.0 .AND. IERR.LE.N ) THEN + INFO = IERR + ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN + INFO = IERR - N + ELSE + INFO = N + 1 + END IF + GO TO 110 + END IF +* +* Compute Eigenvectors +* (Workspace: need 6*N) +* + IF( ILV ) THEN + IF( ILVL ) THEN + IF( ILVR ) THEN + CHTEMP = 'B' + ELSE + CHTEMP = 'L' + END IF + ELSE + CHTEMP = 'R' + END IF + CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, + $ VR, LDVR, N, IN, WORK( IWRK ), IERR ) + IF( IERR.NE.0 ) THEN + INFO = N + 2 + GO TO 110 + END IF +* +* Undo balancing on VL and VR and normalization +* (Workspace: none needed) +* + IF( ILVL ) THEN + CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ), + $ WORK( IRIGHT ), N, VL, LDVL, IERR ) + DO 50 JC = 1, N + IF( ALPHAI( JC ).LT.ZERO ) + $ GO TO 50 + TEMP = ZERO + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 10 JR = 1, N + TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) ) + 10 CONTINUE + ELSE + DO 20 JR = 1, N + TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+ + $ ABS( VL( JR, JC+1 ) ) ) + 20 CONTINUE + END IF + IF( TEMP.LT.SMLNUM ) + $ GO TO 50 + TEMP = ONE / TEMP + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 30 JR = 1, N + VL( JR, JC ) = VL( JR, JC )*TEMP + 30 CONTINUE + ELSE + DO 40 JR = 1, N + VL( JR, JC ) = VL( JR, JC )*TEMP + VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP + 40 CONTINUE + END IF + 50 CONTINUE + END IF + IF( ILVR ) THEN + CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ), + $ WORK( IRIGHT ), N, VR, LDVR, IERR ) + DO 100 JC = 1, N + IF( ALPHAI( JC ).LT.ZERO ) + $ GO TO 100 + TEMP = ZERO + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 60 JR = 1, N + TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) ) + 60 CONTINUE + ELSE + DO 70 JR = 1, N + TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+ + $ ABS( VR( JR, JC+1 ) ) ) + 70 CONTINUE + END IF + IF( TEMP.LT.SMLNUM ) + $ GO TO 100 + TEMP = ONE / TEMP + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 80 JR = 1, N + VR( JR, JC ) = VR( JR, JC )*TEMP + 80 CONTINUE + ELSE + DO 90 JR = 1, N + VR( JR, JC ) = VR( JR, JC )*TEMP + VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP + 90 CONTINUE + END IF + 100 CONTINUE + END IF +* +* End of eigenvector calculation +* + END IF +* +* Undo scaling if necessary +* + IF( ILASCL ) THEN + CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR ) + CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR ) + END IF +* + IF( ILBSCL ) THEN + CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) + END IF +* + 110 CONTINUE +* + WORK( 1 ) = MAXWRK +* + RETURN +* +* End of DGGEV +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/sggev.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,489 @@ + SUBROUTINE SGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, + $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) +* +* -- LAPACK driver routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOBVL, JOBVR + INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ), + $ B( LDB, * ), BETA( * ), VL( LDVL, * ), + $ VR( LDVR, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* SGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B) +* the generalized eigenvalues, and optionally, the left and/or right +* generalized eigenvectors. +* +* A generalized eigenvalue for a pair of matrices (A,B) is a scalar +* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is +* singular. It is usually represented as the pair (alpha,beta), as +* there is a reasonable interpretation for beta=0, and even for both +* being zero. +* +* The right eigenvector v(j) corresponding to the eigenvalue lambda(j) +* of (A,B) satisfies +* +* A * v(j) = lambda(j) * B * v(j). +* +* The left eigenvector u(j) corresponding to the eigenvalue lambda(j) +* of (A,B) satisfies +* +* u(j)**H * A = lambda(j) * u(j)**H * B . +* +* where u(j)**H is the conjugate-transpose of u(j). +* +* +* Arguments +* ========= +* +* JOBVL (input) CHARACTER*1 +* = 'N': do not compute the left generalized eigenvectors; +* = 'V': compute the left generalized eigenvectors. +* +* JOBVR (input) CHARACTER*1 +* = 'N': do not compute the right generalized eigenvectors; +* = 'V': compute the right generalized eigenvectors. +* +* N (input) INTEGER +* The order of the matrices A, B, VL, and VR. N >= 0. +* +* A (input/output) REAL array, dimension (LDA, N) +* On entry, the matrix A in the pair (A,B). +* On exit, A has been overwritten. +* +* LDA (input) INTEGER +* The leading dimension of A. LDA >= max(1,N). +* +* B (input/output) REAL array, dimension (LDB, N) +* On entry, the matrix B in the pair (A,B). +* On exit, B has been overwritten. +* +* LDB (input) INTEGER +* The leading dimension of B. LDB >= max(1,N). +* +* ALPHAR (output) REAL array, dimension (N) +* ALPHAI (output) REAL array, dimension (N) +* BETA (output) REAL array, dimension (N) +* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will +* be the generalized eigenvalues. If ALPHAI(j) is zero, then +* the j-th eigenvalue is real; if positive, then the j-th and +* (j+1)-st eigenvalues are a complex conjugate pair, with +* ALPHAI(j+1) negative. +* +* Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) +* may easily over- or underflow, and BETA(j) may even be zero. +* Thus, the user should avoid naively computing the ratio +* alpha/beta. However, ALPHAR and ALPHAI will be always less +* than and usually comparable with norm(A) in magnitude, and +* BETA always less than and usually comparable with norm(B). +* +* VL (output) REAL array, dimension (LDVL,N) +* If JOBVL = 'V', the left eigenvectors u(j) are stored one +* after another in the columns of VL, in the same order as +* their eigenvalues. If the j-th eigenvalue is real, then +* u(j) = VL(:,j), the j-th column of VL. If the j-th and +* (j+1)-th eigenvalues form a complex conjugate pair, then +* u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1). +* Each eigenvector is scaled so the largest component has +* abs(real part)+abs(imag. part)=1. +* Not referenced if JOBVL = 'N'. +* +* LDVL (input) INTEGER +* The leading dimension of the matrix VL. LDVL >= 1, and +* if JOBVL = 'V', LDVL >= N. +* +* VR (output) REAL array, dimension (LDVR,N) +* If JOBVR = 'V', the right eigenvectors v(j) are stored one +* after another in the columns of VR, in the same order as +* their eigenvalues. If the j-th eigenvalue is real, then +* v(j) = VR(:,j), the j-th column of VR. If the j-th and +* (j+1)-th eigenvalues form a complex conjugate pair, then +* v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1). +* Each eigenvector is scaled so the largest component has +* abs(real part)+abs(imag. part)=1. +* Not referenced if JOBVR = 'N'. +* +* LDVR (input) INTEGER +* The leading dimension of the matrix VR. LDVR >= 1, and +* if JOBVR = 'V', LDVR >= N. +* +* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,8*N). +* For good performance, LWORK must generally be larger. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value. +* = 1,...,N: +* The QZ iteration failed. No eigenvectors have been +* calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) +* should be correct for j=INFO+1,...,N. +* > N: =N+1: other than QZ iteration failed in SHGEQZ. +* =N+2: error return from STGEVC. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY + CHARACTER CHTEMP + INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO, + $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK, + $ MINWRK + REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, + $ SMLNUM, TEMP +* .. +* .. Local Arrays .. + LOGICAL LDUMMA( 1 ) +* .. +* .. External Subroutines .. + EXTERNAL SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLABAD, + $ SLACPY, SLASCL, SLASET, SORGQR, SORMQR, STGEVC, + $ XERBLA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + REAL SLAMCH, SLANGE + EXTERNAL LSAME, ILAENV, SLAMCH, SLANGE +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SQRT +* .. +* .. Executable Statements .. +* +* Decode the input arguments +* + IF( LSAME( JOBVL, 'N' ) ) THEN + IJOBVL = 1 + ILVL = .FALSE. + ELSE IF( LSAME( JOBVL, 'V' ) ) THEN + IJOBVL = 2 + ILVL = .TRUE. + ELSE + IJOBVL = -1 + ILVL = .FALSE. + END IF +* + IF( LSAME( JOBVR, 'N' ) ) THEN + IJOBVR = 1 + ILVR = .FALSE. + ELSE IF( LSAME( JOBVR, 'V' ) ) THEN + IJOBVR = 2 + ILVR = .TRUE. + ELSE + IJOBVR = -1 + ILVR = .FALSE. + END IF + ILV = ILVL .OR. ILVR +* +* Test the input arguments +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( IJOBVL.LE.0 ) THEN + INFO = -1 + ELSE IF( IJOBVR.LE.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN + INFO = -12 + ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN + INFO = -14 + END IF +* +* Compute workspace +* (Note: Comments in the code beginning "Workspace:" describe the +* minimal amount of workspace needed at that point in the code, +* as well as the preferred amount for good performance. +* NB refers to the optimal block size for the immediately +* following subroutine, as returned by ILAENV. The workspace is +* computed assuming ILO = 1 and IHI = N, the worst case.) +* + IF( INFO.EQ.0 ) THEN + MINWRK = MAX( 1, 8*N ) + MAXWRK = MAX( 1, N*( 7 + + $ ILAENV( 1, 'SGEQRF', ' ', N, 1, N, 0 ) ) ) + MAXWRK = MAX( MAXWRK, N*( 7 + + $ ILAENV( 1, 'SORMQR', ' ', N, 1, N, 0 ) ) ) + IF( ILVL ) THEN + MAXWRK = MAX( MAXWRK, N*( 7 + + $ ILAENV( 1, 'SORGQR', ' ', N, 1, N, -1 ) ) ) + END IF + WORK( 1 ) = MAXWRK +* + IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) + $ INFO = -16 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SGGEV ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Get machine constants +* + EPS = SLAMCH( 'P' ) + SMLNUM = SLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM + CALL SLABAD( SMLNUM, BIGNUM ) + SMLNUM = SQRT( SMLNUM ) / EPS + BIGNUM = ONE / SMLNUM +* +* Scale A if max element outside range [SMLNUM,BIGNUM] +* + ANRM = SLANGE( 'M', N, N, A, LDA, WORK ) + ILASCL = .FALSE. + IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN + ANRMTO = SMLNUM + ILASCL = .TRUE. + ELSE IF( ANRM.GT.BIGNUM ) THEN + ANRMTO = BIGNUM + ILASCL = .TRUE. + END IF + IF( ILASCL ) + $ CALL SLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) +* +* Scale B if max element outside range [SMLNUM,BIGNUM] +* + BNRM = SLANGE( 'M', N, N, B, LDB, WORK ) + ILBSCL = .FALSE. + IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN + BNRMTO = SMLNUM + ILBSCL = .TRUE. + ELSE IF( BNRM.GT.BIGNUM ) THEN + BNRMTO = BIGNUM + ILBSCL = .TRUE. + END IF + IF( ILBSCL ) + $ CALL SLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) +* +* Permute the matrices A, B to isolate eigenvalues if possible +* (Workspace: need 6*N) +* + ILEFT = 1 + IRIGHT = N + 1 + IWRK = IRIGHT + N + CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ), + $ WORK( IRIGHT ), WORK( IWRK ), IERR ) +* +* Reduce B to triangular form (QR decomposition of B) +* (Workspace: need N, prefer N*NB) +* + IROWS = IHI + 1 - ILO + IF( ILV ) THEN + ICOLS = N + 1 - ILO + ELSE + ICOLS = IROWS + END IF + ITAU = IWRK + IWRK = ITAU + IROWS + CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) +* +* Apply the orthogonal transformation to matrix A +* (Workspace: need N, prefer N*NB) +* + CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, + $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), + $ LWORK+1-IWRK, IERR ) +* +* Initialize VL +* (Workspace: need N, prefer N*NB) +* + IF( ILVL ) THEN + CALL SLASET( 'Full', N, N, ZERO, ONE, VL, LDVL ) + IF( IROWS.GT.1 ) THEN + CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, + $ VL( ILO+1, ILO ), LDVL ) + END IF + CALL SORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, + $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) + END IF +* +* Initialize VR +* + IF( ILVR ) + $ CALL SLASET( 'Full', N, N, ZERO, ONE, VR, LDVR ) +* +* Reduce to generalized Hessenberg form +* (Workspace: none needed) +* + IF( ILV ) THEN +* +* Eigenvectors requested -- work on whole matrix. +* + CALL SGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, + $ LDVL, VR, LDVR, IERR ) + ELSE + CALL SGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, + $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR ) + END IF +* +* Perform QZ algorithm (Compute eigenvalues, and optionally, the +* Schur forms and Schur vectors) +* (Workspace: need N) +* + IWRK = ITAU + IF( ILV ) THEN + CHTEMP = 'S' + ELSE + CHTEMP = 'E' + END IF + CALL SHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, + $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) + IF( IERR.NE.0 ) THEN + IF( IERR.GT.0 .AND. IERR.LE.N ) THEN + INFO = IERR + ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN + INFO = IERR - N + ELSE + INFO = N + 1 + END IF + GO TO 110 + END IF +* +* Compute Eigenvectors +* (Workspace: need 6*N) +* + IF( ILV ) THEN + IF( ILVL ) THEN + IF( ILVR ) THEN + CHTEMP = 'B' + ELSE + CHTEMP = 'L' + END IF + ELSE + CHTEMP = 'R' + END IF + CALL STGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, + $ VR, LDVR, N, IN, WORK( IWRK ), IERR ) + IF( IERR.NE.0 ) THEN + INFO = N + 2 + GO TO 110 + END IF +* +* Undo balancing on VL and VR and normalization +* (Workspace: none needed) +* + IF( ILVL ) THEN + CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ), + $ WORK( IRIGHT ), N, VL, LDVL, IERR ) + DO 50 JC = 1, N + IF( ALPHAI( JC ).LT.ZERO ) + $ GO TO 50 + TEMP = ZERO + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 10 JR = 1, N + TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) ) + 10 CONTINUE + ELSE + DO 20 JR = 1, N + TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+ + $ ABS( VL( JR, JC+1 ) ) ) + 20 CONTINUE + END IF + IF( TEMP.LT.SMLNUM ) + $ GO TO 50 + TEMP = ONE / TEMP + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 30 JR = 1, N + VL( JR, JC ) = VL( JR, JC )*TEMP + 30 CONTINUE + ELSE + DO 40 JR = 1, N + VL( JR, JC ) = VL( JR, JC )*TEMP + VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP + 40 CONTINUE + END IF + 50 CONTINUE + END IF + IF( ILVR ) THEN + CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ), + $ WORK( IRIGHT ), N, VR, LDVR, IERR ) + DO 100 JC = 1, N + IF( ALPHAI( JC ).LT.ZERO ) + $ GO TO 100 + TEMP = ZERO + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 60 JR = 1, N + TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) ) + 60 CONTINUE + ELSE + DO 70 JR = 1, N + TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+ + $ ABS( VR( JR, JC+1 ) ) ) + 70 CONTINUE + END IF + IF( TEMP.LT.SMLNUM ) + $ GO TO 100 + TEMP = ONE / TEMP + IF( ALPHAI( JC ).EQ.ZERO ) THEN + DO 80 JR = 1, N + VR( JR, JC ) = VR( JR, JC )*TEMP + 80 CONTINUE + ELSE + DO 90 JR = 1, N + VR( JR, JC ) = VR( JR, JC )*TEMP + VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP + 90 CONTINUE + END IF + 100 CONTINUE + END IF +* +* End of eigenvector calculation +* + END IF +* +* Undo scaling if necessary +* + IF( ILASCL ) THEN + CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR ) + CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR ) + END IF +* + IF( ILBSCL ) THEN + CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) + END IF +* + 110 CONTINUE +* + WORK( 1 ) = MAXWRK +* + RETURN +* +* End of SGGEV +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zggbak.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,220 @@ + SUBROUTINE ZGGBAK( JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, + $ LDV, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOB, SIDE + INTEGER IHI, ILO, INFO, LDV, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION LSCALE( * ), RSCALE( * ) + COMPLEX*16 V( LDV, * ) +* .. +* +* Purpose +* ======= +* +* ZGGBAK forms the right or left eigenvectors of a complex generalized +* eigenvalue problem A*x = lambda*B*x, by backward transformation on +* the computed eigenvectors of the balanced pair of matrices output by +* ZGGBAL. +* +* Arguments +* ========= +* +* JOB (input) CHARACTER*1 +* Specifies the type of backward transformation required: +* = 'N': do nothing, return immediately; +* = 'P': do backward transformation for permutation only; +* = 'S': do backward transformation for scaling only; +* = 'B': do backward transformations for both permutation and +* scaling. +* JOB must be the same as the argument JOB supplied to ZGGBAL. +* +* SIDE (input) CHARACTER*1 +* = 'R': V contains right eigenvectors; +* = 'L': V contains left eigenvectors. +* +* N (input) INTEGER +* The number of rows of the matrix V. N >= 0. +* +* ILO (input) INTEGER +* IHI (input) INTEGER +* The integers ILO and IHI determined by ZGGBAL. +* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. +* +* LSCALE (input) DOUBLE PRECISION array, dimension (N) +* Details of the permutations and/or scaling factors applied +* to the left side of A and B, as returned by ZGGBAL. +* +* RSCALE (input) DOUBLE PRECISION array, dimension (N) +* Details of the permutations and/or scaling factors applied +* to the right side of A and B, as returned by ZGGBAL. +* +* M (input) INTEGER +* The number of columns of the matrix V. M >= 0. +* +* V (input/output) COMPLEX*16 array, dimension (LDV,M) +* On entry, the matrix of right or left eigenvectors to be +* transformed, as returned by ZTGEVC. +* On exit, V is overwritten by the transformed eigenvectors. +* +* LDV (input) INTEGER +* The leading dimension of the matrix V. LDV >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* Further Details +* =============== +* +* See R.C. Ward, Balancing the generalized eigenvalue problem, +* SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LEFTV, RIGHTV + INTEGER I, K +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZDSCAL, ZSWAP +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + RIGHTV = LSAME( SIDE, 'R' ) + LEFTV = LSAME( SIDE, 'L' ) +* + INFO = 0 + IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. + $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN + INFO = -1 + ELSE IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( ILO.LT.1 ) THEN + INFO = -4 + ELSE IF( N.EQ.0 .AND. IHI.EQ.0 .AND. ILO.NE.1 ) THEN + INFO = -4 + ELSE IF( N.GT.0 .AND. ( IHI.LT.ILO .OR. IHI.GT.MAX( 1, N ) ) ) + $ THEN + INFO = -5 + ELSE IF( N.EQ.0 .AND. ILO.EQ.1 .AND. IHI.NE.0 ) THEN + INFO = -5 + ELSE IF( M.LT.0 ) THEN + INFO = -8 + ELSE IF( LDV.LT.MAX( 1, N ) ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGGBAK', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN + IF( M.EQ.0 ) + $ RETURN + IF( LSAME( JOB, 'N' ) ) + $ RETURN +* + IF( ILO.EQ.IHI ) + $ GO TO 30 +* +* Backward balance +* + IF( LSAME( JOB, 'S' ) .OR. LSAME( JOB, 'B' ) ) THEN +* +* Backward transformation on right eigenvectors +* + IF( RIGHTV ) THEN + DO 10 I = ILO, IHI + CALL ZDSCAL( M, RSCALE( I ), V( I, 1 ), LDV ) + 10 CONTINUE + END IF +* +* Backward transformation on left eigenvectors +* + IF( LEFTV ) THEN + DO 20 I = ILO, IHI + CALL ZDSCAL( M, LSCALE( I ), V( I, 1 ), LDV ) + 20 CONTINUE + END IF + END IF +* +* Backward permutation +* + 30 CONTINUE + IF( LSAME( JOB, 'P' ) .OR. LSAME( JOB, 'B' ) ) THEN +* +* Backward permutation on right eigenvectors +* + IF( RIGHTV ) THEN + IF( ILO.EQ.1 ) + $ GO TO 50 + DO 40 I = ILO - 1, 1, -1 + K = RSCALE( I ) + IF( K.EQ.I ) + $ GO TO 40 + CALL ZSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 40 CONTINUE +* + 50 CONTINUE + IF( IHI.EQ.N ) + $ GO TO 70 + DO 60 I = IHI + 1, N + K = RSCALE( I ) + IF( K.EQ.I ) + $ GO TO 60 + CALL ZSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 60 CONTINUE + END IF +* +* Backward permutation on left eigenvectors +* + 70 CONTINUE + IF( LEFTV ) THEN + IF( ILO.EQ.1 ) + $ GO TO 90 + DO 80 I = ILO - 1, 1, -1 + K = LSCALE( I ) + IF( K.EQ.I ) + $ GO TO 80 + CALL ZSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 80 CONTINUE +* + 90 CONTINUE + IF( IHI.EQ.N ) + $ GO TO 110 + DO 100 I = IHI + 1, N + K = LSCALE( I ) + IF( K.EQ.I ) + $ GO TO 100 + CALL ZSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV ) + 100 CONTINUE + END IF + END IF +* + 110 CONTINUE +* + RETURN +* +* End of ZGGBAK +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zggev.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,454 @@ + SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, + $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO ) +* +* -- LAPACK driver routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOBVL, JOBVR + INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), + $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), + $ WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices +* (A,B), the generalized eigenvalues, and optionally, the left and/or +* right generalized eigenvectors. +* +* A generalized eigenvalue for a pair of matrices (A,B) is a scalar +* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is +* singular. It is usually represented as the pair (alpha,beta), as +* there is a reasonable interpretation for beta=0, and even for both +* being zero. +* +* The right generalized eigenvector v(j) corresponding to the +* generalized eigenvalue lambda(j) of (A,B) satisfies +* +* A * v(j) = lambda(j) * B * v(j). +* +* The left generalized eigenvector u(j) corresponding to the +* generalized eigenvalues lambda(j) of (A,B) satisfies +* +* u(j)**H * A = lambda(j) * u(j)**H * B +* +* where u(j)**H is the conjugate-transpose of u(j). +* +* Arguments +* ========= +* +* JOBVL (input) CHARACTER*1 +* = 'N': do not compute the left generalized eigenvectors; +* = 'V': compute the left generalized eigenvectors. +* +* JOBVR (input) CHARACTER*1 +* = 'N': do not compute the right generalized eigenvectors; +* = 'V': compute the right generalized eigenvectors. +* +* N (input) INTEGER +* The order of the matrices A, B, VL, and VR. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA, N) +* On entry, the matrix A in the pair (A,B). +* On exit, A has been overwritten. +* +* LDA (input) INTEGER +* The leading dimension of A. LDA >= max(1,N). +* +* B (input/output) COMPLEX*16 array, dimension (LDB, N) +* On entry, the matrix B in the pair (A,B). +* On exit, B has been overwritten. +* +* LDB (input) INTEGER +* The leading dimension of B. LDB >= max(1,N). +* +* ALPHA (output) COMPLEX*16 array, dimension (N) +* BETA (output) COMPLEX*16 array, dimension (N) +* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the +* generalized eigenvalues. +* +* Note: the quotients ALPHA(j)/BETA(j) may easily over- or +* underflow, and BETA(j) may even be zero. Thus, the user +* should avoid naively computing the ratio alpha/beta. +* However, ALPHA will be always less than and usually +* comparable with norm(A) in magnitude, and BETA always less +* than and usually comparable with norm(B). +* +* VL (output) COMPLEX*16 array, dimension (LDVL,N) +* If JOBVL = 'V', the left generalized eigenvectors u(j) are +* stored one after another in the columns of VL, in the same +* order as their eigenvalues. +* Each eigenvector is scaled so the largest component has +* abs(real part) + abs(imag. part) = 1. +* Not referenced if JOBVL = 'N'. +* +* LDVL (input) INTEGER +* The leading dimension of the matrix VL. LDVL >= 1, and +* if JOBVL = 'V', LDVL >= N. +* +* VR (output) COMPLEX*16 array, dimension (LDVR,N) +* If JOBVR = 'V', the right generalized eigenvectors v(j) are +* stored one after another in the columns of VR, in the same +* order as their eigenvalues. +* Each eigenvector is scaled so the largest component has +* abs(real part) + abs(imag. part) = 1. +* Not referenced if JOBVR = 'N'. +* +* LDVR (input) INTEGER +* The leading dimension of the matrix VR. LDVR >= 1, and +* if JOBVR = 'V', LDVR >= N. +* +* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,2*N). +* For good performance, LWORK must generally be larger. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* RWORK (workspace/output) DOUBLE PRECISION array, dimension (8*N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value. +* =1,...,N: +* The QZ iteration failed. No eigenvectors have been +* calculated, but ALPHA(j) and BETA(j) should be +* correct for j=INFO+1,...,N. +* > N: =N+1: other then QZ iteration failed in DHGEQZ, +* =N+2: error return from DTGEVC. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), + $ CONE = ( 1.0D0, 0.0D0 ) ) +* .. +* .. Local Scalars .. + LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY + CHARACTER CHTEMP + INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO, + $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR, + $ LWKMIN, LWKOPT + DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, + $ SMLNUM, TEMP + COMPLEX*16 X +* .. +* .. Local Arrays .. + LOGICAL LDUMMA( 1 ) +* .. +* .. External Subroutines .. + EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, + $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR, + $ ZUNMQR +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + DOUBLE PRECISION DLAMCH, ZLANGE + EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT +* .. +* .. Statement Functions .. + DOUBLE PRECISION ABS1 +* .. +* .. Statement Function definitions .. + ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) +* .. +* .. Executable Statements .. +* +* Decode the input arguments +* + IF( LSAME( JOBVL, 'N' ) ) THEN + IJOBVL = 1 + ILVL = .FALSE. + ELSE IF( LSAME( JOBVL, 'V' ) ) THEN + IJOBVL = 2 + ILVL = .TRUE. + ELSE + IJOBVL = -1 + ILVL = .FALSE. + END IF +* + IF( LSAME( JOBVR, 'N' ) ) THEN + IJOBVR = 1 + ILVR = .FALSE. + ELSE IF( LSAME( JOBVR, 'V' ) ) THEN + IJOBVR = 2 + ILVR = .TRUE. + ELSE + IJOBVR = -1 + ILVR = .FALSE. + END IF + ILV = ILVL .OR. ILVR +* +* Test the input arguments +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( IJOBVL.LE.0 ) THEN + INFO = -1 + ELSE IF( IJOBVR.LE.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN + INFO = -11 + ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN + INFO = -13 + END IF +* +* Compute workspace +* (Note: Comments in the code beginning "Workspace:" describe the +* minimal amount of workspace needed at that point in the code, +* as well as the preferred amount for good performance. +* NB refers to the optimal block size for the immediately +* following subroutine, as returned by ILAENV. The workspace is +* computed assuming ILO = 1 and IHI = N, the worst case.) +* + IF( INFO.EQ.0 ) THEN + LWKMIN = MAX( 1, 2*N ) + LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) ) + LWKOPT = MAX( LWKOPT, N + + $ N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) ) + IF( ILVL ) THEN + LWKOPT = MAX( LWKOPT, N + + $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) + END IF + WORK( 1 ) = LWKOPT +* + IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) + $ INFO = -15 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGGEV ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Get machine constants +* + EPS = DLAMCH( 'E' )*DLAMCH( 'B' ) + SMLNUM = DLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) + SMLNUM = SQRT( SMLNUM ) / EPS + BIGNUM = ONE / SMLNUM +* +* Scale A if max element outside range [SMLNUM,BIGNUM] +* + ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK ) + ILASCL = .FALSE. + IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN + ANRMTO = SMLNUM + ILASCL = .TRUE. + ELSE IF( ANRM.GT.BIGNUM ) THEN + ANRMTO = BIGNUM + ILASCL = .TRUE. + END IF + IF( ILASCL ) + $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) +* +* Scale B if max element outside range [SMLNUM,BIGNUM] +* + BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK ) + ILBSCL = .FALSE. + IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN + BNRMTO = SMLNUM + ILBSCL = .TRUE. + ELSE IF( BNRM.GT.BIGNUM ) THEN + BNRMTO = BIGNUM + ILBSCL = .TRUE. + END IF + IF( ILBSCL ) + $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) +* +* Permute the matrices A, B to isolate eigenvalues if possible +* (Real Workspace: need 6*N) +* + ILEFT = 1 + IRIGHT = N + 1 + IRWRK = IRIGHT + N + CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR ) +* +* Reduce B to triangular form (QR decomposition of B) +* (Complex Workspace: need N, prefer N*NB) +* + IROWS = IHI + 1 - ILO + IF( ILV ) THEN + ICOLS = N + 1 - ILO + ELSE + ICOLS = IROWS + END IF + ITAU = 1 + IWRK = ITAU + IROWS + CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) +* +* Apply the orthogonal transformation to matrix A +* (Complex Workspace: need N, prefer N*NB) +* + CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, + $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), + $ LWORK+1-IWRK, IERR ) +* +* Initialize VL +* (Complex Workspace: need N, prefer N*NB) +* + IF( ILVL ) THEN + CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL ) + IF( IROWS.GT.1 ) THEN + CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, + $ VL( ILO+1, ILO ), LDVL ) + END IF + CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, + $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) + END IF +* +* Initialize VR +* + IF( ILVR ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR ) +* +* Reduce to generalized Hessenberg form +* + IF( ILV ) THEN +* +* Eigenvectors requested -- work on whole matrix. +* + CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, + $ LDVL, VR, LDVR, IERR ) + ELSE + CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, + $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR ) + END IF +* +* Perform QZ algorithm (Compute eigenvalues, and optionally, the +* Schur form and Schur vectors) +* (Complex Workspace: need N) +* (Real Workspace: need N) +* + IWRK = ITAU + IF( ILV ) THEN + CHTEMP = 'S' + ELSE + CHTEMP = 'E' + END IF + CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, + $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ), + $ LWORK+1-IWRK, RWORK( IRWRK ), IERR ) + IF( IERR.NE.0 ) THEN + IF( IERR.GT.0 .AND. IERR.LE.N ) THEN + INFO = IERR + ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN + INFO = IERR - N + ELSE + INFO = N + 1 + END IF + GO TO 70 + END IF +* +* Compute Eigenvectors +* (Real Workspace: need 2*N) +* (Complex Workspace: need 2*N) +* + IF( ILV ) THEN + IF( ILVL ) THEN + IF( ILVR ) THEN + CHTEMP = 'B' + ELSE + CHTEMP = 'L' + END IF + ELSE + CHTEMP = 'R' + END IF +* + CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, + $ VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ), + $ IERR ) + IF( IERR.NE.0 ) THEN + INFO = N + 2 + GO TO 70 + END IF +* +* Undo balancing on VL and VR and normalization +* (Workspace: none needed) +* + IF( ILVL ) THEN + CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), N, VL, LDVL, IERR ) + DO 30 JC = 1, N + TEMP = ZERO + DO 10 JR = 1, N + TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) ) + 10 CONTINUE + IF( TEMP.LT.SMLNUM ) + $ GO TO 30 + TEMP = ONE / TEMP + DO 20 JR = 1, N + VL( JR, JC ) = VL( JR, JC )*TEMP + 20 CONTINUE + 30 CONTINUE + END IF + IF( ILVR ) THEN + CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), N, VR, LDVR, IERR ) + DO 60 JC = 1, N + TEMP = ZERO + DO 40 JR = 1, N + TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) ) + 40 CONTINUE + IF( TEMP.LT.SMLNUM ) + $ GO TO 60 + TEMP = ONE / TEMP + DO 50 JR = 1, N + VR( JR, JC ) = VR( JR, JC )*TEMP + 50 CONTINUE + 60 CONTINUE + END IF + END IF +* +* Undo scaling if necessary +* + IF( ILASCL ) + $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) +* + IF( ILBSCL ) + $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) +* + 70 CONTINUE + WORK( 1 ) = LWKOPT +* + RETURN +* +* End of ZGGEV +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zgghrd.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,264 @@ + SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, + $ LDQ, Z, LDZ, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER COMPQ, COMPZ + INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ Z( LDZ, * ) +* .. +* +* Purpose +* ======= +* +* ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper +* Hessenberg form using unitary transformations, where A is a +* general matrix and B is upper triangular. The form of the +* generalized eigenvalue problem is +* A*x = lambda*B*x, +* and B is typically made upper triangular by computing its QR +* factorization and moving the unitary matrix Q to the left side +* of the equation. +* +* This subroutine simultaneously reduces A to a Hessenberg matrix H: +* Q**H*A*Z = H +* and transforms B to another upper triangular matrix T: +* Q**H*B*Z = T +* in order to reduce the problem to its standard form +* H*y = lambda*T*y +* where y = Z**H*x. +* +* The unitary matrices Q and Z are determined as products of Givens +* rotations. They may either be formed explicitly, or they may be +* postmultiplied into input matrices Q1 and Z1, so that +* Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H +* Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H +* If Q1 is the unitary matrix from the QR factorization of B in the +* original equation A*x = lambda*B*x, then ZGGHRD reduces the original +* problem to generalized Hessenberg form. +* +* Arguments +* ========= +* +* COMPQ (input) CHARACTER*1 +* = 'N': do not compute Q; +* = 'I': Q is initialized to the unit matrix, and the +* unitary matrix Q is returned; +* = 'V': Q must contain a unitary matrix Q1 on entry, +* and the product Q1*Q is returned. +* +* COMPZ (input) CHARACTER*1 +* = 'N': do not compute Q; +* = 'I': Q is initialized to the unit matrix, and the +* unitary matrix Q is returned; +* = 'V': Q must contain a unitary matrix Q1 on entry, +* and the product Q1*Q is returned. +* +* N (input) INTEGER +* The order of the matrices A and B. N >= 0. +* +* ILO (input) INTEGER +* IHI (input) INTEGER +* ILO and IHI mark the rows and columns of A which are to be +* reduced. It is assumed that A is already upper triangular +* in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are +* normally set by a previous call to ZGGBAL; otherwise they +* should be set to 1 and N respectively. +* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA, N) +* On entry, the N-by-N general matrix to be reduced. +* On exit, the upper triangle and the first subdiagonal of A +* are overwritten with the upper Hessenberg matrix H, and the +* rest is set to zero. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* B (input/output) COMPLEX*16 array, dimension (LDB, N) +* On entry, the N-by-N upper triangular matrix B. +* On exit, the upper triangular matrix T = Q**H B Z. The +* elements below the diagonal are set to zero. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,N). +* +* Q (input/output) COMPLEX*16 array, dimension (LDQ, N) +* On entry, if COMPQ = 'V', the unitary matrix Q1, typically +* from the QR factorization of B. +* On exit, if COMPQ='I', the unitary matrix Q, and if +* COMPQ = 'V', the product Q1*Q. +* Not referenced if COMPQ='N'. +* +* LDQ (input) INTEGER +* The leading dimension of the array Q. +* LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. +* +* Z (input/output) COMPLEX*16 array, dimension (LDZ, N) +* On entry, if COMPZ = 'V', the unitary matrix Z1. +* On exit, if COMPZ='I', the unitary matrix Z, and if +* COMPZ = 'V', the product Z1*Z. +* Not referenced if COMPZ='N'. +* +* LDZ (input) INTEGER +* The leading dimension of the array Z. +* LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* Further Details +* =============== +* +* This routine reduces A to Hessenberg and B to triangular form by +* an unblocked reduction, as described in _Matrix_Computations_, +* by Golub and van Loan (Johns Hopkins Press). +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 CONE, CZERO + PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ), + $ CZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL ILQ, ILZ + INTEGER ICOMPQ, ICOMPZ, JCOL, JROW + DOUBLE PRECISION C + COMPLEX*16 CTEMP, S +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT +* .. +* .. Intrinsic Functions .. + INTRINSIC DCONJG, MAX +* .. +* .. Executable Statements .. +* +* Decode COMPQ +* + IF( LSAME( COMPQ, 'N' ) ) THEN + ILQ = .FALSE. + ICOMPQ = 1 + ELSE IF( LSAME( COMPQ, 'V' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 2 + ELSE IF( LSAME( COMPQ, 'I' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 3 + ELSE + ICOMPQ = 0 + END IF +* +* Decode COMPZ +* + IF( LSAME( COMPZ, 'N' ) ) THEN + ILZ = .FALSE. + ICOMPZ = 1 + ELSE IF( LSAME( COMPZ, 'V' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 2 + ELSE IF( LSAME( COMPZ, 'I' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 3 + ELSE + ICOMPZ = 0 + END IF +* +* Test the input parameters. +* + INFO = 0 + IF( ICOMPQ.LE.0 ) THEN + INFO = -1 + ELSE IF( ICOMPZ.LE.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( ILO.LT.1 ) THEN + INFO = -4 + ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -9 + ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN + INFO = -11 + ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN + INFO = -13 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGGHRD', -INFO ) + RETURN + END IF +* +* Initialize Q and Z if desired. +* + IF( ICOMPQ.EQ.3 ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) + IF( ICOMPZ.EQ.3 ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) +* +* Quick return if possible +* + IF( N.LE.1 ) + $ RETURN +* +* Zero out lower triangle of B +* + DO 20 JCOL = 1, N - 1 + DO 10 JROW = JCOL + 1, N + B( JROW, JCOL ) = CZERO + 10 CONTINUE + 20 CONTINUE +* +* Reduce A and B +* + DO 40 JCOL = ILO, IHI - 2 +* + DO 30 JROW = IHI, JCOL + 2, -1 +* +* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) +* + CTEMP = A( JROW-1, JCOL ) + CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S, + $ A( JROW-1, JCOL ) ) + A( JROW, JCOL ) = CZERO + CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA, + $ A( JROW, JCOL+1 ), LDA, C, S ) + CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB, + $ B( JROW, JROW-1 ), LDB, C, S ) + IF( ILQ ) + $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, + $ DCONJG( S ) ) +* +* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) +* + CTEMP = B( JROW, JROW ) + CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S, + $ B( JROW, JROW ) ) + B( JROW, JROW-1 ) = CZERO + CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S ) + CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C, + $ S ) + IF( ILZ ) + $ CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S ) + 30 CONTINUE + 40 CONTINUE +* + RETURN +* +* End of ZGGHRD +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zhgeqz.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,759 @@ + SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, + $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, + $ RWORK, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER COMPQ, COMPZ, JOB + INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ), + $ Q( LDQ, * ), T( LDT, * ), WORK( * ), + $ Z( LDZ, * ) +* .. +* +* Purpose +* ======= +* +* ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T), +* where H is an upper Hessenberg matrix and T is upper triangular, +* using the single-shift QZ method. +* Matrix pairs of this type are produced by the reduction to +* generalized upper Hessenberg form of a complex matrix pair (A,B): +* +* A = Q1*H*Z1**H, B = Q1*T*Z1**H, +* +* as computed by ZGGHRD. +* +* If JOB='S', then the Hessenberg-triangular pair (H,T) is +* also reduced to generalized Schur form, +* +* H = Q*S*Z**H, T = Q*P*Z**H, +* +* where Q and Z are unitary matrices and S and P are upper triangular. +* +* Optionally, the unitary matrix Q from the generalized Schur +* factorization may be postmultiplied into an input matrix Q1, and the +* unitary matrix Z may be postmultiplied into an input matrix Z1. +* If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced +* the matrix pair (A,B) to generalized Hessenberg form, then the output +* matrices Q1*Q and Z1*Z are the unitary factors from the generalized +* Schur factorization of (A,B): +* +* A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. +* +* To avoid overflow, eigenvalues of the matrix pair (H,T) +* (equivalently, of (A,B)) are computed as a pair of complex values +* (alpha,beta). If beta is nonzero, lambda = alpha / beta is an +* eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP) +* A*x = lambda*B*x +* and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the +* alternate form of the GNEP +* mu*A*y = B*y. +* The values of alpha and beta for the i-th eigenvalue can be read +* directly from the generalized Schur form: alpha = S(i,i), +* beta = P(i,i). +* +* Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix +* Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), +* pp. 241--256. +* +* Arguments +* ========= +* +* JOB (input) CHARACTER*1 +* = 'E': Compute eigenvalues only; +* = 'S': Computer eigenvalues and the Schur form. +* +* COMPQ (input) CHARACTER*1 +* = 'N': Left Schur vectors (Q) are not computed; +* = 'I': Q is initialized to the unit matrix and the matrix Q +* of left Schur vectors of (H,T) is returned; +* = 'V': Q must contain a unitary matrix Q1 on entry and +* the product Q1*Q is returned. +* +* COMPZ (input) CHARACTER*1 +* = 'N': Right Schur vectors (Z) are not computed; +* = 'I': Q is initialized to the unit matrix and the matrix Z +* of right Schur vectors of (H,T) is returned; +* = 'V': Z must contain a unitary matrix Z1 on entry and +* the product Z1*Z is returned. +* +* N (input) INTEGER +* The order of the matrices H, T, Q, and Z. N >= 0. +* +* ILO (input) INTEGER +* IHI (input) INTEGER +* ILO and IHI mark the rows and columns of H which are in +* Hessenberg form. It is assumed that A is already upper +* triangular in rows and columns 1:ILO-1 and IHI+1:N. +* If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. +* +* H (input/output) COMPLEX*16 array, dimension (LDH, N) +* On entry, the N-by-N upper Hessenberg matrix H. +* On exit, if JOB = 'S', H contains the upper triangular +* matrix S from the generalized Schur factorization. +* If JOB = 'E', the diagonal of H matches that of S, but +* the rest of H is unspecified. +* +* LDH (input) INTEGER +* The leading dimension of the array H. LDH >= max( 1, N ). +* +* T (input/output) COMPLEX*16 array, dimension (LDT, N) +* On entry, the N-by-N upper triangular matrix T. +* On exit, if JOB = 'S', T contains the upper triangular +* matrix P from the generalized Schur factorization. +* If JOB = 'E', the diagonal of T matches that of P, but +* the rest of T is unspecified. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= max( 1, N ). +* +* ALPHA (output) COMPLEX*16 array, dimension (N) +* The complex scalars alpha that define the eigenvalues of +* GNEP. ALPHA(i) = S(i,i) in the generalized Schur +* factorization. +* +* BETA (output) COMPLEX*16 array, dimension (N) +* The real non-negative scalars beta that define the +* eigenvalues of GNEP. BETA(i) = P(i,i) in the generalized +* Schur factorization. +* +* Together, the quantities alpha = ALPHA(j) and beta = BETA(j) +* represent the j-th eigenvalue of the matrix pair (A,B), in +* one of the forms lambda = alpha/beta or mu = beta/alpha. +* Since either lambda or mu may overflow, they should not, +* in general, be computed. +* +* Q (input/output) COMPLEX*16 array, dimension (LDQ, N) +* On entry, if COMPZ = 'V', the unitary matrix Q1 used in the +* reduction of (A,B) to generalized Hessenberg form. +* On exit, if COMPZ = 'I', the unitary matrix of left Schur +* vectors of (H,T), and if COMPZ = 'V', the unitary matrix of +* left Schur vectors of (A,B). +* Not referenced if COMPZ = 'N'. +* +* LDQ (input) INTEGER +* The leading dimension of the array Q. LDQ >= 1. +* If COMPQ='V' or 'I', then LDQ >= N. +* +* Z (input/output) COMPLEX*16 array, dimension (LDZ, N) +* On entry, if COMPZ = 'V', the unitary matrix Z1 used in the +* reduction of (A,B) to generalized Hessenberg form. +* On exit, if COMPZ = 'I', the unitary matrix of right Schur +* vectors of (H,T), and if COMPZ = 'V', the unitary matrix of +* right Schur vectors of (A,B). +* Not referenced if COMPZ = 'N'. +* +* LDZ (input) INTEGER +* The leading dimension of the array Z. LDZ >= 1. +* If COMPZ='V' or 'I', then LDZ >= N. +* +* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) +* On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* RWORK (workspace) DOUBLE PRECISION array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* = 1,...,N: the QZ iteration did not converge. (H,T) is not +* in Schur form, but ALPHA(i) and BETA(i), +* i=INFO+1,...,N should be correct. +* = N+1,...,2*N: the shift calculation failed. (H,T) is not +* in Schur form, but ALPHA(i) and BETA(i), +* i=INFO-N+1,...,N should be correct. +* +* Further Details +* =============== +* +* We assume that complex ABS works as long as its value is less than +* overflow. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), + $ CONE = ( 1.0D+0, 0.0D+0 ) ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + DOUBLE PRECISION HALF + PARAMETER ( HALF = 0.5D+0 ) +* .. +* .. Local Scalars .. + LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY + INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, + $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, + $ JR, MAXIT + DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, + $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP + COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, + $ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1, + $ U12, X +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, ZLANHS + EXTERNAL LSAME, DLAMCH, ZLANHS +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN, + $ SQRT +* .. +* .. Statement Functions .. + DOUBLE PRECISION ABS1 +* .. +* .. Statement Function definitions .. + ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) +* .. +* .. Executable Statements .. +* +* Decode JOB, COMPQ, COMPZ +* + IF( LSAME( JOB, 'E' ) ) THEN + ILSCHR = .FALSE. + ISCHUR = 1 + ELSE IF( LSAME( JOB, 'S' ) ) THEN + ILSCHR = .TRUE. + ISCHUR = 2 + ELSE + ISCHUR = 0 + END IF +* + IF( LSAME( COMPQ, 'N' ) ) THEN + ILQ = .FALSE. + ICOMPQ = 1 + ELSE IF( LSAME( COMPQ, 'V' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 2 + ELSE IF( LSAME( COMPQ, 'I' ) ) THEN + ILQ = .TRUE. + ICOMPQ = 3 + ELSE + ICOMPQ = 0 + END IF +* + IF( LSAME( COMPZ, 'N' ) ) THEN + ILZ = .FALSE. + ICOMPZ = 1 + ELSE IF( LSAME( COMPZ, 'V' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 2 + ELSE IF( LSAME( COMPZ, 'I' ) ) THEN + ILZ = .TRUE. + ICOMPZ = 3 + ELSE + ICOMPZ = 0 + END IF +* +* Check Argument Values +* + INFO = 0 + WORK( 1 ) = MAX( 1, N ) + LQUERY = ( LWORK.EQ.-1 ) + IF( ISCHUR.EQ.0 ) THEN + INFO = -1 + ELSE IF( ICOMPQ.EQ.0 ) THEN + INFO = -2 + ELSE IF( ICOMPZ.EQ.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( ILO.LT.1 ) THEN + INFO = -5 + ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN + INFO = -6 + ELSE IF( LDH.LT.N ) THEN + INFO = -8 + ELSE IF( LDT.LT.N ) THEN + INFO = -10 + ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN + INFO = -14 + ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN + INFO = -16 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -18 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZHGEQZ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* +* WORK( 1 ) = CMPLX( 1 ) + IF( N.LE.0 ) THEN + WORK( 1 ) = DCMPLX( 1 ) + RETURN + END IF +* +* Initialize Q and Z +* + IF( ICOMPQ.EQ.3 ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) + IF( ICOMPZ.EQ.3 ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) +* +* Machine Constants +* + IN = IHI + 1 - ILO + SAFMIN = DLAMCH( 'S' ) + ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) + ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK ) + BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK ) + ATOL = MAX( SAFMIN, ULP*ANORM ) + BTOL = MAX( SAFMIN, ULP*BNORM ) + ASCALE = ONE / MAX( SAFMIN, ANORM ) + BSCALE = ONE / MAX( SAFMIN, BNORM ) +* +* +* Set Eigenvalues IHI+1:N +* + DO 10 J = IHI + 1, N + ABSB = ABS( T( J, J ) ) + IF( ABSB.GT.SAFMIN ) THEN + SIGNBC = DCONJG( T( J, J ) / ABSB ) + T( J, J ) = ABSB + IF( ILSCHR ) THEN + CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) + CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) + ELSE + H( J, J ) = H( J, J )*SIGNBC + END IF + IF( ILZ ) + $ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) + ELSE + T( J, J ) = CZERO + END IF + ALPHA( J ) = H( J, J ) + BETA( J ) = T( J, J ) + 10 CONTINUE +* +* If IHI < ILO, skip QZ steps +* + IF( IHI.LT.ILO ) + $ GO TO 190 +* +* MAIN QZ ITERATION LOOP +* +* Initialize dynamic indices +* +* Eigenvalues ILAST+1:N have been found. +* Column operations modify rows IFRSTM:whatever +* Row operations modify columns whatever:ILASTM +* +* If only eigenvalues are being computed, then +* IFRSTM is the row of the last splitting row above row ILAST; +* this is always at least ILO. +* IITER counts iterations since the last eigenvalue was found, +* to tell when to use an extraordinary shift. +* MAXIT is the maximum number of QZ sweeps allowed. +* + ILAST = IHI + IF( ILSCHR ) THEN + IFRSTM = 1 + ILASTM = N + ELSE + IFRSTM = ILO + ILASTM = IHI + END IF + IITER = 0 + ESHIFT = CZERO + MAXIT = 30*( IHI-ILO+1 ) +* + DO 170 JITER = 1, MAXIT +* +* Check for too many iterations. +* + IF( JITER.GT.MAXIT ) + $ GO TO 180 +* +* Split the matrix if possible. +* +* Two tests: +* 1: H(j,j-1)=0 or j=ILO +* 2: T(j,j)=0 +* +* Special case: j=ILAST +* + IF( ILAST.EQ.ILO ) THEN + GO TO 60 + ELSE + IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN + H( ILAST, ILAST-1 ) = CZERO + GO TO 60 + END IF + END IF +* + IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN + T( ILAST, ILAST ) = CZERO + GO TO 50 + END IF +* +* General case: j<ILAST +* + DO 40 J = ILAST - 1, ILO, -1 +* +* Test 1: for H(j,j-1)=0 or j=ILO +* + IF( J.EQ.ILO ) THEN + ILAZRO = .TRUE. + ELSE + IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN + H( J, J-1 ) = CZERO + ILAZRO = .TRUE. + ELSE + ILAZRO = .FALSE. + END IF + END IF +* +* Test 2: for T(j,j)=0 +* + IF( ABS( T( J, J ) ).LT.BTOL ) THEN + T( J, J ) = CZERO +* +* Test 1a: Check for 2 consecutive small subdiagonals in A +* + ILAZR2 = .FALSE. + IF( .NOT.ILAZRO ) THEN + IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1, + $ J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) ) + $ ILAZR2 = .TRUE. + END IF +* +* If both tests pass (1 & 2), i.e., the leading diagonal +* element of B in the block is zero, split a 1x1 block off +* at the top. (I.e., at the J-th row/column) The leading +* diagonal element of the remainder can also be zero, so +* this may have to be done repeatedly. +* + IF( ILAZRO .OR. ILAZR2 ) THEN + DO 20 JCH = J, ILAST - 1 + CTEMP = H( JCH, JCH ) + CALL ZLARTG( CTEMP, H( JCH+1, JCH ), C, S, + $ H( JCH, JCH ) ) + H( JCH+1, JCH ) = CZERO + CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, + $ H( JCH+1, JCH+1 ), LDH, C, S ) + CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, + $ T( JCH+1, JCH+1 ), LDT, C, S ) + IF( ILQ ) + $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, + $ C, DCONJG( S ) ) + IF( ILAZR2 ) + $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C + ILAZR2 = .FALSE. + IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN + IF( JCH+1.GE.ILAST ) THEN + GO TO 60 + ELSE + IFIRST = JCH + 1 + GO TO 70 + END IF + END IF + T( JCH+1, JCH+1 ) = CZERO + 20 CONTINUE + GO TO 50 + ELSE +* +* Only test 2 passed -- chase the zero to T(ILAST,ILAST) +* Then process as in the case T(ILAST,ILAST)=0 +* + DO 30 JCH = J, ILAST - 1 + CTEMP = T( JCH, JCH+1 ) + CALL ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S, + $ T( JCH, JCH+1 ) ) + T( JCH+1, JCH+1 ) = CZERO + IF( JCH.LT.ILASTM-1 ) + $ CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, + $ T( JCH+1, JCH+2 ), LDT, C, S ) + CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, + $ H( JCH+1, JCH-1 ), LDH, C, S ) + IF( ILQ ) + $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, + $ C, DCONJG( S ) ) + CTEMP = H( JCH+1, JCH ) + CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S, + $ H( JCH+1, JCH ) ) + H( JCH+1, JCH-1 ) = CZERO + CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, + $ H( IFRSTM, JCH-1 ), 1, C, S ) + CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, + $ T( IFRSTM, JCH-1 ), 1, C, S ) + IF( ILZ ) + $ CALL ZROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, + $ C, S ) + 30 CONTINUE + GO TO 50 + END IF + ELSE IF( ILAZRO ) THEN +* +* Only test 1 passed -- work on J:ILAST +* + IFIRST = J + GO TO 70 + END IF +* +* Neither test passed -- try next J +* + 40 CONTINUE +* +* (Drop-through is "impossible") +* + INFO = 2*N + 1 + GO TO 210 +* +* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a +* 1x1 block. +* + 50 CONTINUE + CTEMP = H( ILAST, ILAST ) + CALL ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S, + $ H( ILAST, ILAST ) ) + H( ILAST, ILAST-1 ) = CZERO + CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, + $ H( IFRSTM, ILAST-1 ), 1, C, S ) + CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, + $ T( IFRSTM, ILAST-1 ), 1, C, S ) + IF( ILZ ) + $ CALL ZROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) +* +* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA +* + 60 CONTINUE + ABSB = ABS( T( ILAST, ILAST ) ) + IF( ABSB.GT.SAFMIN ) THEN + SIGNBC = DCONJG( T( ILAST, ILAST ) / ABSB ) + T( ILAST, ILAST ) = ABSB + IF( ILSCHR ) THEN + CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 ) + CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ), + $ 1 ) + ELSE + H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC + END IF + IF( ILZ ) + $ CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 ) + ELSE + T( ILAST, ILAST ) = CZERO + END IF + ALPHA( ILAST ) = H( ILAST, ILAST ) + BETA( ILAST ) = T( ILAST, ILAST ) +* +* Go to next block -- exit if finished. +* + ILAST = ILAST - 1 + IF( ILAST.LT.ILO ) + $ GO TO 190 +* +* Reset counters +* + IITER = 0 + ESHIFT = CZERO + IF( .NOT.ILSCHR ) THEN + ILASTM = ILAST + IF( IFRSTM.GT.ILAST ) + $ IFRSTM = ILO + END IF + GO TO 160 +* +* QZ step +* +* This iteration only involves rows/columns IFIRST:ILAST. We +* assume IFIRST < ILAST, and that the diagonal of B is non-zero. +* + 70 CONTINUE + IITER = IITER + 1 + IF( .NOT.ILSCHR ) THEN + IFRSTM = IFIRST + END IF +* +* Compute the Shift. +* +* At this point, IFIRST < ILAST, and the diagonal elements of +* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in +* magnitude) +* + IF( ( IITER / 10 )*10.NE.IITER ) THEN +* +* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of +* the bottom-right 2x2 block of A inv(B) which is nearest to +* the bottom-right element. +* +* We factor B as U*D, where U has unit diagonals, and +* compute (A*inv(D))*inv(U). +* + U12 = ( BSCALE*T( ILAST-1, ILAST ) ) / + $ ( BSCALE*T( ILAST, ILAST ) ) + AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / + $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) + AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / + $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) + AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / + $ ( BSCALE*T( ILAST, ILAST ) ) + AD22 = ( ASCALE*H( ILAST, ILAST ) ) / + $ ( BSCALE*T( ILAST, ILAST ) ) + ABI22 = AD22 - U12*AD21 +* + T1 = HALF*( AD11+ABI22 ) + RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 ) + TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) + + $ DIMAG( T1-ABI22 )*DIMAG( RTDISC ) + IF( TEMP.LE.ZERO ) THEN + SHIFT = T1 + RTDISC + ELSE + SHIFT = T1 - RTDISC + END IF + ELSE +* +* Exceptional shift. Chosen for no particularly good reason. +* + ESHIFT = ESHIFT + DCONJG( ( ASCALE*H( ILAST-1, ILAST ) ) / + $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) ) + SHIFT = ESHIFT + END IF +* +* Now check for two consecutive small subdiagonals. +* + DO 80 J = ILAST - 1, IFIRST + 1, -1 + ISTART = J + CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) ) + TEMP = ABS1( CTEMP ) + TEMP2 = ASCALE*ABS1( H( J+1, J ) ) + TEMPR = MAX( TEMP, TEMP2 ) + IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN + TEMP = TEMP / TEMPR + TEMP2 = TEMP2 / TEMPR + END IF + IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL ) + $ GO TO 90 + 80 CONTINUE +* + ISTART = IFIRST + CTEMP = ASCALE*H( IFIRST, IFIRST ) - + $ SHIFT*( BSCALE*T( IFIRST, IFIRST ) ) + 90 CONTINUE +* +* Do an implicit-shift QZ sweep. +* +* Initial Q +* + CTEMP2 = ASCALE*H( ISTART+1, ISTART ) + CALL ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) +* +* Sweep +* + DO 150 J = ISTART, ILAST - 1 + IF( J.GT.ISTART ) THEN + CTEMP = H( J, J-1 ) + CALL ZLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) + H( J+1, J-1 ) = CZERO + END IF +* + DO 100 JC = J, ILASTM + CTEMP = C*H( J, JC ) + S*H( J+1, JC ) + H( J+1, JC ) = -DCONJG( S )*H( J, JC ) + C*H( J+1, JC ) + H( J, JC ) = CTEMP + CTEMP2 = C*T( J, JC ) + S*T( J+1, JC ) + T( J+1, JC ) = -DCONJG( S )*T( J, JC ) + C*T( J+1, JC ) + T( J, JC ) = CTEMP2 + 100 CONTINUE + IF( ILQ ) THEN + DO 110 JR = 1, N + CTEMP = C*Q( JR, J ) + DCONJG( S )*Q( JR, J+1 ) + Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) + Q( JR, J ) = CTEMP + 110 CONTINUE + END IF +* + CTEMP = T( J+1, J+1 ) + CALL ZLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) + T( J+1, J ) = CZERO +* + DO 120 JR = IFRSTM, MIN( J+2, ILAST ) + CTEMP = C*H( JR, J+1 ) + S*H( JR, J ) + H( JR, J ) = -DCONJG( S )*H( JR, J+1 ) + C*H( JR, J ) + H( JR, J+1 ) = CTEMP + 120 CONTINUE + DO 130 JR = IFRSTM, J + CTEMP = C*T( JR, J+1 ) + S*T( JR, J ) + T( JR, J ) = -DCONJG( S )*T( JR, J+1 ) + C*T( JR, J ) + T( JR, J+1 ) = CTEMP + 130 CONTINUE + IF( ILZ ) THEN + DO 140 JR = 1, N + CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) + Z( JR, J ) = -DCONJG( S )*Z( JR, J+1 ) + C*Z( JR, J ) + Z( JR, J+1 ) = CTEMP + 140 CONTINUE + END IF + 150 CONTINUE +* + 160 CONTINUE +* + 170 CONTINUE +* +* Drop-through = non-convergence +* + 180 CONTINUE + INFO = ILAST + GO TO 210 +* +* Successful completion of all QZ steps +* + 190 CONTINUE +* +* Set Eigenvalues 1:ILO-1 +* + DO 200 J = 1, ILO - 1 + ABSB = ABS( T( J, J ) ) + IF( ABSB.GT.SAFMIN ) THEN + SIGNBC = DCONJG( T( J, J ) / ABSB ) + T( J, J ) = ABSB + IF( ILSCHR ) THEN + CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) + CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) + ELSE + H( J, J ) = H( J, J )*SIGNBC + END IF + IF( ILZ ) + $ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) + ELSE + T( J, J ) = CZERO + END IF + ALPHA( J ) = H( J, J ) + BETA( J ) = T( J, J ) + 200 CONTINUE +* +* Normal Termination +* + INFO = 0 +* +* Exit (other than argument error) -- return optimal workspace size +* + 210 CONTINUE + WORK( 1 ) = DCMPLX( N ) + RETURN +* +* End of ZHGEQZ +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/ztgevc.f Mon Dec 15 13:49:16 2008 +0100 @@ -0,0 +1,633 @@ + SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, + $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER HOWMNY, SIDE + INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N +* .. +* .. Array Arguments .. + LOGICAL SELECT( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ), + $ VR( LDVR, * ), WORK( * ) +* .. +* +* +* Purpose +* ======= +* +* ZTGEVC computes some or all of the right and/or left eigenvectors of +* a pair of complex matrices (S,P), where S and P are upper triangular. +* Matrix pairs of this type are produced by the generalized Schur +* factorization of a complex matrix pair (A,B): +* +* A = Q*S*Z**H, B = Q*P*Z**H +* +* as computed by ZGGHRD + ZHGEQZ. +* +* The right eigenvector x and the left eigenvector y of (S,P) +* corresponding to an eigenvalue w are defined by: +* +* S*x = w*P*x, (y**H)*S = w*(y**H)*P, +* +* where y**H denotes the conjugate tranpose of y. +* The eigenvalues are not input to this routine, but are computed +* directly from the diagonal elements of S and P. +* +* This routine returns the matrices X and/or Y of right and left +* eigenvectors of (S,P), or the products Z*X and/or Q*Y, +* where Z and Q are input matrices. +* If Q and Z are the unitary factors from the generalized Schur +* factorization of a matrix pair (A,B), then Z*X and Q*Y +* are the matrices of right and left eigenvectors of (A,B). +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'R': compute right eigenvectors only; +* = 'L': compute left eigenvectors only; +* = 'B': compute both right and left eigenvectors. +* +* HOWMNY (input) CHARACTER*1 +* = 'A': compute all right and/or left eigenvectors; +* = 'B': compute all right and/or left eigenvectors, +* backtransformed by the matrices in VR and/or VL; +* = 'S': compute selected right and/or left eigenvectors, +* specified by the logical array SELECT. +* +* SELECT (input) LOGICAL array, dimension (N) +* If HOWMNY='S', SELECT specifies the eigenvectors to be +* computed. The eigenvector corresponding to the j-th +* eigenvalue is computed if SELECT(j) = .TRUE.. +* Not referenced if HOWMNY = 'A' or 'B'. +* +* N (input) INTEGER +* The order of the matrices S and P. N >= 0. +* +* S (input) COMPLEX*16 array, dimension (LDS,N) +* The upper triangular matrix S from a generalized Schur +* factorization, as computed by ZHGEQZ. +* +* LDS (input) INTEGER +* The leading dimension of array S. LDS >= max(1,N). +* +* P (input) COMPLEX*16 array, dimension (LDP,N) +* The upper triangular matrix P from a generalized Schur +* factorization, as computed by ZHGEQZ. P must have real +* diagonal elements. +* +* LDP (input) INTEGER +* The leading dimension of array P. LDP >= max(1,N). +* +* VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) +* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must +* contain an N-by-N matrix Q (usually the unitary matrix Q +* of left Schur vectors returned by ZHGEQZ). +* On exit, if SIDE = 'L' or 'B', VL contains: +* if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); +* if HOWMNY = 'B', the matrix Q*Y; +* if HOWMNY = 'S', the left eigenvectors of (S,P) specified by +* SELECT, stored consecutively in the columns of +* VL, in the same order as their eigenvalues. +* Not referenced if SIDE = 'R'. +* +* LDVL (input) INTEGER +* The leading dimension of array VL. LDVL >= 1, and if +* SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N. +* +* VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) +* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must +* contain an N-by-N matrix Q (usually the unitary matrix Z +* of right Schur vectors returned by ZHGEQZ). +* On exit, if SIDE = 'R' or 'B', VR contains: +* if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); +* if HOWMNY = 'B', the matrix Z*X; +* if HOWMNY = 'S', the right eigenvectors of (S,P) specified by +* SELECT, stored consecutively in the columns of +* VR, in the same order as their eigenvalues. +* Not referenced if SIDE = 'L'. +* +* LDVR (input) INTEGER +* The leading dimension of the array VR. LDVR >= 1, and if +* SIDE = 'R' or 'B', LDVR >= N. +* +* MM (input) INTEGER +* The number of columns in the arrays VL and/or VR. MM >= M. +* +* M (output) INTEGER +* The number of columns in the arrays VL and/or VR actually +* used to store the eigenvectors. If HOWMNY = 'A' or 'B', M +* is set to N. Each selected eigenvector occupies one column. +* +* WORK (workspace) COMPLEX*16 array, dimension (2*N) +* +* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), + $ CONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP, + $ LSA, LSB + INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC, + $ J, JE, JR + DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG, + $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA, + $ SCALE, SMALL, TEMP, ULP, XMAX + COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH + COMPLEX*16 ZLADIV + EXTERNAL LSAME, DLAMCH, ZLADIV +* .. +* .. External Subroutines .. + EXTERNAL DLABAD, XERBLA, ZGEMV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN +* .. +* .. Statement Functions .. + DOUBLE PRECISION ABS1 +* .. +* .. Statement Function definitions .. + ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) +* .. +* .. Executable Statements .. +* +* Decode and Test the input parameters +* + IF( LSAME( HOWMNY, 'A' ) ) THEN + IHWMNY = 1 + ILALL = .TRUE. + ILBACK = .FALSE. + ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN + IHWMNY = 2 + ILALL = .FALSE. + ILBACK = .FALSE. + ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN + IHWMNY = 3 + ILALL = .TRUE. + ILBACK = .TRUE. + ELSE + IHWMNY = -1 + END IF +* + IF( LSAME( SIDE, 'R' ) ) THEN + ISIDE = 1 + COMPL = .FALSE. + COMPR = .TRUE. + ELSE IF( LSAME( SIDE, 'L' ) ) THEN + ISIDE = 2 + COMPL = .TRUE. + COMPR = .FALSE. + ELSE IF( LSAME( SIDE, 'B' ) ) THEN + ISIDE = 3 + COMPL = .TRUE. + COMPR = .TRUE. + ELSE + ISIDE = -1 + END IF +* + INFO = 0 + IF( ISIDE.LT.0 ) THEN + INFO = -1 + ELSE IF( IHWMNY.LT.0 ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( LDS.LT.MAX( 1, N ) ) THEN + INFO = -6 + ELSE IF( LDP.LT.MAX( 1, N ) ) THEN + INFO = -8 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZTGEVC', -INFO ) + RETURN + END IF +* +* Count the number of eigenvectors +* + IF( .NOT.ILALL ) THEN + IM = 0 + DO 10 J = 1, N + IF( SELECT( J ) ) + $ IM = IM + 1 + 10 CONTINUE + ELSE + IM = N + END IF +* +* Check diagonal of B +* + ILBBAD = .FALSE. + DO 20 J = 1, N + IF( DIMAG( P( J, J ) ).NE.ZERO ) + $ ILBBAD = .TRUE. + 20 CONTINUE +* + IF( ILBBAD ) THEN + INFO = -7 + ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN + INFO = -10 + ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN + INFO = -12 + ELSE IF( MM.LT.IM ) THEN + INFO = -13 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZTGEVC', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + M = IM + IF( N.EQ.0 ) + $ RETURN +* +* Machine Constants +* + SAFMIN = DLAMCH( 'Safe minimum' ) + BIG = ONE / SAFMIN + CALL DLABAD( SAFMIN, BIG ) + ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) + SMALL = SAFMIN*N / ULP + BIG = ONE / SMALL + BIGNUM = ONE / ( SAFMIN*N ) +* +* Compute the 1-norm of each column of the strictly upper triangular +* part of A and B to check for possible overflow in the triangular +* solver. +* + ANORM = ABS1( S( 1, 1 ) ) + BNORM = ABS1( P( 1, 1 ) ) + RWORK( 1 ) = ZERO + RWORK( N+1 ) = ZERO + DO 40 J = 2, N + RWORK( J ) = ZERO + RWORK( N+J ) = ZERO + DO 30 I = 1, J - 1 + RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) ) + RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) ) + 30 CONTINUE + ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) ) + BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) ) + 40 CONTINUE +* + ASCALE = ONE / MAX( ANORM, SAFMIN ) + BSCALE = ONE / MAX( BNORM, SAFMIN ) +* +* Left eigenvectors +* + IF( COMPL ) THEN + IEIG = 0 +* +* Main loop over eigenvalues +* + DO 140 JE = 1, N + IF( ILALL ) THEN + ILCOMP = .TRUE. + ELSE + ILCOMP = SELECT( JE ) + END IF + IF( ILCOMP ) THEN + IEIG = IEIG + 1 +* + IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND. + $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN +* +* Singular matrix pencil -- return unit eigenvector +* + DO 50 JR = 1, N + VL( JR, IEIG ) = CZERO + 50 CONTINUE + VL( IEIG, IEIG ) = CONE + GO TO 140 + END IF +* +* Non-singular eigenvalue: +* Compute coefficients a and b in +* H +* y ( a A - b B ) = 0 +* + TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE, + $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN ) + SALPHA = ( TEMP*S( JE, JE ) )*ASCALE + SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE + ACOEFF = SBETA*ASCALE + BCOEFF = SALPHA*BSCALE +* +* Scale to avoid underflow +* + LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL + LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. + $ SMALL +* + SCALE = ONE + IF( LSA ) + $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) + IF( LSB ) + $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* + $ MIN( BNORM, BIG ) ) + IF( LSA .OR. LSB ) THEN + SCALE = MIN( SCALE, ONE / + $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), + $ ABS1( BCOEFF ) ) ) ) + IF( LSA ) THEN + ACOEFF = ASCALE*( SCALE*SBETA ) + ELSE + ACOEFF = SCALE*ACOEFF + END IF + IF( LSB ) THEN + BCOEFF = BSCALE*( SCALE*SALPHA ) + ELSE + BCOEFF = SCALE*BCOEFF + END IF + END IF +* + ACOEFA = ABS( ACOEFF ) + BCOEFA = ABS1( BCOEFF ) + XMAX = ONE + DO 60 JR = 1, N + WORK( JR ) = CZERO + 60 CONTINUE + WORK( JE ) = CONE + DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) +* +* H +* Triangular solve of (a A - b B) y = 0 +* +* H +* (rowwise in (a A - b B) , or columnwise in a A - b B) +* + DO 100 J = JE + 1, N +* +* Compute +* j-1 +* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) +* k=je +* (Scale if necessary) +* + TEMP = ONE / XMAX + IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM* + $ TEMP ) THEN + DO 70 JR = JE, J - 1 + WORK( JR ) = TEMP*WORK( JR ) + 70 CONTINUE + XMAX = ONE + END IF + SUMA = CZERO + SUMB = CZERO +* + DO 80 JR = JE, J - 1 + SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR ) + SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR ) + 80 CONTINUE + SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB +* +* Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) ) +* +* with scaling and perturbation of the denominator +* + D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) ) + IF( ABS1( D ).LE.DMIN ) + $ D = DCMPLX( DMIN ) +* + IF( ABS1( D ).LT.ONE ) THEN + IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN + TEMP = ONE / ABS1( SUM ) + DO 90 JR = JE, J - 1 + WORK( JR ) = TEMP*WORK( JR ) + 90 CONTINUE + XMAX = TEMP*XMAX + SUM = TEMP*SUM + END IF + END IF + WORK( J ) = ZLADIV( -SUM, D ) + XMAX = MAX( XMAX, ABS1( WORK( J ) ) ) + 100 CONTINUE +* +* Back transform eigenvector if HOWMNY='B'. +* + IF( ILBACK ) THEN + CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL, + $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 ) + ISRC = 2 + IBEG = 1 + ELSE + ISRC = 1 + IBEG = JE + END IF +* +* Copy and scale eigenvector into column of VL +* + XMAX = ZERO + DO 110 JR = IBEG, N + XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) + 110 CONTINUE +* + IF( XMAX.GT.SAFMIN ) THEN + TEMP = ONE / XMAX + DO 120 JR = IBEG, N + VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) + 120 CONTINUE + ELSE + IBEG = N + 1 + END IF +* + DO 130 JR = 1, IBEG - 1 + VL( JR, IEIG ) = CZERO + 130 CONTINUE +* + END IF + 140 CONTINUE + END IF +* +* Right eigenvectors +* + IF( COMPR ) THEN + IEIG = IM + 1 +* +* Main loop over eigenvalues +* + DO 250 JE = N, 1, -1 + IF( ILALL ) THEN + ILCOMP = .TRUE. + ELSE + ILCOMP = SELECT( JE ) + END IF + IF( ILCOMP ) THEN + IEIG = IEIG - 1 +* + IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND. + $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN +* +* Singular matrix pencil -- return unit eigenvector +* + DO 150 JR = 1, N + VR( JR, IEIG ) = CZERO + 150 CONTINUE + VR( IEIG, IEIG ) = CONE + GO TO 250 + END IF +* +* Non-singular eigenvalue: +* Compute coefficients a and b in +* +* ( a A - b B ) x = 0 +* + TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE, + $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN ) + SALPHA = ( TEMP*S( JE, JE ) )*ASCALE + SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE + ACOEFF = SBETA*ASCALE + BCOEFF = SALPHA*BSCALE +* +* Scale to avoid underflow +* + LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL + LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. + $ SMALL +* + SCALE = ONE + IF( LSA ) + $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) + IF( LSB ) + $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* + $ MIN( BNORM, BIG ) ) + IF( LSA .OR. LSB ) THEN + SCALE = MIN( SCALE, ONE / + $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), + $ ABS1( BCOEFF ) ) ) ) + IF( LSA ) THEN + ACOEFF = ASCALE*( SCALE*SBETA ) + ELSE + ACOEFF = SCALE*ACOEFF + END IF + IF( LSB ) THEN + BCOEFF = BSCALE*( SCALE*SALPHA ) + ELSE + BCOEFF = SCALE*BCOEFF + END IF + END IF +* + ACOEFA = ABS( ACOEFF ) + BCOEFA = ABS1( BCOEFF ) + XMAX = ONE + DO 160 JR = 1, N + WORK( JR ) = CZERO + 160 CONTINUE + WORK( JE ) = CONE + DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) +* +* Triangular solve of (a A - b B) x = 0 (columnwise) +* +* WORK(1:j-1) contains sums w, +* WORK(j+1:JE) contains x +* + DO 170 JR = 1, JE - 1 + WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE ) + 170 CONTINUE + WORK( JE ) = CONE +* + DO 210 J = JE - 1, 1, -1 +* +* Form x(j) := - w(j) / d +* with scaling and perturbation of the denominator +* + D = ACOEFF*S( J, J ) - BCOEFF*P( J, J ) + IF( ABS1( D ).LE.DMIN ) + $ D = DCMPLX( DMIN ) +* + IF( ABS1( D ).LT.ONE ) THEN + IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN + TEMP = ONE / ABS1( WORK( J ) ) + DO 180 JR = 1, JE + WORK( JR ) = TEMP*WORK( JR ) + 180 CONTINUE + END IF + END IF +* + WORK( J ) = ZLADIV( -WORK( J ), D ) +* + IF( J.GT.1 ) THEN +* +* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling +* + IF( ABS1( WORK( J ) ).GT.ONE ) THEN + TEMP = ONE / ABS1( WORK( J ) ) + IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE. + $ BIGNUM*TEMP ) THEN + DO 190 JR = 1, JE + WORK( JR ) = TEMP*WORK( JR ) + 190 CONTINUE + END IF + END IF +* + CA = ACOEFF*WORK( J ) + CB = BCOEFF*WORK( J ) + DO 200 JR = 1, J - 1 + WORK( JR ) = WORK( JR ) + CA*S( JR, J ) - + $ CB*P( JR, J ) + 200 CONTINUE + END IF + 210 CONTINUE +* +* Back transform eigenvector if HOWMNY='B'. +* + IF( ILBACK ) THEN + CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1, + $ CZERO, WORK( N+1 ), 1 ) + ISRC = 2 + IEND = N + ELSE + ISRC = 1 + IEND = JE + END IF +* +* Copy and scale eigenvector into column of VR +* + XMAX = ZERO + DO 220 JR = 1, IEND + XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) + 220 CONTINUE +* + IF( XMAX.GT.SAFMIN ) THEN + TEMP = ONE / XMAX + DO 230 JR = 1, IEND + VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) + 230 CONTINUE + ELSE + IEND = 0 + END IF +* + DO 240 JR = IEND + 1, N + VR( JR, IEIG ) = CZERO + 240 CONTINUE +* + END IF + 250 CONTINUE + END IF +* + RETURN +* +* End of ZTGEVC +* + END