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