Annotation of rpl/lapack/lapack/dgghrd.f, revision 1.10
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