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