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