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: *> \ingroup doubleOTHERcomputational
193: *
194: *> \par Further Details:
195: * =====================
196: *>
197: *> \verbatim
198: *>
199: *> This routine reduces A to Hessenberg and B to triangular form by
200: *> an unblocked reduction, as described in _Matrix_Computations_,
201: *> by Golub and Van Loan (Johns Hopkins Press.)
202: *> \endverbatim
203: *>
204: * =====================================================================
205: SUBROUTINE DGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
206: $ LDQ, Z, LDZ, INFO )
207: *
208: * -- LAPACK computational routine --
209: * -- LAPACK is a software package provided by Univ. of Tennessee, --
210: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211: *
212: * .. Scalar Arguments ..
213: CHARACTER COMPQ, COMPZ
214: INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
215: * ..
216: * .. Array Arguments ..
217: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
218: $ Z( LDZ, * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: DOUBLE PRECISION ONE, ZERO
225: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
226: * ..
227: * .. Local Scalars ..
228: LOGICAL ILQ, ILZ
229: INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
230: DOUBLE PRECISION C, S, TEMP
231: * ..
232: * .. External Functions ..
233: LOGICAL LSAME
234: EXTERNAL LSAME
235: * ..
236: * .. External Subroutines ..
237: EXTERNAL DLARTG, DLASET, DROT, XERBLA
238: * ..
239: * .. Intrinsic Functions ..
240: INTRINSIC MAX
241: * ..
242: * .. Executable Statements ..
243: *
244: * Decode COMPQ
245: *
246: IF( LSAME( COMPQ, 'N' ) ) THEN
247: ILQ = .FALSE.
248: ICOMPQ = 1
249: ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
250: ILQ = .TRUE.
251: ICOMPQ = 2
252: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
253: ILQ = .TRUE.
254: ICOMPQ = 3
255: ELSE
256: ICOMPQ = 0
257: END IF
258: *
259: * Decode COMPZ
260: *
261: IF( LSAME( COMPZ, 'N' ) ) THEN
262: ILZ = .FALSE.
263: ICOMPZ = 1
264: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
265: ILZ = .TRUE.
266: ICOMPZ = 2
267: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
268: ILZ = .TRUE.
269: ICOMPZ = 3
270: ELSE
271: ICOMPZ = 0
272: END IF
273: *
274: * Test the input parameters.
275: *
276: INFO = 0
277: IF( ICOMPQ.LE.0 ) THEN
278: INFO = -1
279: ELSE IF( ICOMPZ.LE.0 ) THEN
280: INFO = -2
281: ELSE IF( N.LT.0 ) THEN
282: INFO = -3
283: ELSE IF( ILO.LT.1 ) THEN
284: INFO = -4
285: ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
286: INFO = -5
287: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
288: INFO = -7
289: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
290: INFO = -9
291: ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
292: INFO = -11
293: ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
294: INFO = -13
295: END IF
296: IF( INFO.NE.0 ) THEN
297: CALL XERBLA( 'DGGHRD', -INFO )
298: RETURN
299: END IF
300: *
301: * Initialize Q and Z if desired.
302: *
303: IF( ICOMPQ.EQ.3 )
304: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
305: IF( ICOMPZ.EQ.3 )
306: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
307: *
308: * Quick return if possible
309: *
310: IF( N.LE.1 )
311: $ RETURN
312: *
313: * Zero out lower triangle of B
314: *
315: DO 20 JCOL = 1, N - 1
316: DO 10 JROW = JCOL + 1, N
317: B( JROW, JCOL ) = ZERO
318: 10 CONTINUE
319: 20 CONTINUE
320: *
321: * Reduce A and B
322: *
323: DO 40 JCOL = ILO, IHI - 2
324: *
325: DO 30 JROW = IHI, JCOL + 2, -1
326: *
327: * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
328: *
329: TEMP = A( JROW-1, JCOL )
330: CALL DLARTG( TEMP, A( JROW, JCOL ), C, S,
331: $ A( JROW-1, JCOL ) )
332: A( JROW, JCOL ) = ZERO
333: CALL DROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
334: $ A( JROW, JCOL+1 ), LDA, C, S )
335: CALL DROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
336: $ B( JROW, JROW-1 ), LDB, C, S )
337: IF( ILQ )
338: $ CALL DROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, S )
339: *
340: * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
341: *
342: TEMP = B( JROW, JROW )
343: CALL DLARTG( TEMP, B( JROW, JROW-1 ), C, S,
344: $ B( JROW, JROW ) )
345: B( JROW, JROW-1 ) = ZERO
346: CALL DROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
347: CALL DROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
348: $ S )
349: IF( ILZ )
350: $ CALL DROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
351: 30 CONTINUE
352: 40 CONTINUE
353: *
354: RETURN
355: *
356: * End of DGGHRD
357: *
358: 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