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