Annotation of rpl/lapack/lapack/zgghrd.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b ZGGHRD
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZGGHRD + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgghrd.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgghrd.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgghrd.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
! 22: * LDQ, Z, LDZ, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER COMPQ, COMPZ
! 26: * INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
! 30: * $ Z( LDZ, * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper
! 40: *> Hessenberg form using unitary transformations, where A is a
! 41: *> general matrix and B is upper triangular. The form of the
! 42: *> generalized eigenvalue problem is
! 43: *> A*x = lambda*B*x,
! 44: *> and B is typically made upper triangular by computing its QR
! 45: *> factorization and moving the unitary matrix Q to the left side
! 46: *> of the equation.
! 47: *>
! 48: *> This subroutine simultaneously reduces A to a Hessenberg matrix H:
! 49: *> Q**H*A*Z = H
! 50: *> and transforms B to another upper triangular matrix T:
! 51: *> Q**H*B*Z = T
! 52: *> in order to reduce the problem to its standard form
! 53: *> H*y = lambda*T*y
! 54: *> where y = Z**H*x.
! 55: *>
! 56: *> The unitary matrices Q and Z are determined as products of Givens
! 57: *> rotations. They may either be formed explicitly, or they may be
! 58: *> postmultiplied into input matrices Q1 and Z1, so that
! 59: *> Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
! 60: *> Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
! 61: *> If Q1 is the unitary matrix from the QR factorization of B in the
! 62: *> original equation A*x = lambda*B*x, then ZGGHRD reduces the original
! 63: *> problem to generalized Hessenberg form.
! 64: *> \endverbatim
! 65: *
! 66: * Arguments:
! 67: * ==========
! 68: *
! 69: *> \param[in] COMPQ
! 70: *> \verbatim
! 71: *> COMPQ is CHARACTER*1
! 72: *> = 'N': do not compute Q;
! 73: *> = 'I': Q is initialized to the unit matrix, and the
! 74: *> unitary matrix Q is returned;
! 75: *> = 'V': Q must contain a unitary matrix Q1 on entry,
! 76: *> and the product Q1*Q is returned.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] COMPZ
! 80: *> \verbatim
! 81: *> COMPZ is CHARACTER*1
! 82: *> = 'N': do not compute Q;
! 83: *> = 'I': Q is initialized to the unit matrix, and the
! 84: *> unitary matrix Q is returned;
! 85: *> = 'V': Q must contain a unitary matrix Q1 on entry,
! 86: *> and the product Q1*Q is returned.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] N
! 90: *> \verbatim
! 91: *> N is INTEGER
! 92: *> The order of the matrices A and B. N >= 0.
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[in] ILO
! 96: *> \verbatim
! 97: *> ILO is INTEGER
! 98: *> \endverbatim
! 99: *>
! 100: *> \param[in] IHI
! 101: *> \verbatim
! 102: *> IHI is INTEGER
! 103: *>
! 104: *> ILO and IHI mark the rows and columns of A which are to be
! 105: *> reduced. It is assumed that A is already upper triangular
! 106: *> in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
! 107: *> normally set by a previous call to ZGGBAL; otherwise they
! 108: *> should be set to 1 and N respectively.
! 109: *> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
! 110: *> \endverbatim
! 111: *>
! 112: *> \param[in,out] A
! 113: *> \verbatim
! 114: *> A is COMPLEX*16 array, dimension (LDA, N)
! 115: *> On entry, the N-by-N general matrix to be reduced.
! 116: *> On exit, the upper triangle and the first subdiagonal of A
! 117: *> are overwritten with the upper Hessenberg matrix H, and the
! 118: *> rest is set to zero.
! 119: *> \endverbatim
! 120: *>
! 121: *> \param[in] LDA
! 122: *> \verbatim
! 123: *> LDA is INTEGER
! 124: *> The leading dimension of the array A. LDA >= max(1,N).
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in,out] B
! 128: *> \verbatim
! 129: *> B is COMPLEX*16 array, dimension (LDB, N)
! 130: *> On entry, the N-by-N upper triangular matrix B.
! 131: *> On exit, the upper triangular matrix T = Q**H B Z. The
! 132: *> elements below the diagonal are set to zero.
! 133: *> \endverbatim
! 134: *>
! 135: *> \param[in] LDB
! 136: *> \verbatim
! 137: *> LDB is INTEGER
! 138: *> The leading dimension of the array B. LDB >= max(1,N).
! 139: *> \endverbatim
! 140: *>
! 141: *> \param[in,out] Q
! 142: *> \verbatim
! 143: *> Q is COMPLEX*16 array, dimension (LDQ, N)
! 144: *> On entry, if COMPQ = 'V', the unitary matrix Q1, typically
! 145: *> from the QR factorization of B.
! 146: *> On exit, if COMPQ='I', the unitary matrix Q, and if
! 147: *> COMPQ = 'V', the product Q1*Q.
! 148: *> Not referenced if COMPQ='N'.
! 149: *> \endverbatim
! 150: *>
! 151: *> \param[in] LDQ
! 152: *> \verbatim
! 153: *> LDQ is INTEGER
! 154: *> The leading dimension of the array Q.
! 155: *> LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
! 156: *> \endverbatim
! 157: *>
! 158: *> \param[in,out] Z
! 159: *> \verbatim
! 160: *> Z is COMPLEX*16 array, dimension (LDZ, N)
! 161: *> On entry, if COMPZ = 'V', the unitary matrix Z1.
! 162: *> On exit, if COMPZ='I', the unitary matrix Z, and if
! 163: *> COMPZ = 'V', the product Z1*Z.
! 164: *> Not referenced if COMPZ='N'.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[in] LDZ
! 168: *> \verbatim
! 169: *> LDZ is INTEGER
! 170: *> The leading dimension of the array Z.
! 171: *> LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[out] INFO
! 175: *> \verbatim
! 176: *> INFO is INTEGER
! 177: *> = 0: successful exit.
! 178: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 179: *> \endverbatim
! 180: *
! 181: * Authors:
! 182: * ========
! 183: *
! 184: *> \author Univ. of Tennessee
! 185: *> \author Univ. of California Berkeley
! 186: *> \author Univ. of Colorado Denver
! 187: *> \author NAG Ltd.
! 188: *
! 189: *> \date November 2011
! 190: *
! 191: *> \ingroup complex16OTHERcomputational
! 192: *
! 193: *> \par Further Details:
! 194: * =====================
! 195: *>
! 196: *> \verbatim
! 197: *>
! 198: *> This routine reduces A to Hessenberg and B to triangular form by
! 199: *> an unblocked reduction, as described in _Matrix_Computations_,
! 200: *> by Golub and van Loan (Johns Hopkins Press).
! 201: *> \endverbatim
! 202: *>
! 203: * =====================================================================
1.1 bertrand 204: SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
205: $ LDQ, Z, LDZ, INFO )
206: *
1.8 ! bertrand 207: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 208: * -- LAPACK is a software package provided by Univ. of Tennessee, --
209: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 210: * November 2011
1.1 bertrand 211: *
212: * .. Scalar Arguments ..
213: CHARACTER COMPQ, COMPZ
214: INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
215: * ..
216: * .. Array Arguments ..
217: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
218: $ Z( LDZ, * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: COMPLEX*16 CONE, CZERO
225: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
226: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
227: * ..
228: * .. Local Scalars ..
229: LOGICAL ILQ, ILZ
230: INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
231: DOUBLE PRECISION C
232: COMPLEX*16 CTEMP, S
233: * ..
234: * .. External Functions ..
235: LOGICAL LSAME
236: EXTERNAL LSAME
237: * ..
238: * .. External Subroutines ..
239: EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT
240: * ..
241: * .. Intrinsic Functions ..
242: INTRINSIC DCONJG, MAX
243: * ..
244: * .. Executable Statements ..
245: *
246: * Decode COMPQ
247: *
248: IF( LSAME( COMPQ, 'N' ) ) THEN
249: ILQ = .FALSE.
250: ICOMPQ = 1
251: ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
252: ILQ = .TRUE.
253: ICOMPQ = 2
254: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
255: ILQ = .TRUE.
256: ICOMPQ = 3
257: ELSE
258: ICOMPQ = 0
259: END IF
260: *
261: * Decode COMPZ
262: *
263: IF( LSAME( COMPZ, 'N' ) ) THEN
264: ILZ = .FALSE.
265: ICOMPZ = 1
266: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
267: ILZ = .TRUE.
268: ICOMPZ = 2
269: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
270: ILZ = .TRUE.
271: ICOMPZ = 3
272: ELSE
273: ICOMPZ = 0
274: END IF
275: *
276: * Test the input parameters.
277: *
278: INFO = 0
279: IF( ICOMPQ.LE.0 ) THEN
280: INFO = -1
281: ELSE IF( ICOMPZ.LE.0 ) THEN
282: INFO = -2
283: ELSE IF( N.LT.0 ) THEN
284: INFO = -3
285: ELSE IF( ILO.LT.1 ) THEN
286: INFO = -4
287: ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
288: INFO = -5
289: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
290: INFO = -7
291: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
292: INFO = -9
293: ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
294: INFO = -11
295: ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
296: INFO = -13
297: END IF
298: IF( INFO.NE.0 ) THEN
299: CALL XERBLA( 'ZGGHRD', -INFO )
300: RETURN
301: END IF
302: *
303: * Initialize Q and Z if desired.
304: *
305: IF( ICOMPQ.EQ.3 )
306: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
307: IF( ICOMPZ.EQ.3 )
308: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
309: *
310: * Quick return if possible
311: *
312: IF( N.LE.1 )
313: $ RETURN
314: *
315: * Zero out lower triangle of B
316: *
317: DO 20 JCOL = 1, N - 1
318: DO 10 JROW = JCOL + 1, N
319: B( JROW, JCOL ) = CZERO
320: 10 CONTINUE
321: 20 CONTINUE
322: *
323: * Reduce A and B
324: *
325: DO 40 JCOL = ILO, IHI - 2
326: *
327: DO 30 JROW = IHI, JCOL + 2, -1
328: *
329: * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
330: *
331: CTEMP = A( JROW-1, JCOL )
332: CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
333: $ A( JROW-1, JCOL ) )
334: A( JROW, JCOL ) = CZERO
335: CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
336: $ A( JROW, JCOL+1 ), LDA, C, S )
337: CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
338: $ B( JROW, JROW-1 ), LDB, C, S )
339: IF( ILQ )
340: $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
341: $ DCONJG( S ) )
342: *
343: * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
344: *
345: CTEMP = B( JROW, JROW )
346: CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
347: $ B( JROW, JROW ) )
348: B( JROW, JROW-1 ) = CZERO
349: CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
350: CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
351: $ S )
352: IF( ILZ )
353: $ CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
354: 30 CONTINUE
355: 40 CONTINUE
356: *
357: RETURN
358: *
359: * End of ZGGHRD
360: *
361: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgghrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack