NAME

SYNOPSIS

  SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, 
 *      ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
  CHARACTER * 1 JOB, COMPQ, COMPZ
  INTEGER N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO
  DOUBLE PRECISION A(LDA,*), B(LDB,*), ALPHAR(*), ALPHAI(*), BETA(*), Q(LDQ,*), Z(LDZ,*), WORK(*)
 
  SUBROUTINE DHGEQZ_64( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, 
 *      LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
  CHARACTER * 1 JOB, COMPQ, COMPZ
  INTEGER*8 N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO
  DOUBLE PRECISION A(LDA,*), B(LDB,*), ALPHAR(*), ALPHAI(*), BETA(*), Q(LDQ,*), Z(LDZ,*), WORK(*)
 

F95 INTERFACE

  SUBROUTINE HGEQZ( JOB, COMPQ, COMPZ, [N], ILO, IHI, A, [LDA], B, 
 *       [LDB], ALPHAR, ALPHAI, BETA, Q, [LDQ], Z, [LDZ], [WORK], [LWORK], 
 *       [INFO])
  CHARACTER(LEN=1) :: JOB, COMPQ, COMPZ
  INTEGER :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO
  REAL(8), DIMENSION(:) :: ALPHAR, ALPHAI, BETA, WORK
  REAL(8), DIMENSION(:,:) :: A, B, Q, Z
 
  SUBROUTINE HGEQZ_64( JOB, COMPQ, COMPZ, [N], ILO, IHI, A, [LDA], B, 
 *       [LDB], ALPHAR, ALPHAI, BETA, Q, [LDQ], Z, [LDZ], [WORK], [LWORK], 
 *       [INFO])
  CHARACTER(LEN=1) :: JOB, COMPQ, COMPZ
  INTEGER(8) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, INFO
  REAL(8), DIMENSION(:) :: ALPHAR, ALPHAI, BETA, WORK
  REAL(8), DIMENSION(:,:) :: A, B, Q, Z
 

C INTERFACE

void dhgeqz(char job, char compq, char compz, int n, int ilo, int ihi, double *a, int lda, double *b, int ldb, double *alphar, double *alphai, double *beta, double *q, int ldq, double *z, int ldz, int *info);

void dhgeqz_64(char job, char compq, char compz, long n, long ilo, long ihi, double *a, long lda, double *b, long ldb, double *alphar, double *alphai, double *beta, double *q, long ldq, double *z, long ldz, long *info);

PURPOSE

If JOB='S', then the pair (A,B) is simultaneously reduced to Schur form by applying one orthogonal tranformation (usually called Q) on the left and another (usually called Z) on the right. The 2-by-2 upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks of A will be reduced to positive diagonal matrices. (I.e., if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be positive.)

If JOB='E', then at each iteration, the same transformations are computed, but they are only applied to those parts of A and B which are needed to compute ALPHAR, ALPHAI, and BETAR.

If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal transformations used to reduce (A,B) are accumulated into the arrays Q and Z s.t.:

(in) A(in) Z(in)* = Q(out) A(out) Z(out)* (in) B(in) Z(in)* = Q(out) B(out) Z(out)*

Ref: C.B. Moler & G.W. Stewart, ``An Algorithm for Generalized Matrixigenvalue Problems'', SIAM J. Numer. Anal., 10(1973),p. 241--256.

ARGUMENTS

* JOB (input)

* COMPQ (input)

* COMPZ (input)

* N (input)

The order of the matrices A, B, Q, and Z. N >= 0.

* ILO (input)

It is assumed that A is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

* IHI (input)

See the description of ILO.

* A (input)

On entry, the N-by-N upper Hessenberg matrix A. Elements below the subdiagonal must be zero. If JOB='S', then on exit A and B will have been simultaneously reduced to generalized Schur form. If JOB='E', then on exit A will have been destroyed. The diagonal blocks will be correct, but the off-diagonal portion will be meaningless.

* LDA (input)

The leading dimension of the array A. LDA >= max( 1, N ).

* B (input)

On entry, the N-by-N upper triangular matrix B. Elements below the diagonal must be zero. 2-by-2 blocks in B corresponding to 2-by-2 blocks in A will be reduced to positive diagonal form. (I.e., if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be positive.) If JOB='S', then on exit A and B will have been simultaneously reduced to Schur form. If JOB='E', then on exit B will have been destroyed. Elements corresponding to diagonal blocks of A will be correct, but the off-diagonal portion will be meaningless.

* LDB (input)

The leading dimension of the array B. LDB >= max( 1, N ).

* ALPHAR (output)

ALPHAR(1:N) will be set to real parts of the diagonal elements of A that would result from reducing A and B to Schur form and then further reducing them both to triangular form using unitary transformations s.t. the diagonal of B was non-negative real. Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j). Note that the (real or complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized eigenvalues of the matrix pencil A - wB.

* ALPHAI (output)

ALPHAI(1:N) will be set to imaginary parts of the diagonal elements of A that would result from reducing A and B to Schur form and then further reducing them both to triangular form using unitary transformations s.t. the diagonal of B was non-negative real. Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0. Note that the (real or complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized eigenvalues of the matrix pencil A - wB.

* BETA (output)

BETA(1:N) will be set to the (real) diagonal elements of B that would result from reducing A and B to Schur form and then further reducing them both to triangular form using unitary transformations s.t. the diagonal of B was non-negative real. Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j). Note that the (real or complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized eigenvalues of the matrix pencil A - wB. (Note that BETA(1:N) will always be non-negative, and no BETAI is necessary.)

* Q (input/output)

If COMPQ='N', then Q will not be referenced. If COMPQ='V' or 'I', then the transpose of the orthogonal transformations which are applied to A and B on the left will be applied to the array Q on the right.

* LDQ (input)

The leading dimension of the array Q. LDQ >= 1. If COMPQ='V' or 'I', then LDQ >= N.

* Z (input/output)

If COMPZ='N', then Z will not be referenced. If COMPZ='V' or 'I', then the orthogonal transformations which are applied to A and B on the right will be applied to the array Z on the right.

* LDZ (input)

The leading dimension of the array Z. LDZ >= 1. If COMPZ='V' or 'I', then LDZ >= N.

* WORK (workspace)

On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.

* LWORK (input)

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.

* INFO (output)