001:       SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
002:      $                   LDQ, Z, LDZ, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          COMPQ, COMPZ
010:       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
014:      $                   Z( LDZ, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper
021: *  Hessenberg form using unitary transformations, where A is a
022: *  general matrix and B is upper triangular.  The form of the
023: *  generalized eigenvalue problem is
024: *     A*x = lambda*B*x,
025: *  and B is typically made upper triangular by computing its QR
026: *  factorization and moving the unitary matrix Q to the left side
027: *  of the equation.
028: *
029: *  This subroutine simultaneously reduces A to a Hessenberg matrix H:
030: *     Q**H*A*Z = H
031: *  and transforms B to another upper triangular matrix T:
032: *     Q**H*B*Z = T
033: *  in order to reduce the problem to its standard form
034: *     H*y = lambda*T*y
035: *  where y = Z**H*x.
036: *
037: *  The unitary matrices Q and Z are determined as products of Givens
038: *  rotations.  They may either be formed explicitly, or they may be
039: *  postmultiplied into input matrices Q1 and Z1, so that
040: *       Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
041: *       Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
042: *  If Q1 is the unitary matrix from the QR factorization of B in the
043: *  original equation A*x = lambda*B*x, then ZGGHRD reduces the original
044: *  problem to generalized Hessenberg form.
045: *
046: *  Arguments
047: *  =========
048: *
049: *  COMPQ   (input) CHARACTER*1
050: *          = 'N': do not compute Q;
051: *          = 'I': Q is initialized to the unit matrix, and the
052: *                 unitary matrix Q is returned;
053: *          = 'V': Q must contain a unitary matrix Q1 on entry,
054: *                 and the product Q1*Q is returned.
055: *
056: *  COMPZ   (input) CHARACTER*1
057: *          = 'N': do not compute Q;
058: *          = 'I': Q is initialized to the unit matrix, and the
059: *                 unitary matrix Q is returned;
060: *          = 'V': Q must contain a unitary matrix Q1 on entry,
061: *                 and the product Q1*Q is returned.
062: *
063: *  N       (input) INTEGER
064: *          The order of the matrices A and B.  N >= 0.
065: *
066: *  ILO     (input) INTEGER
067: *  IHI     (input) INTEGER
068: *          ILO and IHI mark the rows and columns of A which are to be
069: *          reduced.  It is assumed that A is already upper triangular
070: *          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
071: *          normally set by a previous call to ZGGBAL; otherwise they
072: *          should be set to 1 and N respectively.
073: *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
074: *
075: *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
076: *          On entry, the N-by-N general matrix to be reduced.
077: *          On exit, the upper triangle and the first subdiagonal of A
078: *          are overwritten with the upper Hessenberg matrix H, and the
079: *          rest is set to zero.
080: *
081: *  LDA     (input) INTEGER
082: *          The leading dimension of the array A.  LDA >= max(1,N).
083: *
084: *  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
085: *          On entry, the N-by-N upper triangular matrix B.
086: *          On exit, the upper triangular matrix T = Q**H B Z.  The
087: *          elements below the diagonal are set to zero.
088: *
089: *  LDB     (input) INTEGER
090: *          The leading dimension of the array B.  LDB >= max(1,N).
091: *
092: *  Q       (input/output) COMPLEX*16 array, dimension (LDQ, N)
093: *          On entry, if COMPQ = 'V', the unitary matrix Q1, typically
094: *          from the QR factorization of B.
095: *          On exit, if COMPQ='I', the unitary matrix Q, and if
096: *          COMPQ = 'V', the product Q1*Q.
097: *          Not referenced if COMPQ='N'.
098: *
099: *  LDQ     (input) INTEGER
100: *          The leading dimension of the array Q.
101: *          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
102: *
103: *  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N)
104: *          On entry, if COMPZ = 'V', the unitary matrix Z1.
105: *          On exit, if COMPZ='I', the unitary matrix Z, and if
106: *          COMPZ = 'V', the product Z1*Z.
107: *          Not referenced if COMPZ='N'.
108: *
109: *  LDZ     (input) INTEGER
110: *          The leading dimension of the array Z.
111: *          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
112: *
113: *  INFO    (output) INTEGER
114: *          = 0:  successful exit.
115: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
116: *
117: *  Further Details
118: *  ===============
119: *
120: *  This routine reduces A to Hessenberg and B to triangular form by
121: *  an unblocked reduction, as described in _Matrix_Computations_,
122: *  by Golub and van Loan (Johns Hopkins Press).
123: *
124: *  =====================================================================
125: *
126: *     .. Parameters ..
127:       COMPLEX*16         CONE, CZERO
128:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
129:      $                   CZERO = ( 0.0D+0, 0.0D+0 ) )
130: *     ..
131: *     .. Local Scalars ..
132:       LOGICAL            ILQ, ILZ
133:       INTEGER            ICOMPQ, ICOMPZ, JCOL, JROW
134:       DOUBLE PRECISION   C
135:       COMPLEX*16         CTEMP, S
136: *     ..
137: *     .. External Functions ..
138:       LOGICAL            LSAME
139:       EXTERNAL           LSAME
140: *     ..
141: *     .. External Subroutines ..
142:       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT
143: *     ..
144: *     .. Intrinsic Functions ..
145:       INTRINSIC          DCONJG, MAX
146: *     ..
147: *     .. Executable Statements ..
148: *
149: *     Decode COMPQ
150: *
151:       IF( LSAME( COMPQ, 'N' ) ) THEN
152:          ILQ = .FALSE.
153:          ICOMPQ = 1
154:       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
155:          ILQ = .TRUE.
156:          ICOMPQ = 2
157:       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
158:          ILQ = .TRUE.
159:          ICOMPQ = 3
160:       ELSE
161:          ICOMPQ = 0
162:       END IF
163: *
164: *     Decode COMPZ
165: *
166:       IF( LSAME( COMPZ, 'N' ) ) THEN
167:          ILZ = .FALSE.
168:          ICOMPZ = 1
169:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
170:          ILZ = .TRUE.
171:          ICOMPZ = 2
172:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
173:          ILZ = .TRUE.
174:          ICOMPZ = 3
175:       ELSE
176:          ICOMPZ = 0
177:       END IF
178: *
179: *     Test the input parameters.
180: *
181:       INFO = 0
182:       IF( ICOMPQ.LE.0 ) THEN
183:          INFO = -1
184:       ELSE IF( ICOMPZ.LE.0 ) THEN
185:          INFO = -2
186:       ELSE IF( N.LT.0 ) THEN
187:          INFO = -3
188:       ELSE IF( ILO.LT.1 ) THEN
189:          INFO = -4
190:       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
191:          INFO = -5
192:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
193:          INFO = -7
194:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
195:          INFO = -9
196:       ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
197:          INFO = -11
198:       ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
199:          INFO = -13
200:       END IF
201:       IF( INFO.NE.0 ) THEN
202:          CALL XERBLA( 'ZGGHRD', -INFO )
203:          RETURN
204:       END IF
205: *
206: *     Initialize Q and Z if desired.
207: *
208:       IF( ICOMPQ.EQ.3 )
209:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
210:       IF( ICOMPZ.EQ.3 )
211:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
212: *
213: *     Quick return if possible
214: *
215:       IF( N.LE.1 )
216:      $   RETURN
217: *
218: *     Zero out lower triangle of B
219: *
220:       DO 20 JCOL = 1, N - 1
221:          DO 10 JROW = JCOL + 1, N
222:             B( JROW, JCOL ) = CZERO
223:    10    CONTINUE
224:    20 CONTINUE
225: *
226: *     Reduce A and B
227: *
228:       DO 40 JCOL = ILO, IHI - 2
229: *
230:          DO 30 JROW = IHI, JCOL + 2, -1
231: *
232: *           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
233: *
234:             CTEMP = A( JROW-1, JCOL )
235:             CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
236:      $                   A( JROW-1, JCOL ) )
237:             A( JROW, JCOL ) = CZERO
238:             CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
239:      $                 A( JROW, JCOL+1 ), LDA, C, S )
240:             CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
241:      $                 B( JROW, JROW-1 ), LDB, C, S )
242:             IF( ILQ )
243:      $         CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
244:      $                    DCONJG( S ) )
245: *
246: *           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
247: *
248:             CTEMP = B( JROW, JROW )
249:             CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
250:      $                   B( JROW, JROW ) )
251:             B( JROW, JROW-1 ) = CZERO
252:             CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
253:             CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
254:      $                 S )
255:             IF( ILZ )
256:      $         CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
257:    30    CONTINUE
258:    40 CONTINUE
259: *
260:       RETURN
261: *
262: *     End of ZGGHRD
263: *
264:       END
265: