LAPACK  3.9.0
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.8.0) --
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, cgemm,
275  $ cgemv, ctrmv, clacpy, xerbla
276 * ..
277 * .. Intrinsic Functions ..
278  INTRINSIC real, cmplx, conjg, max
279 * ..
280 * .. Executable Statements ..
281 *
282 * Decode and test the input parameters.
283 *
284  info = 0
285  nb = ilaenv( 1, 'CGGHD3', ' ', n, ilo, ihi, -1 )
286  lwkopt = max( 6*n*nb, 1 )
287  work( 1 ) = cmplx( lwkopt )
288  initq = lsame( compq, 'I' )
289  wantq = initq .OR. lsame( compq, 'V' )
290  initz = lsame( compz, 'I' )
291  wantz = initz .OR. lsame( compz, 'V' )
292  lquery = ( lwork.EQ.-1 )
293 *
294  IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
295  info = -1
296  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
297  info = -2
298  ELSE IF( n.LT.0 ) THEN
299  info = -3
300  ELSE IF( ilo.LT.1 ) THEN
301  info = -4
302  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
303  info = -5
304  ELSE IF( lda.LT.max( 1, n ) ) THEN
305  info = -7
306  ELSE IF( ldb.LT.max( 1, n ) ) THEN
307  info = -9
308  ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 ) THEN
309  info = -11
310  ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 ) THEN
311  info = -13
312  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
313  info = -15
314  END IF
315  IF( info.NE.0 ) THEN
316  CALL xerbla( 'CGGHD3', -info )
317  RETURN
318  ELSE IF( lquery ) THEN
319  RETURN
320  END IF
321 *
322 * Initialize Q and Z if desired.
323 *
324  IF( initq )
325  $ CALL claset( 'All', n, n, czero, cone, q, ldq )
326  IF( initz )
327  $ CALL claset( 'All', n, n, czero, cone, z, ldz )
328 *
329 * Zero out lower triangle of B.
330 *
331  IF( n.GT.1 )
332  $ CALL claset( 'Lower', n-1, n-1, czero, czero, b(2, 1), ldb )
333 *
334 * Quick return if possible
335 *
336  nh = ihi - ilo + 1
337  IF( nh.LE.1 ) THEN
338  work( 1 ) = cone
339  RETURN
340  END IF
341 *
342 * Determine the blocksize.
343 *
344  nbmin = ilaenv( 2, 'CGGHD3', ' ', n, ilo, ihi, -1 )
345  IF( nb.GT.1 .AND. nb.LT.nh ) THEN
346 *
347 * Determine when to use unblocked instead of blocked code.
348 *
349  nx = max( nb, ilaenv( 3, 'CGGHD3', ' ', n, ilo, ihi, -1 ) )
350  IF( nx.LT.nh ) THEN
351 *
352 * Determine if workspace is large enough for blocked code.
353 *
354  IF( lwork.LT.lwkopt ) THEN
355 *
356 * Not enough workspace to use optimal NB: determine the
357 * minimum value of NB, and reduce NB or force use of
358 * unblocked code.
359 *
360  nbmin = max( 2, ilaenv( 2, 'CGGHD3', ' ', n, ilo, ihi,
361  $ -1 ) )
362  IF( lwork.GE.6*n*nbmin ) THEN
363  nb = lwork / ( 6*n )
364  ELSE
365  nb = 1
366  END IF
367  END IF
368  END IF
369  END IF
370 *
371  IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
372 *
373 * Use unblocked code below
374 *
375  jcol = ilo
376 *
377  ELSE
378 *
379 * Use blocked code
380 *
381  kacc22 = ilaenv( 16, 'CGGHD3', ' ', n, ilo, ihi, -1 )
382  blk22 = kacc22.EQ.2
383  DO jcol = ilo, ihi-2, nb
384  nnb = min( nb, ihi-jcol-1 )
385 *
386 * Initialize small unitary factors that will hold the
387 * accumulated Givens rotations in workspace.
388 * N2NB denotes the number of 2*NNB-by-2*NNB factors
389 * NBLST denotes the (possibly smaller) order of the last
390 * factor.
391 *
392  n2nb = ( ihi-jcol-1 ) / nnb - 1
393  nblst = ihi - jcol - n2nb*nnb
394  CALL claset( 'All', nblst, nblst, czero, cone, work, nblst )
395  pw = nblst * nblst + 1
396  DO i = 1, n2nb
397  CALL claset( 'All', 2*nnb, 2*nnb, czero, cone,
398  $ work( pw ), 2*nnb )
399  pw = pw + 4*nnb*nnb
400  END DO
401 *
402 * Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
403 *
404  DO j = jcol, jcol+nnb-1
405 *
406 * Reduce Jth column of A. Store cosines and sines in Jth
407 * column of A and B, respectively.
408 *
409  DO i = ihi, j+2, -1
410  temp = a( i-1, j )
411  CALL clartg( temp, a( i, j ), c, s, a( i-1, j ) )
412  a( i, j ) = cmplx( c )
413  b( i, j ) = s
414  END DO
415 *
416 * Accumulate Givens rotations into workspace array.
417 *
418  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
419  len = 2 + j - jcol
420  jrow = j + n2nb*nnb + 2
421  DO i = ihi, jrow, -1
422  ctemp = a( i, j )
423  s = b( i, j )
424  DO jj = ppw, ppw+len-1
425  temp = work( jj + nblst )
426  work( jj + nblst ) = ctemp*temp - s*work( jj )
427  work( jj ) = conjg( s )*temp + ctemp*work( jj )
428  END DO
429  len = len + 1
430  ppw = ppw - nblst - 1
431  END DO
432 *
433  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
434  j0 = jrow - nnb
435  DO jrow = j0, j+2, -nnb
436  ppw = ppwo
437  len = 2 + j - jcol
438  DO i = jrow+nnb-1, jrow, -1
439  ctemp = a( i, j )
440  s = b( i, j )
441  DO jj = ppw, ppw+len-1
442  temp = work( jj + 2*nnb )
443  work( jj + 2*nnb ) = ctemp*temp - s*work( jj )
444  work( jj ) = conjg( s )*temp + ctemp*work( jj )
445  END DO
446  len = len + 1
447  ppw = ppw - 2*nnb - 1
448  END DO
449  ppwo = ppwo + 4*nnb*nnb
450  END DO
451 *
452 * TOP denotes the number of top rows in A and B that will
453 * not be updated during the next steps.
454 *
455  IF( jcol.LE.2 ) THEN
456  top = 0
457  ELSE
458  top = jcol
459  END IF
460 *
461 * Propagate transformations through B and replace stored
462 * left sines/cosines by right sines/cosines.
463 *
464  DO jj = n, j+1, -1
465 *
466 * Update JJth column of B.
467 *
468  DO i = min( jj+1, ihi ), j+2, -1
469  ctemp = a( i, j )
470  s = b( i, j )
471  temp = b( i, jj )
472  b( i, jj ) = ctemp*temp - conjg( s )*b( i-1, jj )
473  b( i-1, jj ) = s*temp + ctemp*b( i-1, jj )
474  END DO
475 *
476 * Annihilate B( JJ+1, JJ ).
477 *
478  IF( jj.LT.ihi ) THEN
479  temp = b( jj+1, jj+1 )
480  CALL clartg( temp, b( jj+1, jj ), c, s,
481  $ b( jj+1, jj+1 ) )
482  b( jj+1, jj ) = czero
483  CALL crot( jj-top, b( top+1, jj+1 ), 1,
484  $ b( top+1, jj ), 1, c, s )
485  a( jj+1, j ) = cmplx( c )
486  b( jj+1, j ) = -conjg( s )
487  END IF
488  END DO
489 *
490 * Update A by transformations from right.
491 *
492  jj = mod( ihi-j-1, 3 )
493  DO i = ihi-j-3, jj+1, -3
494  ctemp = 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 + conjg( s2 )*temp2
507  temp2 = -s2*temp3 + c2*temp2
508  a( k, j+i+2 ) = c1*temp2 + conjg( s1 )*temp1
509  temp1 = -s1*temp2 + c1*temp1
510  a( k, j+i+1 ) = ctemp*temp1 + conjg( s )*temp
511  a( k, j+i ) = -s*temp1 + ctemp*temp
512  END DO
513  END DO
514 *
515  IF( jj.GT.0 ) THEN
516  DO i = jj, 1, -1
517  c = dble( a( j+1+i, j ) )
518  CALL crot( ihi-top, a( top+1, j+i+1 ), 1,
519  $ a( top+1, j+i ), 1, c,
520  $ -conjg( b( j+1+i, j ) ) )
521  END DO
522  END IF
523 *
524 * Update (J+1)th column of A by transformations from left.
525 *
526  IF ( j .LT. jcol + nnb - 1 ) THEN
527  len = 1 + j - jcol
528 *
529 * Multiply with the trailing accumulated unitary
530 * matrix, which takes the form
531 *
532 * [ U11 U12 ]
533 * U = [ ],
534 * [ U21 U22 ]
535 *
536 * where U21 is a LEN-by-LEN matrix and U12 is lower
537 * triangular.
538 *
539  jrow = ihi - nblst + 1
540  CALL cgemv( 'Conjugate', nblst, len, cone, work,
541  $ nblst, a( jrow, j+1 ), 1, czero,
542  $ work( pw ), 1 )
543  ppw = pw + len
544  DO i = jrow, jrow+nblst-len-1
545  work( ppw ) = a( i, j+1 )
546  ppw = ppw + 1
547  END DO
548  CALL ctrmv( 'Lower', 'Conjugate', 'Non-unit',
549  $ nblst-len, work( len*nblst + 1 ), nblst,
550  $ work( pw+len ), 1 )
551  CALL cgemv( 'Conjugate', len, nblst-len, cone,
552  $ work( (len+1)*nblst - len + 1 ), nblst,
553  $ a( jrow+nblst-len, j+1 ), 1, cone,
554  $ work( pw+len ), 1 )
555  ppw = pw
556  DO i = jrow, jrow+nblst-1
557  a( i, j+1 ) = work( ppw )
558  ppw = ppw + 1
559  END DO
560 *
561 * Multiply with the other accumulated unitary
562 * matrices, which take the form
563 *
564 * [ U11 U12 0 ]
565 * [ ]
566 * U = [ U21 U22 0 ],
567 * [ ]
568 * [ 0 0 I ]
569 *
570 * where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
571 * matrix, U21 is a LEN-by-LEN upper triangular matrix
572 * and U12 is an NNB-by-NNB lower triangular matrix.
573 *
574  ppwo = 1 + nblst*nblst
575  j0 = jrow - nnb
576  DO jrow = j0, jcol+1, -nnb
577  ppw = pw + len
578  DO i = jrow, jrow+nnb-1
579  work( ppw ) = a( i, j+1 )
580  ppw = ppw + 1
581  END DO
582  ppw = pw
583  DO i = jrow+nnb, jrow+nnb+len-1
584  work( ppw ) = a( i, j+1 )
585  ppw = ppw + 1
586  END DO
587  CALL ctrmv( 'Upper', 'Conjugate', 'Non-unit', len,
588  $ work( ppwo + nnb ), 2*nnb, work( pw ),
589  $ 1 )
590  CALL ctrmv( 'Lower', 'Conjugate', 'Non-unit', nnb,
591  $ work( ppwo + 2*len*nnb ),
592  $ 2*nnb, work( pw + len ), 1 )
593  CALL cgemv( 'Conjugate', nnb, len, cone,
594  $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
595  $ cone, work( pw ), 1 )
596  CALL cgemv( 'Conjugate', len, nnb, cone,
597  $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
598  $ a( jrow+nnb, j+1 ), 1, cone,
599  $ work( pw+len ), 1 )
600  ppw = pw
601  DO i = jrow, jrow+len+nnb-1
602  a( i, j+1 ) = work( ppw )
603  ppw = ppw + 1
604  END DO
605  ppwo = ppwo + 4*nnb*nnb
606  END DO
607  END IF
608  END DO
609 *
610 * Apply accumulated unitary matrices to A.
611 *
612  cola = n - jcol - nnb + 1
613  j = ihi - nblst + 1
614  CALL cgemm( 'Conjugate', 'No Transpose', nblst,
615  $ cola, nblst, cone, work, nblst,
616  $ a( j, jcol+nnb ), lda, czero, work( pw ),
617  $ nblst )
618  CALL clacpy( 'All', nblst, cola, work( pw ), nblst,
619  $ a( j, jcol+nnb ), lda )
620  ppwo = nblst*nblst + 1
621  j0 = j - nnb
622  DO j = j0, jcol+1, -nnb
623  IF ( blk22 ) THEN
624 *
625 * Exploit the structure of
626 *
627 * [ U11 U12 ]
628 * U = [ ]
629 * [ U21 U22 ],
630 *
631 * where all blocks are NNB-by-NNB, U21 is upper
632 * triangular and U12 is lower triangular.
633 *
634  CALL cunm22( 'Left', 'Conjugate', 2*nnb, cola, nnb,
635  $ nnb, work( ppwo ), 2*nnb,
636  $ a( j, jcol+nnb ), lda, work( pw ),
637  $ lwork-pw+1, ierr )
638  ELSE
639 *
640 * Ignore the structure of U.
641 *
642  CALL cgemm( 'Conjugate', 'No Transpose', 2*nnb,
643  $ cola, 2*nnb, cone, work( ppwo ), 2*nnb,
644  $ a( j, jcol+nnb ), lda, czero, work( pw ),
645  $ 2*nnb )
646  CALL clacpy( 'All', 2*nnb, cola, work( pw ), 2*nnb,
647  $ a( j, jcol+nnb ), lda )
648  END IF
649  ppwo = ppwo + 4*nnb*nnb
650  END DO
651 *
652 * Apply accumulated unitary matrices to Q.
653 *
654  IF( wantq ) THEN
655  j = ihi - nblst + 1
656  IF ( initq ) THEN
657  topq = max( 2, j - jcol + 1 )
658  nh = ihi - topq + 1
659  ELSE
660  topq = 1
661  nh = n
662  END IF
663  CALL cgemm( 'No Transpose', 'No Transpose', nh,
664  $ nblst, nblst, cone, q( topq, j ), ldq,
665  $ work, nblst, czero, work( pw ), nh )
666  CALL clacpy( 'All', nh, nblst, work( pw ), nh,
667  $ q( topq, j ), ldq )
668  ppwo = nblst*nblst + 1
669  j0 = j - nnb
670  DO j = j0, jcol+1, -nnb
671  IF ( initq ) THEN
672  topq = max( 2, j - jcol + 1 )
673  nh = ihi - topq + 1
674  END IF
675  IF ( blk22 ) THEN
676 *
677 * Exploit the structure of U.
678 *
679  CALL cunm22( 'Right', 'No Transpose', nh, 2*nnb,
680  $ nnb, nnb, work( ppwo ), 2*nnb,
681  $ q( topq, j ), ldq, work( pw ),
682  $ lwork-pw+1, ierr )
683  ELSE
684 *
685 * Ignore the structure of U.
686 *
687  CALL cgemm( 'No Transpose', 'No Transpose', nh,
688  $ 2*nnb, 2*nnb, cone, q( topq, j ), ldq,
689  $ work( ppwo ), 2*nnb, czero, work( pw ),
690  $ nh )
691  CALL clacpy( 'All', nh, 2*nnb, work( pw ), nh,
692  $ q( topq, j ), ldq )
693  END IF
694  ppwo = ppwo + 4*nnb*nnb
695  END DO
696  END IF
697 *
698 * Accumulate right Givens rotations if required.
699 *
700  IF ( wantz .OR. top.GT.0 ) THEN
701 *
702 * Initialize small unitary factors that will hold the
703 * accumulated Givens rotations in workspace.
704 *
705  CALL claset( 'All', nblst, nblst, czero, cone, work,
706  $ nblst )
707  pw = nblst * nblst + 1
708  DO i = 1, n2nb
709  CALL claset( 'All', 2*nnb, 2*nnb, czero, cone,
710  $ work( pw ), 2*nnb )
711  pw = pw + 4*nnb*nnb
712  END DO
713 *
714 * Accumulate Givens rotations into workspace array.
715 *
716  DO j = jcol, jcol+nnb-1
717  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
718  len = 2 + j - jcol
719  jrow = j + n2nb*nnb + 2
720  DO i = ihi, jrow, -1
721  ctemp = a( i, j )
722  a( i, j ) = czero
723  s = b( i, j )
724  b( i, j ) = czero
725  DO jj = ppw, ppw+len-1
726  temp = work( jj + nblst )
727  work( jj + nblst ) = ctemp*temp -
728  $ conjg( s )*work( jj )
729  work( jj ) = s*temp + ctemp*work( jj )
730  END DO
731  len = len + 1
732  ppw = ppw - nblst - 1
733  END DO
734 *
735  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
736  j0 = jrow - nnb
737  DO jrow = j0, j+2, -nnb
738  ppw = ppwo
739  len = 2 + j - jcol
740  DO i = jrow+nnb-1, jrow, -1
741  ctemp = a( i, j )
742  a( i, j ) = czero
743  s = b( i, j )
744  b( i, j ) = czero
745  DO jj = ppw, ppw+len-1
746  temp = work( jj + 2*nnb )
747  work( jj + 2*nnb ) = ctemp*temp -
748  $ conjg( s )*work( jj )
749  work( jj ) = s*temp + ctemp*work( jj )
750  END DO
751  len = len + 1
752  ppw = ppw - 2*nnb - 1
753  END DO
754  ppwo = ppwo + 4*nnb*nnb
755  END DO
756  END DO
757  ELSE
758 *
759  CALL claset( 'Lower', ihi - jcol - 1, nnb, czero, czero,
760  $ a( jcol + 2, jcol ), lda )
761  CALL claset( 'Lower', ihi - jcol - 1, nnb, czero, czero,
762  $ b( jcol + 2, jcol ), ldb )
763  END IF
764 *
765 * Apply accumulated unitary matrices to A and B.
766 *
767  IF ( top.GT.0 ) THEN
768  j = ihi - nblst + 1
769  CALL cgemm( 'No Transpose', 'No Transpose', top,
770  $ nblst, nblst, cone, a( 1, j ), lda,
771  $ work, nblst, czero, work( pw ), top )
772  CALL clacpy( 'All', top, nblst, work( pw ), top,
773  $ a( 1, j ), lda )
774  ppwo = nblst*nblst + 1
775  j0 = j - nnb
776  DO j = j0, jcol+1, -nnb
777  IF ( blk22 ) THEN
778 *
779 * Exploit the structure of U.
780 *
781  CALL cunm22( 'Right', 'No Transpose', top, 2*nnb,
782  $ nnb, nnb, work( ppwo ), 2*nnb,
783  $ a( 1, j ), lda, work( pw ),
784  $ lwork-pw+1, ierr )
785  ELSE
786 *
787 * Ignore the structure of U.
788 *
789  CALL cgemm( 'No Transpose', 'No Transpose', top,
790  $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
791  $ work( ppwo ), 2*nnb, czero,
792  $ work( pw ), top )
793  CALL clacpy( 'All', top, 2*nnb, work( pw ), top,
794  $ a( 1, j ), lda )
795  END IF
796  ppwo = ppwo + 4*nnb*nnb
797  END DO
798 *
799  j = ihi - nblst + 1
800  CALL cgemm( 'No Transpose', 'No Transpose', top,
801  $ nblst, nblst, cone, b( 1, j ), ldb,
802  $ work, nblst, czero, work( pw ), top )
803  CALL clacpy( 'All', top, nblst, work( pw ), top,
804  $ b( 1, j ), ldb )
805  ppwo = nblst*nblst + 1
806  j0 = j - nnb
807  DO j = j0, jcol+1, -nnb
808  IF ( blk22 ) THEN
809 *
810 * Exploit the structure of U.
811 *
812  CALL cunm22( 'Right', 'No Transpose', top, 2*nnb,
813  $ nnb, nnb, work( ppwo ), 2*nnb,
814  $ b( 1, j ), ldb, work( pw ),
815  $ lwork-pw+1, ierr )
816  ELSE
817 *
818 * Ignore the structure of U.
819 *
820  CALL cgemm( 'No Transpose', 'No Transpose', top,
821  $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
822  $ work( ppwo ), 2*nnb, czero,
823  $ work( pw ), top )
824  CALL clacpy( 'All', top, 2*nnb, work( pw ), top,
825  $ b( 1, j ), ldb )
826  END IF
827  ppwo = ppwo + 4*nnb*nnb
828  END DO
829  END IF
830 *
831 * Apply accumulated unitary matrices to Z.
832 *
833  IF( wantz ) THEN
834  j = ihi - nblst + 1
835  IF ( initq ) THEN
836  topq = max( 2, j - jcol + 1 )
837  nh = ihi - topq + 1
838  ELSE
839  topq = 1
840  nh = n
841  END IF
842  CALL cgemm( 'No Transpose', 'No Transpose', nh,
843  $ nblst, nblst, cone, z( topq, j ), ldz,
844  $ work, nblst, czero, work( pw ), nh )
845  CALL clacpy( 'All', nh, nblst, work( pw ), nh,
846  $ z( topq, j ), ldz )
847  ppwo = nblst*nblst + 1
848  j0 = j - nnb
849  DO j = j0, jcol+1, -nnb
850  IF ( initq ) THEN
851  topq = max( 2, j - jcol + 1 )
852  nh = ihi - topq + 1
853  END IF
854  IF ( blk22 ) THEN
855 *
856 * Exploit the structure of U.
857 *
858  CALL cunm22( 'Right', 'No Transpose', nh, 2*nnb,
859  $ nnb, nnb, work( ppwo ), 2*nnb,
860  $ z( topq, j ), ldz, work( pw ),
861  $ lwork-pw+1, ierr )
862  ELSE
863 *
864 * Ignore the structure of U.
865 *
866  CALL cgemm( 'No Transpose', 'No Transpose', nh,
867  $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
868  $ work( ppwo ), 2*nnb, czero, work( pw ),
869  $ nh )
870  CALL clacpy( 'All', nh, 2*nnb, work( pw ), nh,
871  $ z( topq, j ), ldz )
872  END IF
873  ppwo = ppwo + 4*nnb*nnb
874  END DO
875  END IF
876  END DO
877  END IF
878 *
879 * Use unblocked code to reduce the rest of the matrix
880 * Avoid re-initialization of modified Q and Z.
881 *
882  compq2 = compq
883  compz2 = compz
884  IF ( jcol.NE.ilo ) THEN
885  IF ( wantq )
886  $ compq2 = 'V'
887  IF ( wantz )
888  $ compz2 = 'V'
889  END IF
890 *
891  IF ( jcol.LT.ihi )
892  $ CALL cgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
893  $ ldq, z, ldz, ierr )
894  work( 1 ) = cmplx( lwkopt )
895 *
896  RETURN
897 *
898 * End of CGGHD3
899 *
900  END
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cgghd3
subroutine cgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
CGGHD3
Definition: cgghd3.f:233
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
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
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
ctrmv
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:149
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
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
clartg
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f:105