1: *> \brief \b ZGGHD3
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGGHD3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgghd3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgghd3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgghd3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGGHD3( 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: * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30: * $ Z( LDZ, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZGGHD3 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 ZGGHD3 reduces the original
63: *> problem to generalized Hessenberg form.
64: *>
65: *> This is a blocked variant of CGGHRD, using matrix-matrix
66: *> multiplications for parts of the computation to enhance performance.
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: *> unitary matrix Q is returned;
78: *> = 'V': Q must contain a unitary 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: *> unitary matrix Z is returned;
88: *> = 'V': Z must contain a unitary 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 ZGGBAL; 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDB, N)
133: *> On entry, the N-by-N upper triangular matrix B.
134: *> On exit, the upper triangular matrix T = Q**H 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 COMPLEX*16 array, dimension (LDQ, N)
147: *> On entry, if COMPQ = 'V', the unitary matrix Q1, typically
148: *> from the QR factorization of B.
149: *> On exit, if COMPQ='I', the unitary 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 COMPLEX*16 array, dimension (LDZ, N)
164: *> On entry, if COMPZ = 'V', the unitary matrix Z1.
165: *> On exit, if COMPZ='I', the unitary 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] WORK
178: *> \verbatim
179: *> WORK is COMPLEX*16 array, dimension (LWORK)
180: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
181: *> \endverbatim
182: *>
183: *> \param[in] LWORK
184: *> \verbatim
185: *> LWORK is INTEGER
186: *> The length of the array WORK. LWORK >= 1.
187: *> For optimum performance LWORK >= 6*N*NB, where NB is the
188: *> optimal blocksize.
189: *>
190: *> If LWORK = -1, then a workspace query is assumed; the routine
191: *> only calculates the optimal size of the WORK array, returns
192: *> this value as the first entry of the WORK array, and no error
193: *> message related to LWORK is issued by XERBLA.
194: *> \endverbatim
195: *>
196: *> \param[out] INFO
197: *> \verbatim
198: *> INFO is INTEGER
199: *> = 0: successful exit.
200: *> < 0: if INFO = -i, the i-th argument had an illegal value.
201: *> \endverbatim
202: *
203: * Authors:
204: * ========
205: *
206: *> \author Univ. of Tennessee
207: *> \author Univ. of California Berkeley
208: *> \author Univ. of Colorado Denver
209: *> \author NAG Ltd.
210: *
211: *> \date January 2015
212: *
213: *> \ingroup complex16OTHERcomputational
214: *
215: *> \par Further Details:
216: * =====================
217: *>
218: *> \verbatim
219: *>
220: *> This routine reduces A to Hessenberg form and maintains B in
221: *> using a blocked variant of Moler and Stewart's original algorithm,
222: *> as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
223: *> (BIT 2008).
224: *> \endverbatim
225: *>
226: * =====================================================================
227: SUBROUTINE ZGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
228: $ LDQ, Z, LDZ, WORK, LWORK, INFO )
229: *
230: * -- LAPACK computational routine (version 3.8.0) --
231: * -- LAPACK is a software package provided by Univ. of Tennessee, --
232: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233: * January 2015
234: *
235: IMPLICIT NONE
236: *
237: * .. Scalar Arguments ..
238: CHARACTER COMPQ, COMPZ
239: INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
240: * ..
241: * .. Array Arguments ..
242: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
243: $ Z( LDZ, * ), WORK( * )
244: * ..
245: *
246: * =====================================================================
247: *
248: * .. Parameters ..
249: COMPLEX*16 CONE, CZERO
250: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
251: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
252: * ..
253: * .. Local Scalars ..
254: LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
255: CHARACTER*1 COMPQ2, COMPZ2
256: INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
257: $ KACC22, LEN, LWKOPT, N2NB, NB, NBLST, NBMIN,
258: $ NH, NNB, NX, PPW, PPWO, PW, TOP, TOPQ
259: DOUBLE PRECISION C
260: COMPLEX*16 C1, C2, CTEMP, S, S1, S2, TEMP, TEMP1, TEMP2,
261: $ TEMP3
262: * ..
263: * .. External Functions ..
264: LOGICAL LSAME
265: INTEGER ILAENV
266: EXTERNAL ILAENV, LSAME
267: * ..
268: * .. External Subroutines ..
269: EXTERNAL ZGGHRD, ZLARTG, ZLASET, ZUNM22, ZROT, ZGEMM,
270: $ ZGEMV, ZTRMV, ZLACPY, XERBLA
271: * ..
272: * .. Intrinsic Functions ..
273: INTRINSIC DBLE, DCMPLX, DCONJG, MAX
274: * ..
275: * .. Executable Statements ..
276: *
277: * Decode and test the input parameters.
278: *
279: INFO = 0
280: NB = ILAENV( 1, 'ZGGHD3', ' ', N, ILO, IHI, -1 )
281: LWKOPT = MAX( 6*N*NB, 1 )
282: WORK( 1 ) = DCMPLX( 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( 'ZGGHD3', -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 ZLASET( 'All', N, N, CZERO, CONE, Q, LDQ )
321: IF( INITZ )
322: $ CALL ZLASET( 'All', N, N, CZERO, CONE, Z, LDZ )
323: *
324: * Zero out lower triangle of B.
325: *
326: IF( N.GT.1 )
327: $ CALL ZLASET( 'Lower', N-1, N-1, CZERO, CZERO, 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 ) = CONE
334: RETURN
335: END IF
336: *
337: * Determine the blocksize.
338: *
339: NBMIN = ILAENV( 2, 'ZGGHD3', ' ', 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, 'ZGGHD3', ' ', 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, 'ZGGHD3', ' ', 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, 'ZGGHD3', ' ', 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 unitary 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 ZLASET( 'All', NBLST, NBLST, CZERO, CONE, WORK, NBLST )
390: PW = NBLST * NBLST + 1
391: DO I = 1, N2NB
392: CALL ZLASET( 'All', 2*NNB, 2*NNB, CZERO, CONE,
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 ZLARTG( TEMP, A( I, J ), C, S, A( I-1, J ) )
407: A( I, J ) = DCMPLX( 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: CTEMP = A( I, J )
418: S = B( I, J )
419: DO JJ = PPW, PPW+LEN-1
420: TEMP = WORK( JJ + NBLST )
421: WORK( JJ + NBLST ) = CTEMP*TEMP - S*WORK( JJ )
422: WORK( JJ ) = DCONJG( S )*TEMP + CTEMP*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: CTEMP = 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 ) = CTEMP*TEMP - S*WORK( JJ )
439: WORK( JJ ) = DCONJG( S )*TEMP + CTEMP*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: CTEMP = A( I, J )
465: S = B( I, J )
466: TEMP = B( I, JJ )
467: B( I, JJ ) = CTEMP*TEMP - DCONJG( S )*B( I-1, JJ )
468: B( I-1, JJ ) = S*TEMP + CTEMP*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 ZLARTG( TEMP, B( JJ+1, JJ ), C, S,
476: $ B( JJ+1, JJ+1 ) )
477: B( JJ+1, JJ ) = CZERO
478: CALL ZROT( JJ-TOP, B( TOP+1, JJ+1 ), 1,
479: $ B( TOP+1, JJ ), 1, C, S )
480: A( JJ+1, J ) = DCMPLX( C )
481: B( JJ+1, J ) = -DCONJG( S )
482: END IF
483: END DO
484: *
485: * Update A by transformations from right.
486: *
487: JJ = MOD( IHI-J-1, 3 )
488: DO I = IHI-J-3, JJ+1, -3
489: CTEMP = A( J+1+I, J )
490: S = -B( J+1+I, J )
491: C1 = A( J+2+I, J )
492: S1 = -B( J+2+I, J )
493: C2 = A( J+3+I, J )
494: S2 = -B( J+3+I, J )
495: *
496: DO K = TOP+1, IHI
497: TEMP = A( K, J+I )
498: TEMP1 = A( K, J+I+1 )
499: TEMP2 = A( K, J+I+2 )
500: TEMP3 = A( K, J+I+3 )
501: A( K, J+I+3 ) = C2*TEMP3 + DCONJG( S2 )*TEMP2
502: TEMP2 = -S2*TEMP3 + C2*TEMP2
503: A( K, J+I+2 ) = C1*TEMP2 + DCONJG( S1 )*TEMP1
504: TEMP1 = -S1*TEMP2 + C1*TEMP1
505: A( K, J+I+1 ) = CTEMP*TEMP1 + DCONJG( S )*TEMP
506: A( K, J+I ) = -S*TEMP1 + CTEMP*TEMP
507: END DO
508: END DO
509: *
510: IF( JJ.GT.0 ) THEN
511: DO I = JJ, 1, -1
512: C = DBLE( A( J+1+I, J ) )
513: CALL ZROT( IHI-TOP, A( TOP+1, J+I+1 ), 1,
514: $ A( TOP+1, J+I ), 1, C,
515: $ -DCONJG( B( J+1+I, J ) ) )
516: END DO
517: END IF
518: *
519: * Update (J+1)th column of A by transformations from left.
520: *
521: IF ( J .LT. JCOL + NNB - 1 ) THEN
522: LEN = 1 + J - JCOL
523: *
524: * Multiply with the trailing accumulated unitary
525: * matrix, which takes the form
526: *
527: * [ U11 U12 ]
528: * U = [ ],
529: * [ U21 U22 ]
530: *
531: * where U21 is a LEN-by-LEN matrix and U12 is lower
532: * triangular.
533: *
534: JROW = IHI - NBLST + 1
535: CALL ZGEMV( 'Conjugate', NBLST, LEN, CONE, WORK,
536: $ NBLST, A( JROW, J+1 ), 1, CZERO,
537: $ WORK( PW ), 1 )
538: PPW = PW + LEN
539: DO I = JROW, JROW+NBLST-LEN-1
540: WORK( PPW ) = A( I, J+1 )
541: PPW = PPW + 1
542: END DO
543: CALL ZTRMV( 'Lower', 'Conjugate', 'Non-unit',
544: $ NBLST-LEN, WORK( LEN*NBLST + 1 ), NBLST,
545: $ WORK( PW+LEN ), 1 )
546: CALL ZGEMV( 'Conjugate', LEN, NBLST-LEN, CONE,
547: $ WORK( (LEN+1)*NBLST - LEN + 1 ), NBLST,
548: $ A( JROW+NBLST-LEN, J+1 ), 1, CONE,
549: $ WORK( PW+LEN ), 1 )
550: PPW = PW
551: DO I = JROW, JROW+NBLST-1
552: A( I, J+1 ) = WORK( PPW )
553: PPW = PPW + 1
554: END DO
555: *
556: * Multiply with the other accumulated unitary
557: * matrices, which take the form
558: *
559: * [ U11 U12 0 ]
560: * [ ]
561: * U = [ U21 U22 0 ],
562: * [ ]
563: * [ 0 0 I ]
564: *
565: * where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
566: * matrix, U21 is a LEN-by-LEN upper triangular matrix
567: * and U12 is an NNB-by-NNB lower triangular matrix.
568: *
569: PPWO = 1 + NBLST*NBLST
570: J0 = JROW - NNB
571: DO JROW = J0, JCOL+1, -NNB
572: PPW = PW + LEN
573: DO I = JROW, JROW+NNB-1
574: WORK( PPW ) = A( I, J+1 )
575: PPW = PPW + 1
576: END DO
577: PPW = PW
578: DO I = JROW+NNB, JROW+NNB+LEN-1
579: WORK( PPW ) = A( I, J+1 )
580: PPW = PPW + 1
581: END DO
582: CALL ZTRMV( 'Upper', 'Conjugate', 'Non-unit', LEN,
583: $ WORK( PPWO + NNB ), 2*NNB, WORK( PW ),
584: $ 1 )
585: CALL ZTRMV( 'Lower', 'Conjugate', 'Non-unit', NNB,
586: $ WORK( PPWO + 2*LEN*NNB ),
587: $ 2*NNB, WORK( PW + LEN ), 1 )
588: CALL ZGEMV( 'Conjugate', NNB, LEN, CONE,
589: $ WORK( PPWO ), 2*NNB, A( JROW, J+1 ), 1,
590: $ CONE, WORK( PW ), 1 )
591: CALL ZGEMV( 'Conjugate', LEN, NNB, CONE,
592: $ WORK( PPWO + 2*LEN*NNB + NNB ), 2*NNB,
593: $ A( JROW+NNB, J+1 ), 1, CONE,
594: $ WORK( PW+LEN ), 1 )
595: PPW = PW
596: DO I = JROW, JROW+LEN+NNB-1
597: A( I, J+1 ) = WORK( PPW )
598: PPW = PPW + 1
599: END DO
600: PPWO = PPWO + 4*NNB*NNB
601: END DO
602: END IF
603: END DO
604: *
605: * Apply accumulated unitary matrices to A.
606: *
607: COLA = N - JCOL - NNB + 1
608: J = IHI - NBLST + 1
609: CALL ZGEMM( 'Conjugate', 'No Transpose', NBLST,
610: $ COLA, NBLST, CONE, WORK, NBLST,
611: $ A( J, JCOL+NNB ), LDA, CZERO, WORK( PW ),
612: $ NBLST )
613: CALL ZLACPY( 'All', NBLST, COLA, WORK( PW ), NBLST,
614: $ A( J, JCOL+NNB ), LDA )
615: PPWO = NBLST*NBLST + 1
616: J0 = J - NNB
617: DO J = J0, JCOL+1, -NNB
618: IF ( BLK22 ) THEN
619: *
620: * Exploit the structure of
621: *
622: * [ U11 U12 ]
623: * U = [ ]
624: * [ U21 U22 ],
625: *
626: * where all blocks are NNB-by-NNB, U21 is upper
627: * triangular and U12 is lower triangular.
628: *
629: CALL ZUNM22( 'Left', 'Conjugate', 2*NNB, COLA, NNB,
630: $ NNB, WORK( PPWO ), 2*NNB,
631: $ A( J, JCOL+NNB ), LDA, WORK( PW ),
632: $ LWORK-PW+1, IERR )
633: ELSE
634: *
635: * Ignore the structure of U.
636: *
637: CALL ZGEMM( 'Conjugate', 'No Transpose', 2*NNB,
638: $ COLA, 2*NNB, CONE, WORK( PPWO ), 2*NNB,
639: $ A( J, JCOL+NNB ), LDA, CZERO, WORK( PW ),
640: $ 2*NNB )
641: CALL ZLACPY( 'All', 2*NNB, COLA, WORK( PW ), 2*NNB,
642: $ A( J, JCOL+NNB ), LDA )
643: END IF
644: PPWO = PPWO + 4*NNB*NNB
645: END DO
646: *
647: * Apply accumulated unitary matrices to Q.
648: *
649: IF( WANTQ ) THEN
650: J = IHI - NBLST + 1
651: IF ( INITQ ) THEN
652: TOPQ = MAX( 2, J - JCOL + 1 )
653: NH = IHI - TOPQ + 1
654: ELSE
655: TOPQ = 1
656: NH = N
657: END IF
658: CALL ZGEMM( 'No Transpose', 'No Transpose', NH,
659: $ NBLST, NBLST, CONE, Q( TOPQ, J ), LDQ,
660: $ WORK, NBLST, CZERO, WORK( PW ), NH )
661: CALL ZLACPY( 'All', NH, NBLST, WORK( PW ), NH,
662: $ Q( TOPQ, J ), LDQ )
663: PPWO = NBLST*NBLST + 1
664: J0 = J - NNB
665: DO J = J0, JCOL+1, -NNB
666: IF ( INITQ ) THEN
667: TOPQ = MAX( 2, J - JCOL + 1 )
668: NH = IHI - TOPQ + 1
669: END IF
670: IF ( BLK22 ) THEN
671: *
672: * Exploit the structure of U.
673: *
674: CALL ZUNM22( 'Right', 'No Transpose', NH, 2*NNB,
675: $ NNB, NNB, WORK( PPWO ), 2*NNB,
676: $ Q( TOPQ, J ), LDQ, WORK( PW ),
677: $ LWORK-PW+1, IERR )
678: ELSE
679: *
680: * Ignore the structure of U.
681: *
682: CALL ZGEMM( 'No Transpose', 'No Transpose', NH,
683: $ 2*NNB, 2*NNB, CONE, Q( TOPQ, J ), LDQ,
684: $ WORK( PPWO ), 2*NNB, CZERO, WORK( PW ),
685: $ NH )
686: CALL ZLACPY( 'All', NH, 2*NNB, WORK( PW ), NH,
687: $ Q( TOPQ, J ), LDQ )
688: END IF
689: PPWO = PPWO + 4*NNB*NNB
690: END DO
691: END IF
692: *
693: * Accumulate right Givens rotations if required.
694: *
695: IF ( WANTZ .OR. TOP.GT.0 ) THEN
696: *
697: * Initialize small unitary factors that will hold the
698: * accumulated Givens rotations in workspace.
699: *
700: CALL ZLASET( 'All', NBLST, NBLST, CZERO, CONE, WORK,
701: $ NBLST )
702: PW = NBLST * NBLST + 1
703: DO I = 1, N2NB
704: CALL ZLASET( 'All', 2*NNB, 2*NNB, CZERO, CONE,
705: $ WORK( PW ), 2*NNB )
706: PW = PW + 4*NNB*NNB
707: END DO
708: *
709: * Accumulate Givens rotations into workspace array.
710: *
711: DO J = JCOL, JCOL+NNB-1
712: PPW = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
713: LEN = 2 + J - JCOL
714: JROW = J + N2NB*NNB + 2
715: DO I = IHI, JROW, -1
716: CTEMP = A( I, J )
717: A( I, J ) = CZERO
718: S = B( I, J )
719: B( I, J ) = CZERO
720: DO JJ = PPW, PPW+LEN-1
721: TEMP = WORK( JJ + NBLST )
722: WORK( JJ + NBLST ) = CTEMP*TEMP -
723: $ DCONJG( S )*WORK( JJ )
724: WORK( JJ ) = S*TEMP + CTEMP*WORK( JJ )
725: END DO
726: LEN = LEN + 1
727: PPW = PPW - NBLST - 1
728: END DO
729: *
730: PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
731: J0 = JROW - NNB
732: DO JROW = J0, J+2, -NNB
733: PPW = PPWO
734: LEN = 2 + J - JCOL
735: DO I = JROW+NNB-1, JROW, -1
736: CTEMP = A( I, J )
737: A( I, J ) = CZERO
738: S = B( I, J )
739: B( I, J ) = CZERO
740: DO JJ = PPW, PPW+LEN-1
741: TEMP = WORK( JJ + 2*NNB )
742: WORK( JJ + 2*NNB ) = CTEMP*TEMP -
743: $ DCONJG( S )*WORK( JJ )
744: WORK( JJ ) = S*TEMP + CTEMP*WORK( JJ )
745: END DO
746: LEN = LEN + 1
747: PPW = PPW - 2*NNB - 1
748: END DO
749: PPWO = PPWO + 4*NNB*NNB
750: END DO
751: END DO
752: ELSE
753: *
754: CALL ZLASET( 'Lower', IHI - JCOL - 1, NNB, CZERO, CZERO,
755: $ A( JCOL + 2, JCOL ), LDA )
756: CALL ZLASET( 'Lower', IHI - JCOL - 1, NNB, CZERO, CZERO,
757: $ B( JCOL + 2, JCOL ), LDB )
758: END IF
759: *
760: * Apply accumulated unitary matrices to A and B.
761: *
762: IF ( TOP.GT.0 ) THEN
763: J = IHI - NBLST + 1
764: CALL ZGEMM( 'No Transpose', 'No Transpose', TOP,
765: $ NBLST, NBLST, CONE, A( 1, J ), LDA,
766: $ WORK, NBLST, CZERO, WORK( PW ), TOP )
767: CALL ZLACPY( 'All', TOP, NBLST, WORK( PW ), TOP,
768: $ A( 1, J ), LDA )
769: PPWO = NBLST*NBLST + 1
770: J0 = J - NNB
771: DO J = J0, JCOL+1, -NNB
772: IF ( BLK22 ) THEN
773: *
774: * Exploit the structure of U.
775: *
776: CALL ZUNM22( 'Right', 'No Transpose', TOP, 2*NNB,
777: $ NNB, NNB, WORK( PPWO ), 2*NNB,
778: $ A( 1, J ), LDA, WORK( PW ),
779: $ LWORK-PW+1, IERR )
780: ELSE
781: *
782: * Ignore the structure of U.
783: *
784: CALL ZGEMM( 'No Transpose', 'No Transpose', TOP,
785: $ 2*NNB, 2*NNB, CONE, A( 1, J ), LDA,
786: $ WORK( PPWO ), 2*NNB, CZERO,
787: $ WORK( PW ), TOP )
788: CALL ZLACPY( 'All', TOP, 2*NNB, WORK( PW ), TOP,
789: $ A( 1, J ), LDA )
790: END IF
791: PPWO = PPWO + 4*NNB*NNB
792: END DO
793: *
794: J = IHI - NBLST + 1
795: CALL ZGEMM( 'No Transpose', 'No Transpose', TOP,
796: $ NBLST, NBLST, CONE, B( 1, J ), LDB,
797: $ WORK, NBLST, CZERO, WORK( PW ), TOP )
798: CALL ZLACPY( 'All', TOP, NBLST, WORK( PW ), TOP,
799: $ B( 1, J ), LDB )
800: PPWO = NBLST*NBLST + 1
801: J0 = J - NNB
802: DO J = J0, JCOL+1, -NNB
803: IF ( BLK22 ) THEN
804: *
805: * Exploit the structure of U.
806: *
807: CALL ZUNM22( 'Right', 'No Transpose', TOP, 2*NNB,
808: $ NNB, NNB, WORK( PPWO ), 2*NNB,
809: $ B( 1, J ), LDB, WORK( PW ),
810: $ LWORK-PW+1, IERR )
811: ELSE
812: *
813: * Ignore the structure of U.
814: *
815: CALL ZGEMM( 'No Transpose', 'No Transpose', TOP,
816: $ 2*NNB, 2*NNB, CONE, B( 1, J ), LDB,
817: $ WORK( PPWO ), 2*NNB, CZERO,
818: $ WORK( PW ), TOP )
819: CALL ZLACPY( 'All', TOP, 2*NNB, WORK( PW ), TOP,
820: $ B( 1, J ), LDB )
821: END IF
822: PPWO = PPWO + 4*NNB*NNB
823: END DO
824: END IF
825: *
826: * Apply accumulated unitary matrices to Z.
827: *
828: IF( WANTZ ) THEN
829: J = IHI - NBLST + 1
830: IF ( INITQ ) THEN
831: TOPQ = MAX( 2, J - JCOL + 1 )
832: NH = IHI - TOPQ + 1
833: ELSE
834: TOPQ = 1
835: NH = N
836: END IF
837: CALL ZGEMM( 'No Transpose', 'No Transpose', NH,
838: $ NBLST, NBLST, CONE, Z( TOPQ, J ), LDZ,
839: $ WORK, NBLST, CZERO, WORK( PW ), NH )
840: CALL ZLACPY( 'All', NH, NBLST, WORK( PW ), NH,
841: $ Z( TOPQ, J ), LDZ )
842: PPWO = NBLST*NBLST + 1
843: J0 = J - NNB
844: DO J = J0, JCOL+1, -NNB
845: IF ( INITQ ) THEN
846: TOPQ = MAX( 2, J - JCOL + 1 )
847: NH = IHI - TOPQ + 1
848: END IF
849: IF ( BLK22 ) THEN
850: *
851: * Exploit the structure of U.
852: *
853: CALL ZUNM22( 'Right', 'No Transpose', NH, 2*NNB,
854: $ NNB, NNB, WORK( PPWO ), 2*NNB,
855: $ Z( TOPQ, J ), LDZ, WORK( PW ),
856: $ LWORK-PW+1, IERR )
857: ELSE
858: *
859: * Ignore the structure of U.
860: *
861: CALL ZGEMM( 'No Transpose', 'No Transpose', NH,
862: $ 2*NNB, 2*NNB, CONE, Z( TOPQ, J ), LDZ,
863: $ WORK( PPWO ), 2*NNB, CZERO, WORK( PW ),
864: $ NH )
865: CALL ZLACPY( 'All', NH, 2*NNB, WORK( PW ), NH,
866: $ Z( TOPQ, J ), LDZ )
867: END IF
868: PPWO = PPWO + 4*NNB*NNB
869: END DO
870: END IF
871: END DO
872: END IF
873: *
874: * Use unblocked code to reduce the rest of the matrix
875: * Avoid re-initialization of modified Q and Z.
876: *
877: COMPQ2 = COMPQ
878: COMPZ2 = COMPZ
879: IF ( JCOL.NE.ILO ) THEN
880: IF ( WANTQ )
881: $ COMPQ2 = 'V'
882: IF ( WANTZ )
883: $ COMPZ2 = 'V'
884: END IF
885: *
886: IF ( JCOL.LT.IHI )
887: $ CALL ZGGHRD( COMPQ2, COMPZ2, N, JCOL, IHI, A, LDA, B, LDB, Q,
888: $ LDQ, Z, LDZ, IERR )
889: WORK( 1 ) = DCMPLX( LWKOPT )
890: *
891: RETURN
892: *
893: * End of ZGGHD3
894: *
895: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgghd3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack