![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgghd3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dgghd3.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: *> \brief \b DGGHD3
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGGHD3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgghd3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgghd3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgghd3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
22: * LDQ, Z, LDZ, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER COMPQ, COMPZ
26: * INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30: * $ Z( LDZ, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DGGHD3 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 DGGHD3 reduces the original
66: *> problem to generalized Hessenberg form.
67: *>
68: *> This is a blocked variant of DGGHRD, using matrix-matrix
69: *> multiplications for parts of the computation to enhance performance.
70: *> \endverbatim
71: *
72: * Arguments:
73: * ==========
74: *
75: *> \param[in] COMPQ
76: *> \verbatim
77: *> COMPQ is CHARACTER*1
78: *> = 'N': do not compute Q;
79: *> = 'I': Q is initialized to the unit matrix, and the
80: *> orthogonal matrix Q is returned;
81: *> = 'V': Q must contain an orthogonal matrix Q1 on entry,
82: *> and the product Q1*Q is returned.
83: *> \endverbatim
84: *>
85: *> \param[in] COMPZ
86: *> \verbatim
87: *> COMPZ is CHARACTER*1
88: *> = 'N': do not compute Z;
89: *> = 'I': Z is initialized to the unit matrix, and the
90: *> orthogonal matrix Z is returned;
91: *> = 'V': Z must contain an orthogonal matrix Z1 on entry,
92: *> and the product Z1*Z is returned.
93: *> \endverbatim
94: *>
95: *> \param[in] N
96: *> \verbatim
97: *> N is INTEGER
98: *> The order of the matrices A and B. N >= 0.
99: *> \endverbatim
100: *>
101: *> \param[in] ILO
102: *> \verbatim
103: *> ILO is INTEGER
104: *> \endverbatim
105: *>
106: *> \param[in] IHI
107: *> \verbatim
108: *> IHI is INTEGER
109: *>
110: *> ILO and IHI mark the rows and columns of A which are to be
111: *> reduced. It is assumed that A is already upper triangular
112: *> in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
113: *> normally set by a previous call to DGGBAL; otherwise they
114: *> should be set to 1 and N respectively.
115: *> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
116: *> \endverbatim
117: *>
118: *> \param[in,out] A
119: *> \verbatim
120: *> A is DOUBLE PRECISION array, dimension (LDA, N)
121: *> On entry, the N-by-N general matrix to be reduced.
122: *> On exit, the upper triangle and the first subdiagonal of A
123: *> are overwritten with the upper Hessenberg matrix H, and the
124: *> rest is set to zero.
125: *> \endverbatim
126: *>
127: *> \param[in] LDA
128: *> \verbatim
129: *> LDA is INTEGER
130: *> The leading dimension of the array A. LDA >= max(1,N).
131: *> \endverbatim
132: *>
133: *> \param[in,out] B
134: *> \verbatim
135: *> B is DOUBLE PRECISION array, dimension (LDB, N)
136: *> On entry, the N-by-N upper triangular matrix B.
137: *> On exit, the upper triangular matrix T = Q**T B Z. The
138: *> elements below the diagonal are set to zero.
139: *> \endverbatim
140: *>
141: *> \param[in] LDB
142: *> \verbatim
143: *> LDB is INTEGER
144: *> The leading dimension of the array B. LDB >= max(1,N).
145: *> \endverbatim
146: *>
147: *> \param[in,out] Q
148: *> \verbatim
149: *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
150: *> On entry, if COMPQ = 'V', the orthogonal matrix Q1,
151: *> typically from the QR factorization of B.
152: *> On exit, if COMPQ='I', the orthogonal matrix Q, and if
153: *> COMPQ = 'V', the product Q1*Q.
154: *> Not referenced if COMPQ='N'.
155: *> \endverbatim
156: *>
157: *> \param[in] LDQ
158: *> \verbatim
159: *> LDQ is INTEGER
160: *> The leading dimension of the array Q.
161: *> LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
162: *> \endverbatim
163: *>
164: *> \param[in,out] Z
165: *> \verbatim
166: *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
167: *> On entry, if COMPZ = 'V', the orthogonal matrix Z1.
168: *> On exit, if COMPZ='I', the orthogonal matrix Z, and if
169: *> COMPZ = 'V', the product Z1*Z.
170: *> Not referenced if COMPZ='N'.
171: *> \endverbatim
172: *>
173: *> \param[in] LDZ
174: *> \verbatim
175: *> LDZ is INTEGER
176: *> The leading dimension of the array Z.
177: *> LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
178: *> \endverbatim
179: *>
180: *> \param[out] WORK
181: *> \verbatim
182: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
183: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
184: *> \endverbatim
185: *>
186: *> \param[in] LWORK
187: *> \verbatim
188: *> LWORK is INTEGER
189: *> The length of the array WORK. LWORK >= 1.
190: *> For optimum performance LWORK >= 6*N*NB, where NB is the
191: *> optimal blocksize.
192: *>
193: *> If LWORK = -1, then a workspace query is assumed; the routine
194: *> only calculates the optimal size of the WORK array, returns
195: *> this value as the first entry of the WORK array, and no error
196: *> message related to LWORK is issued by XERBLA.
197: *> \endverbatim
198: *>
199: *> \param[out] INFO
200: *> \verbatim
201: *> INFO is INTEGER
202: *> = 0: successful exit.
203: *> < 0: if INFO = -i, the i-th argument had an illegal value.
204: *> \endverbatim
205: *
206: * Authors:
207: * ========
208: *
209: *> \author Univ. of Tennessee
210: *> \author Univ. of California Berkeley
211: *> \author Univ. of Colorado Denver
212: *> \author NAG Ltd.
213: *
214: *> \date January 2015
215: *
216: *> \ingroup doubleOTHERcomputational
217: *
218: *> \par Further Details:
219: * =====================
220: *>
221: *> \verbatim
222: *>
223: *> This routine reduces A to Hessenberg form and maintains B in
224: *> using a blocked variant of Moler and Stewart's original algorithm,
225: *> as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
226: *> (BIT 2008).
227: *> \endverbatim
228: *>
229: * =====================================================================
230: SUBROUTINE DGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
231: $ LDQ, Z, LDZ, WORK, LWORK, INFO )
232: *
1.5 bertrand 233: * -- LAPACK computational routine (version 3.8.0) --
1.1 bertrand 234: * -- LAPACK is a software package provided by Univ. of Tennessee, --
235: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
236: * January 2015
237: *
238: IMPLICIT NONE
239: *
240: * .. Scalar Arguments ..
241: CHARACTER COMPQ, COMPZ
242: INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
243: * ..
244: * .. Array Arguments ..
245: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
246: $ Z( LDZ, * ), WORK( * )
247: * ..
248: *
249: * =====================================================================
250: *
251: * .. Parameters ..
252: DOUBLE PRECISION ZERO, ONE
253: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
254: * ..
255: * .. Local Scalars ..
256: LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
257: CHARACTER*1 COMPQ2, COMPZ2
258: INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
259: $ KACC22, LEN, LWKOPT, N2NB, NB, NBLST, NBMIN,
260: $ NH, NNB, NX, PPW, PPWO, PW, TOP, TOPQ
261: DOUBLE PRECISION C, C1, C2, S, S1, S2, TEMP, TEMP1, TEMP2, TEMP3
262: * ..
263: * .. External Functions ..
264: LOGICAL LSAME
265: INTEGER ILAENV
266: EXTERNAL ILAENV, LSAME
267: * ..
268: * .. External Subroutines ..
1.5 bertrand 269: EXTERNAL DGGHRD, DLARTG, DLASET, DORM22, DROT, DGEMM,
270: $ DGEMV, DTRMV, DLACPY, XERBLA
1.1 bertrand 271: * ..
272: * .. Intrinsic Functions ..
273: INTRINSIC DBLE, MAX
274: * ..
275: * .. Executable Statements ..
276: *
277: * Decode and test the input parameters.
278: *
279: INFO = 0
280: NB = ILAENV( 1, 'DGGHD3', ' ', N, ILO, IHI, -1 )
1.2 bertrand 281: LWKOPT = MAX( 6*N*NB, 1 )
1.1 bertrand 282: WORK( 1 ) = DBLE( LWKOPT )
283: INITQ = LSAME( COMPQ, 'I' )
284: WANTQ = INITQ .OR. LSAME( COMPQ, 'V' )
285: INITZ = LSAME( COMPZ, 'I' )
286: WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
287: LQUERY = ( LWORK.EQ.-1 )
288: *
289: IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
290: INFO = -1
291: ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
292: INFO = -2
293: ELSE IF( N.LT.0 ) THEN
294: INFO = -3
295: ELSE IF( ILO.LT.1 ) THEN
296: INFO = -4
297: ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
298: INFO = -5
299: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
300: INFO = -7
301: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
302: INFO = -9
303: ELSE IF( ( WANTQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
304: INFO = -11
305: ELSE IF( ( WANTZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
306: INFO = -13
307: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
308: INFO = -15
309: END IF
310: IF( INFO.NE.0 ) THEN
311: CALL XERBLA( 'DGGHD3', -INFO )
312: RETURN
313: ELSE IF( LQUERY ) THEN
314: RETURN
315: END IF
316: *
317: * Initialize Q and Z if desired.
318: *
319: IF( INITQ )
320: $ CALL DLASET( 'All', N, N, ZERO, ONE, Q, LDQ )
321: IF( INITZ )
322: $ CALL DLASET( 'All', N, N, ZERO, ONE, Z, LDZ )
323: *
324: * Zero out lower triangle of B.
325: *
326: IF( N.GT.1 )
327: $ CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2, 1), LDB )
328: *
329: * Quick return if possible
330: *
331: NH = IHI - ILO + 1
332: IF( NH.LE.1 ) THEN
333: WORK( 1 ) = ONE
334: RETURN
335: END IF
336: *
337: * Determine the blocksize.
338: *
339: NBMIN = ILAENV( 2, 'DGGHD3', ' ', N, ILO, IHI, -1 )
340: IF( NB.GT.1 .AND. NB.LT.NH ) THEN
341: *
342: * Determine when to use unblocked instead of blocked code.
343: *
344: NX = MAX( NB, ILAENV( 3, 'DGGHD3', ' ', N, ILO, IHI, -1 ) )
345: IF( NX.LT.NH ) THEN
346: *
347: * Determine if workspace is large enough for blocked code.
348: *
349: IF( LWORK.LT.LWKOPT ) THEN
350: *
351: * Not enough workspace to use optimal NB: determine the
352: * minimum value of NB, and reduce NB or force use of
353: * unblocked code.
354: *
355: NBMIN = MAX( 2, ILAENV( 2, 'DGGHD3', ' ', N, ILO, IHI,
356: $ -1 ) )
357: IF( LWORK.GE.6*N*NBMIN ) THEN
358: NB = LWORK / ( 6*N )
359: ELSE
360: NB = 1
361: END IF
362: END IF
363: END IF
364: END IF
365: *
366: IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
367: *
368: * Use unblocked code below
369: *
370: JCOL = ILO
371: *
372: ELSE
373: *
374: * Use blocked code
375: *
376: KACC22 = ILAENV( 16, 'DGGHD3', ' ', N, ILO, IHI, -1 )
377: BLK22 = KACC22.EQ.2
378: DO JCOL = ILO, IHI-2, NB
379: NNB = MIN( NB, IHI-JCOL-1 )
380: *
381: * Initialize small orthogonal factors that will hold the
382: * accumulated Givens rotations in workspace.
383: * N2NB denotes the number of 2*NNB-by-2*NNB factors
384: * NBLST denotes the (possibly smaller) order of the last
385: * factor.
386: *
387: N2NB = ( IHI-JCOL-1 ) / NNB - 1
388: NBLST = IHI - JCOL - N2NB*NNB
389: CALL DLASET( 'All', NBLST, NBLST, ZERO, ONE, WORK, NBLST )
390: PW = NBLST * NBLST + 1
391: DO I = 1, N2NB
392: CALL DLASET( 'All', 2*NNB, 2*NNB, ZERO, ONE,
393: $ WORK( PW ), 2*NNB )
394: PW = PW + 4*NNB*NNB
395: END DO
396: *
397: * Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
398: *
399: DO J = JCOL, JCOL+NNB-1
400: *
401: * Reduce Jth column of A. Store cosines and sines in Jth
402: * column of A and B, respectively.
403: *
404: DO I = IHI, J+2, -1
405: TEMP = A( I-1, J )
406: CALL DLARTG( TEMP, A( I, J ), C, S, A( I-1, J ) )
407: A( I, J ) = C
408: B( I, J ) = S
409: END DO
410: *
411: * Accumulate Givens rotations into workspace array.
412: *
413: PPW = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
414: LEN = 2 + J - JCOL
415: JROW = J + N2NB*NNB + 2
416: DO I = IHI, JROW, -1
417: C = A( I, J )
418: S = B( I, J )
419: DO JJ = PPW, PPW+LEN-1
420: TEMP = WORK( JJ + NBLST )
421: WORK( JJ + NBLST ) = C*TEMP - S*WORK( JJ )
422: WORK( JJ ) = S*TEMP + C*WORK( JJ )
423: END DO
424: LEN = LEN + 1
425: PPW = PPW - NBLST - 1
426: END DO
427: *
428: PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
429: J0 = JROW - NNB
430: DO JROW = J0, J+2, -NNB
431: PPW = PPWO
432: LEN = 2 + J - JCOL
433: DO I = JROW+NNB-1, JROW, -1
434: C = A( I, J )
435: S = B( I, J )
436: DO JJ = PPW, PPW+LEN-1
437: TEMP = WORK( JJ + 2*NNB )
438: WORK( JJ + 2*NNB ) = C*TEMP - S*WORK( JJ )
439: WORK( JJ ) = S*TEMP + C*WORK( JJ )
440: END DO
441: LEN = LEN + 1
442: PPW = PPW - 2*NNB - 1
443: END DO
444: PPWO = PPWO + 4*NNB*NNB
445: END DO
446: *
447: * TOP denotes the number of top rows in A and B that will
448: * not be updated during the next steps.
449: *
450: IF( JCOL.LE.2 ) THEN
451: TOP = 0
452: ELSE
453: TOP = JCOL
454: END IF
455: *
456: * Propagate transformations through B and replace stored
457: * left sines/cosines by right sines/cosines.
458: *
459: DO JJ = N, J+1, -1
460: *
461: * Update JJth column of B.
462: *
463: DO I = MIN( JJ+1, IHI ), J+2, -1
464: C = A( I, J )
465: S = B( I, J )
466: TEMP = B( I, JJ )
467: B( I, JJ ) = C*TEMP - S*B( I-1, JJ )
468: B( I-1, JJ ) = S*TEMP + C*B( I-1, JJ )
469: END DO
470: *
471: * Annihilate B( JJ+1, JJ ).
472: *
473: IF( JJ.LT.IHI ) THEN
474: TEMP = B( JJ+1, JJ+1 )
475: CALL DLARTG( TEMP, B( JJ+1, JJ ), C, S,
476: $ B( JJ+1, JJ+1 ) )
477: B( JJ+1, JJ ) = ZERO
478: CALL DROT( JJ-TOP, B( TOP+1, JJ+1 ), 1,
479: $ B( TOP+1, JJ ), 1, C, S )
480: A( JJ+1, J ) = C
481: B( JJ+1, J ) = -S
482: END IF
483: END DO
484: *
485: * Update A by transformations from right.
486: * Explicit loop unrolling provides better performance
487: * compared to DLASR.
488: * CALL DLASR( 'Right', 'Variable', 'Backward', IHI-TOP,
489: * $ IHI-J, A( J+2, J ), B( J+2, J ),
490: * $ A( TOP+1, J+1 ), LDA )
491: *
492: JJ = MOD( IHI-J-1, 3 )
493: DO I = IHI-J-3, JJ+1, -3
494: C = A( J+1+I, J )
495: S = -B( J+1+I, J )
496: C1 = A( J+2+I, J )
497: S1 = -B( J+2+I, J )
498: C2 = A( J+3+I, J )
499: S2 = -B( J+3+I, J )
500: *
501: DO K = TOP+1, IHI
502: TEMP = A( K, J+I )
503: TEMP1 = A( K, J+I+1 )
504: TEMP2 = A( K, J+I+2 )
505: TEMP3 = A( K, J+I+3 )
506: A( K, J+I+3 ) = C2*TEMP3 + S2*TEMP2
507: TEMP2 = -S2*TEMP3 + C2*TEMP2
508: A( K, J+I+2 ) = C1*TEMP2 + S1*TEMP1
509: TEMP1 = -S1*TEMP2 + C1*TEMP1
510: A( K, J+I+1 ) = C*TEMP1 + S*TEMP
511: A( K, J+I ) = -S*TEMP1 + C*TEMP
512: END DO
513: END DO
514: *
515: IF( JJ.GT.0 ) THEN
516: DO I = JJ, 1, -1
517: CALL DROT( IHI-TOP, A( TOP+1, J+I+1 ), 1,
518: $ A( TOP+1, J+I ), 1, A( J+1+I, J ),
519: $ -B( J+1+I, J ) )
520: END DO
521: END IF
522: *
523: * Update (J+1)th column of A by transformations from left.
524: *
525: IF ( J .LT. JCOL + NNB - 1 ) THEN
526: LEN = 1 + J - JCOL
527: *
528: * Multiply with the trailing accumulated orthogonal
529: * matrix, which takes the form
530: *
531: * [ U11 U12 ]
532: * U = [ ],
533: * [ U21 U22 ]
534: *
535: * where U21 is a LEN-by-LEN matrix and U12 is lower
536: * triangular.
537: *
538: JROW = IHI - NBLST + 1
539: CALL DGEMV( 'Transpose', NBLST, LEN, ONE, WORK,
540: $ NBLST, A( JROW, J+1 ), 1, ZERO,
541: $ WORK( PW ), 1 )
542: PPW = PW + LEN
543: DO I = JROW, JROW+NBLST-LEN-1
544: WORK( PPW ) = A( I, J+1 )
545: PPW = PPW + 1
546: END DO
547: CALL DTRMV( 'Lower', 'Transpose', 'Non-unit',
548: $ NBLST-LEN, WORK( LEN*NBLST + 1 ), NBLST,
549: $ WORK( PW+LEN ), 1 )
550: CALL DGEMV( 'Transpose', LEN, NBLST-LEN, ONE,
551: $ WORK( (LEN+1)*NBLST - LEN + 1 ), NBLST,
552: $ A( JROW+NBLST-LEN, J+1 ), 1, ONE,
553: $ WORK( PW+LEN ), 1 )
554: PPW = PW
555: DO I = JROW, JROW+NBLST-1
556: A( I, J+1 ) = WORK( PPW )
557: PPW = PPW + 1
558: END DO
559: *
560: * Multiply with the other accumulated orthogonal
561: * matrices, which take the form
562: *
563: * [ U11 U12 0 ]
564: * [ ]
565: * U = [ U21 U22 0 ],
566: * [ ]
567: * [ 0 0 I ]
568: *
569: * where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
570: * matrix, U21 is a LEN-by-LEN upper triangular matrix
571: * and U12 is an NNB-by-NNB lower triangular matrix.
572: *
573: PPWO = 1 + NBLST*NBLST
574: J0 = JROW - NNB
575: DO JROW = J0, JCOL+1, -NNB
576: PPW = PW + LEN
577: DO I = JROW, JROW+NNB-1
578: WORK( PPW ) = A( I, J+1 )
579: PPW = PPW + 1
580: END DO
581: PPW = PW
582: DO I = JROW+NNB, JROW+NNB+LEN-1
583: WORK( PPW ) = A( I, J+1 )
584: PPW = PPW + 1
585: END DO
586: CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', LEN,
587: $ WORK( PPWO + NNB ), 2*NNB, WORK( PW ),
588: $ 1 )
589: CALL DTRMV( 'Lower', 'Transpose', 'Non-unit', NNB,
590: $ WORK( PPWO + 2*LEN*NNB ),
591: $ 2*NNB, WORK( PW + LEN ), 1 )
592: CALL DGEMV( 'Transpose', NNB, LEN, ONE,
593: $ WORK( PPWO ), 2*NNB, A( JROW, J+1 ), 1,
594: $ ONE, WORK( PW ), 1 )
595: CALL DGEMV( 'Transpose', LEN, NNB, ONE,
596: $ WORK( PPWO + 2*LEN*NNB + NNB ), 2*NNB,
597: $ A( JROW+NNB, J+1 ), 1, ONE,
598: $ WORK( PW+LEN ), 1 )
599: PPW = PW
600: DO I = JROW, JROW+LEN+NNB-1
601: A( I, J+1 ) = WORK( PPW )
602: PPW = PPW + 1
603: END DO
604: PPWO = PPWO + 4*NNB*NNB
605: END DO
606: END IF
607: END DO
608: *
609: * Apply accumulated orthogonal matrices to A.
610: *
611: COLA = N - JCOL - NNB + 1
612: J = IHI - NBLST + 1
613: CALL DGEMM( 'Transpose', 'No Transpose', NBLST,
614: $ COLA, NBLST, ONE, WORK, NBLST,
615: $ A( J, JCOL+NNB ), LDA, ZERO, WORK( PW ),
616: $ NBLST )
617: CALL DLACPY( 'All', NBLST, COLA, WORK( PW ), NBLST,
618: $ A( J, JCOL+NNB ), LDA )
619: PPWO = NBLST*NBLST + 1
620: J0 = J - NNB
621: DO J = J0, JCOL+1, -NNB
622: IF ( BLK22 ) THEN
623: *
624: * Exploit the structure of
625: *
626: * [ U11 U12 ]
627: * U = [ ]
628: * [ U21 U22 ],
629: *
630: * where all blocks are NNB-by-NNB, U21 is upper
631: * triangular and U12 is lower triangular.
632: *
633: CALL DORM22( 'Left', 'Transpose', 2*NNB, COLA, NNB,
634: $ NNB, WORK( PPWO ), 2*NNB,
635: $ A( J, JCOL+NNB ), LDA, WORK( PW ),
636: $ LWORK-PW+1, IERR )
637: ELSE
638: *
639: * Ignore the structure of U.
640: *
641: CALL DGEMM( 'Transpose', 'No Transpose', 2*NNB,
642: $ COLA, 2*NNB, ONE, WORK( PPWO ), 2*NNB,
643: $ A( J, JCOL+NNB ), LDA, ZERO, WORK( PW ),
644: $ 2*NNB )
645: CALL DLACPY( 'All', 2*NNB, COLA, WORK( PW ), 2*NNB,
646: $ A( J, JCOL+NNB ), LDA )
647: END IF
648: PPWO = PPWO + 4*NNB*NNB
649: END DO
650: *
651: * Apply accumulated orthogonal matrices to Q.
652: *
653: IF( WANTQ ) THEN
654: J = IHI - NBLST + 1
655: IF ( INITQ ) THEN
656: TOPQ = MAX( 2, J - JCOL + 1 )
657: NH = IHI - TOPQ + 1
658: ELSE
659: TOPQ = 1
660: NH = N
661: END IF
662: CALL DGEMM( 'No Transpose', 'No Transpose', NH,
663: $ NBLST, NBLST, ONE, Q( TOPQ, J ), LDQ,
664: $ WORK, NBLST, ZERO, WORK( PW ), NH )
665: CALL DLACPY( 'All', NH, NBLST, WORK( PW ), NH,
666: $ Q( TOPQ, J ), LDQ )
667: PPWO = NBLST*NBLST + 1
668: J0 = J - NNB
669: DO J = J0, JCOL+1, -NNB
670: IF ( INITQ ) THEN
671: TOPQ = MAX( 2, J - JCOL + 1 )
672: NH = IHI - TOPQ + 1
673: END IF
674: IF ( BLK22 ) THEN
675: *
676: * Exploit the structure of U.
677: *
678: CALL DORM22( 'Right', 'No Transpose', NH, 2*NNB,
679: $ NNB, NNB, WORK( PPWO ), 2*NNB,
680: $ Q( TOPQ, J ), LDQ, WORK( PW ),
681: $ LWORK-PW+1, IERR )
682: ELSE
683: *
684: * Ignore the structure of U.
685: *
686: CALL DGEMM( 'No Transpose', 'No Transpose', NH,
687: $ 2*NNB, 2*NNB, ONE, Q( TOPQ, J ), LDQ,
688: $ WORK( PPWO ), 2*NNB, ZERO, WORK( PW ),
689: $ NH )
690: CALL DLACPY( 'All', NH, 2*NNB, WORK( PW ), NH,
691: $ Q( TOPQ, J ), LDQ )
692: END IF
693: PPWO = PPWO + 4*NNB*NNB
694: END DO
695: END IF
696: *
697: * Accumulate right Givens rotations if required.
698: *
699: IF ( WANTZ .OR. TOP.GT.0 ) THEN
700: *
701: * Initialize small orthogonal factors that will hold the
702: * accumulated Givens rotations in workspace.
703: *
704: CALL DLASET( 'All', NBLST, NBLST, ZERO, ONE, WORK,
705: $ NBLST )
706: PW = NBLST * NBLST + 1
707: DO I = 1, N2NB
708: CALL DLASET( 'All', 2*NNB, 2*NNB, ZERO, ONE,
709: $ WORK( PW ), 2*NNB )
710: PW = PW + 4*NNB*NNB
711: END DO
712: *
713: * Accumulate Givens rotations into workspace array.
714: *
715: DO J = JCOL, JCOL+NNB-1
716: PPW = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
717: LEN = 2 + J - JCOL
718: JROW = J + N2NB*NNB + 2
719: DO I = IHI, JROW, -1
720: C = A( I, J )
721: A( I, J ) = ZERO
722: S = B( I, J )
723: B( I, J ) = ZERO
724: DO JJ = PPW, PPW+LEN-1
725: TEMP = WORK( JJ + NBLST )
726: WORK( JJ + NBLST ) = C*TEMP - S*WORK( JJ )
727: WORK( JJ ) = S*TEMP + C*WORK( JJ )
728: END DO
729: LEN = LEN + 1
730: PPW = PPW - NBLST - 1
731: END DO
732: *
733: PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
734: J0 = JROW - NNB
735: DO JROW = J0, J+2, -NNB
736: PPW = PPWO
737: LEN = 2 + J - JCOL
738: DO I = JROW+NNB-1, JROW, -1
739: C = A( I, J )
740: A( I, J ) = ZERO
741: S = B( I, J )
742: B( I, J ) = ZERO
743: DO JJ = PPW, PPW+LEN-1
744: TEMP = WORK( JJ + 2*NNB )
745: WORK( JJ + 2*NNB ) = C*TEMP - S*WORK( JJ )
746: WORK( JJ ) = S*TEMP + C*WORK( JJ )
747: END DO
748: LEN = LEN + 1
749: PPW = PPW - 2*NNB - 1
750: END DO
751: PPWO = PPWO + 4*NNB*NNB
752: END DO
753: END DO
754: ELSE
755: *
756: CALL DLASET( 'Lower', IHI - JCOL - 1, NNB, ZERO, ZERO,
757: $ A( JCOL + 2, JCOL ), LDA )
758: CALL DLASET( 'Lower', IHI - JCOL - 1, NNB, ZERO, ZERO,
759: $ B( JCOL + 2, JCOL ), LDB )
760: END IF
761: *
762: * Apply accumulated orthogonal matrices to A and B.
763: *
764: IF ( TOP.GT.0 ) THEN
765: J = IHI - NBLST + 1
766: CALL DGEMM( 'No Transpose', 'No Transpose', TOP,
767: $ NBLST, NBLST, ONE, A( 1, J ), LDA,
768: $ WORK, NBLST, ZERO, WORK( PW ), TOP )
769: CALL DLACPY( 'All', TOP, NBLST, WORK( PW ), TOP,
770: $ A( 1, J ), LDA )
771: PPWO = NBLST*NBLST + 1
772: J0 = J - NNB
773: DO J = J0, JCOL+1, -NNB
774: IF ( BLK22 ) THEN
775: *
776: * Exploit the structure of U.
777: *
778: CALL DORM22( 'Right', 'No Transpose', TOP, 2*NNB,
779: $ NNB, NNB, WORK( PPWO ), 2*NNB,
780: $ A( 1, J ), LDA, WORK( PW ),
781: $ LWORK-PW+1, IERR )
782: ELSE
783: *
784: * Ignore the structure of U.
785: *
786: CALL DGEMM( 'No Transpose', 'No Transpose', TOP,
787: $ 2*NNB, 2*NNB, ONE, A( 1, J ), LDA,
788: $ WORK( PPWO ), 2*NNB, ZERO,
789: $ WORK( PW ), TOP )
790: CALL DLACPY( 'All', TOP, 2*NNB, WORK( PW ), TOP,
791: $ A( 1, J ), LDA )
792: END IF
793: PPWO = PPWO + 4*NNB*NNB
794: END DO
795: *
796: J = IHI - NBLST + 1
797: CALL DGEMM( 'No Transpose', 'No Transpose', TOP,
798: $ NBLST, NBLST, ONE, B( 1, J ), LDB,
799: $ WORK, NBLST, ZERO, WORK( PW ), TOP )
800: CALL DLACPY( 'All', TOP, NBLST, WORK( PW ), TOP,
801: $ B( 1, J ), LDB )
802: PPWO = NBLST*NBLST + 1
803: J0 = J - NNB
804: DO J = J0, JCOL+1, -NNB
805: IF ( BLK22 ) THEN
806: *
807: * Exploit the structure of U.
808: *
809: CALL DORM22( 'Right', 'No Transpose', TOP, 2*NNB,
810: $ NNB, NNB, WORK( PPWO ), 2*NNB,
811: $ B( 1, J ), LDB, WORK( PW ),
812: $ LWORK-PW+1, IERR )
813: ELSE
814: *
815: * Ignore the structure of U.
816: *
817: CALL DGEMM( 'No Transpose', 'No Transpose', TOP,
818: $ 2*NNB, 2*NNB, ONE, B( 1, J ), LDB,
819: $ WORK( PPWO ), 2*NNB, ZERO,
820: $ WORK( PW ), TOP )
821: CALL DLACPY( 'All', TOP, 2*NNB, WORK( PW ), TOP,
822: $ B( 1, J ), LDB )
823: END IF
824: PPWO = PPWO + 4*NNB*NNB
825: END DO
826: END IF
827: *
828: * Apply accumulated orthogonal matrices to Z.
829: *
830: IF( WANTZ ) THEN
831: J = IHI - NBLST + 1
832: IF ( INITQ ) THEN
833: TOPQ = MAX( 2, J - JCOL + 1 )
834: NH = IHI - TOPQ + 1
835: ELSE
836: TOPQ = 1
837: NH = N
838: END IF
839: CALL DGEMM( 'No Transpose', 'No Transpose', NH,
840: $ NBLST, NBLST, ONE, Z( TOPQ, J ), LDZ,
841: $ WORK, NBLST, ZERO, WORK( PW ), NH )
842: CALL DLACPY( 'All', NH, NBLST, WORK( PW ), NH,
843: $ Z( TOPQ, J ), LDZ )
844: PPWO = NBLST*NBLST + 1
845: J0 = J - NNB
846: DO J = J0, JCOL+1, -NNB
847: IF ( INITQ ) THEN
848: TOPQ = MAX( 2, J - JCOL + 1 )
849: NH = IHI - TOPQ + 1
850: END IF
851: IF ( BLK22 ) THEN
852: *
853: * Exploit the structure of U.
854: *
855: CALL DORM22( 'Right', 'No Transpose', NH, 2*NNB,
856: $ NNB, NNB, WORK( PPWO ), 2*NNB,
857: $ Z( TOPQ, J ), LDZ, WORK( PW ),
858: $ LWORK-PW+1, IERR )
859: ELSE
860: *
861: * Ignore the structure of U.
862: *
863: CALL DGEMM( 'No Transpose', 'No Transpose', NH,
864: $ 2*NNB, 2*NNB, ONE, Z( TOPQ, J ), LDZ,
865: $ WORK( PPWO ), 2*NNB, ZERO, WORK( PW ),
866: $ NH )
867: CALL DLACPY( 'All', NH, 2*NNB, WORK( PW ), NH,
868: $ Z( TOPQ, J ), LDZ )
869: END IF
870: PPWO = PPWO + 4*NNB*NNB
871: END DO
872: END IF
873: END DO
874: END IF
875: *
876: * Use unblocked code to reduce the rest of the matrix
877: * Avoid re-initialization of modified Q and Z.
878: *
879: COMPQ2 = COMPQ
880: COMPZ2 = COMPZ
881: IF ( JCOL.NE.ILO ) THEN
882: IF ( WANTQ )
883: $ COMPQ2 = 'V'
884: IF ( WANTZ )
885: $ COMPZ2 = 'V'
886: END IF
887: *
888: IF ( JCOL.LT.IHI )
889: $ CALL DGGHRD( COMPQ2, COMPZ2, N, JCOL, IHI, A, LDA, B, LDB, Q,
890: $ LDQ, Z, LDZ, IERR )
891: WORK( 1 ) = DBLE( LWKOPT )
892: *
893: RETURN
894: *
895: * End of DGGHD3
896: *
897: END