LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zgghd3.f
Go to the documentation of this file.
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 *> \ingroup complex16OTHERcomputational
212 *
213 *> \par Further Details:
214 * =====================
215 *>
216 *> \verbatim
217 *>
218 *> This routine reduces A to Hessenberg form and maintains B in triangular form
219 *> using a blocked variant of Moler and Stewart's original algorithm,
220 *> as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
221 *> (BIT 2008).
222 *> \endverbatim
223 *>
224 * =====================================================================
225  SUBROUTINE zgghd3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
226  $ LDQ, Z, LDZ, WORK, LWORK, INFO )
227 *
228 * -- LAPACK computational routine --
229 * -- LAPACK is a software package provided by Univ. of Tennessee, --
230 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231 *
232  IMPLICIT NONE
233 *
234 * .. Scalar Arguments ..
235  CHARACTER COMPQ, COMPZ
236  INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
237 * ..
238 * .. Array Arguments ..
239  COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
240  $ z( ldz, * ), work( * )
241 * ..
242 *
243 * =====================================================================
244 *
245 * .. Parameters ..
246  COMPLEX*16 CONE, CZERO
247  parameter( cone = ( 1.0d+0, 0.0d+0 ),
248  $ czero = ( 0.0d+0, 0.0d+0 ) )
249 * ..
250 * .. Local Scalars ..
251  LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
252  CHARACTER*1 COMPQ2, COMPZ2
253  INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
254  $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
255  $ nh, nnb, nx, ppw, ppwo, pw, top, topq
256  DOUBLE PRECISION C
257  COMPLEX*16 C1, C2, CTEMP, S, S1, S2, TEMP, TEMP1, TEMP2,
258  $ temp3
259 * ..
260 * .. External Functions ..
261  LOGICAL LSAME
262  INTEGER ILAENV
263  EXTERNAL ilaenv, lsame
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL zgghrd, zlartg, zlaset, zunm22, zrot, zgemm,
267  $ zgemv, ztrmv, zlacpy, xerbla
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC dble, dcmplx, dconjg, max
271 * ..
272 * .. Executable Statements ..
273 *
274 * Decode and test the input parameters.
275 *
276  info = 0
277  nb = ilaenv( 1, 'ZGGHD3', ' ', n, ilo, ihi, -1 )
278  lwkopt = max( 6*n*nb, 1 )
279  work( 1 ) = dcmplx( lwkopt )
280  initq = lsame( compq, 'I' )
281  wantq = initq .OR. lsame( compq, 'V' )
282  initz = lsame( compz, 'I' )
283  wantz = initz .OR. lsame( compz, 'V' )
284  lquery = ( lwork.EQ.-1 )
285 *
286  IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
287  info = -1
288  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
289  info = -2
290  ELSE IF( n.LT.0 ) THEN
291  info = -3
292  ELSE IF( ilo.LT.1 ) THEN
293  info = -4
294  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
295  info = -5
296  ELSE IF( lda.LT.max( 1, n ) ) THEN
297  info = -7
298  ELSE IF( ldb.LT.max( 1, n ) ) THEN
299  info = -9
300  ELSE IF( ( wantq .AND. ldq.LT.n ) .OR. ldq.LT.1 ) THEN
301  info = -11
302  ELSE IF( ( wantz .AND. ldz.LT.n ) .OR. ldz.LT.1 ) THEN
303  info = -13
304  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
305  info = -15
306  END IF
307  IF( info.NE.0 ) THEN
308  CALL xerbla( 'ZGGHD3', -info )
309  RETURN
310  ELSE IF( lquery ) THEN
311  RETURN
312  END IF
313 *
314 * Initialize Q and Z if desired.
315 *
316  IF( initq )
317  $ CALL zlaset( 'All', n, n, czero, cone, q, ldq )
318  IF( initz )
319  $ CALL zlaset( 'All', n, n, czero, cone, z, ldz )
320 *
321 * Zero out lower triangle of B.
322 *
323  IF( n.GT.1 )
324  $ CALL zlaset( 'Lower', n-1, n-1, czero, czero, b(2, 1), ldb )
325 *
326 * Quick return if possible
327 *
328  nh = ihi - ilo + 1
329  IF( nh.LE.1 ) THEN
330  work( 1 ) = cone
331  RETURN
332  END IF
333 *
334 * Determine the blocksize.
335 *
336  nbmin = ilaenv( 2, 'ZGGHD3', ' ', n, ilo, ihi, -1 )
337  IF( nb.GT.1 .AND. nb.LT.nh ) THEN
338 *
339 * Determine when to use unblocked instead of blocked code.
340 *
341  nx = max( nb, ilaenv( 3, 'ZGGHD3', ' ', n, ilo, ihi, -1 ) )
342  IF( nx.LT.nh ) THEN
343 *
344 * Determine if workspace is large enough for blocked code.
345 *
346  IF( lwork.LT.lwkopt ) THEN
347 *
348 * Not enough workspace to use optimal NB: determine the
349 * minimum value of NB, and reduce NB or force use of
350 * unblocked code.
351 *
352  nbmin = max( 2, ilaenv( 2, 'ZGGHD3', ' ', n, ilo, ihi,
353  $ -1 ) )
354  IF( lwork.GE.6*n*nbmin ) THEN
355  nb = lwork / ( 6*n )
356  ELSE
357  nb = 1
358  END IF
359  END IF
360  END IF
361  END IF
362 *
363  IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
364 *
365 * Use unblocked code below
366 *
367  jcol = ilo
368 *
369  ELSE
370 *
371 * Use blocked code
372 *
373  kacc22 = ilaenv( 16, 'ZGGHD3', ' ', n, ilo, ihi, -1 )
374  blk22 = kacc22.EQ.2
375  DO jcol = ilo, ihi-2, nb
376  nnb = min( nb, ihi-jcol-1 )
377 *
378 * Initialize small unitary factors that will hold the
379 * accumulated Givens rotations in workspace.
380 * N2NB denotes the number of 2*NNB-by-2*NNB factors
381 * NBLST denotes the (possibly smaller) order of the last
382 * factor.
383 *
384  n2nb = ( ihi-jcol-1 ) / nnb - 1
385  nblst = ihi - jcol - n2nb*nnb
386  CALL zlaset( 'All', nblst, nblst, czero, cone, work, nblst )
387  pw = nblst * nblst + 1
388  DO i = 1, n2nb
389  CALL zlaset( 'All', 2*nnb, 2*nnb, czero, cone,
390  $ work( pw ), 2*nnb )
391  pw = pw + 4*nnb*nnb
392  END DO
393 *
394 * Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
395 *
396  DO j = jcol, jcol+nnb-1
397 *
398 * Reduce Jth column of A. Store cosines and sines in Jth
399 * column of A and B, respectively.
400 *
401  DO i = ihi, j+2, -1
402  temp = a( i-1, j )
403  CALL zlartg( temp, a( i, j ), c, s, a( i-1, j ) )
404  a( i, j ) = dcmplx( c )
405  b( i, j ) = s
406  END DO
407 *
408 * Accumulate Givens rotations into workspace array.
409 *
410  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
411  len = 2 + j - jcol
412  jrow = j + n2nb*nnb + 2
413  DO i = ihi, jrow, -1
414  ctemp = a( i, j )
415  s = b( i, j )
416  DO jj = ppw, ppw+len-1
417  temp = work( jj + nblst )
418  work( jj + nblst ) = ctemp*temp - s*work( jj )
419  work( jj ) = dconjg( s )*temp + ctemp*work( jj )
420  END DO
421  len = len + 1
422  ppw = ppw - nblst - 1
423  END DO
424 *
425  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
426  j0 = jrow - nnb
427  DO jrow = j0, j+2, -nnb
428  ppw = ppwo
429  len = 2 + j - jcol
430  DO i = jrow+nnb-1, jrow, -1
431  ctemp = a( i, j )
432  s = b( i, j )
433  DO jj = ppw, ppw+len-1
434  temp = work( jj + 2*nnb )
435  work( jj + 2*nnb ) = ctemp*temp - s*work( jj )
436  work( jj ) = dconjg( s )*temp + ctemp*work( jj )
437  END DO
438  len = len + 1
439  ppw = ppw - 2*nnb - 1
440  END DO
441  ppwo = ppwo + 4*nnb*nnb
442  END DO
443 *
444 * TOP denotes the number of top rows in A and B that will
445 * not be updated during the next steps.
446 *
447  IF( jcol.LE.2 ) THEN
448  top = 0
449  ELSE
450  top = jcol
451  END IF
452 *
453 * Propagate transformations through B and replace stored
454 * left sines/cosines by right sines/cosines.
455 *
456  DO jj = n, j+1, -1
457 *
458 * Update JJth column of B.
459 *
460  DO i = min( jj+1, ihi ), j+2, -1
461  ctemp = a( i, j )
462  s = b( i, j )
463  temp = b( i, jj )
464  b( i, jj ) = ctemp*temp - dconjg( s )*b( i-1, jj )
465  b( i-1, jj ) = s*temp + ctemp*b( i-1, jj )
466  END DO
467 *
468 * Annihilate B( JJ+1, JJ ).
469 *
470  IF( jj.LT.ihi ) THEN
471  temp = b( jj+1, jj+1 )
472  CALL zlartg( temp, b( jj+1, jj ), c, s,
473  $ b( jj+1, jj+1 ) )
474  b( jj+1, jj ) = czero
475  CALL zrot( jj-top, b( top+1, jj+1 ), 1,
476  $ b( top+1, jj ), 1, c, s )
477  a( jj+1, j ) = dcmplx( c )
478  b( jj+1, j ) = -dconjg( s )
479  END IF
480  END DO
481 *
482 * Update A by transformations from right.
483 *
484  jj = mod( ihi-j-1, 3 )
485  DO i = ihi-j-3, jj+1, -3
486  ctemp = a( j+1+i, j )
487  s = -b( j+1+i, j )
488  c1 = a( j+2+i, j )
489  s1 = -b( j+2+i, j )
490  c2 = a( j+3+i, j )
491  s2 = -b( j+3+i, j )
492 *
493  DO k = top+1, ihi
494  temp = a( k, j+i )
495  temp1 = a( k, j+i+1 )
496  temp2 = a( k, j+i+2 )
497  temp3 = a( k, j+i+3 )
498  a( k, j+i+3 ) = c2*temp3 + dconjg( s2 )*temp2
499  temp2 = -s2*temp3 + c2*temp2
500  a( k, j+i+2 ) = c1*temp2 + dconjg( s1 )*temp1
501  temp1 = -s1*temp2 + c1*temp1
502  a( k, j+i+1 ) = ctemp*temp1 + dconjg( s )*temp
503  a( k, j+i ) = -s*temp1 + ctemp*temp
504  END DO
505  END DO
506 *
507  IF( jj.GT.0 ) THEN
508  DO i = jj, 1, -1
509  c = dble( a( j+1+i, j ) )
510  CALL zrot( ihi-top, a( top+1, j+i+1 ), 1,
511  $ a( top+1, j+i ), 1, c,
512  $ -dconjg( b( j+1+i, j ) ) )
513  END DO
514  END IF
515 *
516 * Update (J+1)th column of A by transformations from left.
517 *
518  IF ( j .LT. jcol + nnb - 1 ) THEN
519  len = 1 + j - jcol
520 *
521 * Multiply with the trailing accumulated unitary
522 * matrix, which takes the form
523 *
524 * [ U11 U12 ]
525 * U = [ ],
526 * [ U21 U22 ]
527 *
528 * where U21 is a LEN-by-LEN matrix and U12 is lower
529 * triangular.
530 *
531  jrow = ihi - nblst + 1
532  CALL zgemv( 'Conjugate', nblst, len, cone, work,
533  $ nblst, a( jrow, j+1 ), 1, czero,
534  $ work( pw ), 1 )
535  ppw = pw + len
536  DO i = jrow, jrow+nblst-len-1
537  work( ppw ) = a( i, j+1 )
538  ppw = ppw + 1
539  END DO
540  CALL ztrmv( 'Lower', 'Conjugate', 'Non-unit',
541  $ nblst-len, work( len*nblst + 1 ), nblst,
542  $ work( pw+len ), 1 )
543  CALL zgemv( 'Conjugate', len, nblst-len, cone,
544  $ work( (len+1)*nblst - len + 1 ), nblst,
545  $ a( jrow+nblst-len, j+1 ), 1, cone,
546  $ work( pw+len ), 1 )
547  ppw = pw
548  DO i = jrow, jrow+nblst-1
549  a( i, j+1 ) = work( ppw )
550  ppw = ppw + 1
551  END DO
552 *
553 * Multiply with the other accumulated unitary
554 * matrices, which take the form
555 *
556 * [ U11 U12 0 ]
557 * [ ]
558 * U = [ U21 U22 0 ],
559 * [ ]
560 * [ 0 0 I ]
561 *
562 * where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
563 * matrix, U21 is a LEN-by-LEN upper triangular matrix
564 * and U12 is an NNB-by-NNB lower triangular matrix.
565 *
566  ppwo = 1 + nblst*nblst
567  j0 = jrow - nnb
568  DO jrow = j0, jcol+1, -nnb
569  ppw = pw + len
570  DO i = jrow, jrow+nnb-1
571  work( ppw ) = a( i, j+1 )
572  ppw = ppw + 1
573  END DO
574  ppw = pw
575  DO i = jrow+nnb, jrow+nnb+len-1
576  work( ppw ) = a( i, j+1 )
577  ppw = ppw + 1
578  END DO
579  CALL ztrmv( 'Upper', 'Conjugate', 'Non-unit', len,
580  $ work( ppwo + nnb ), 2*nnb, work( pw ),
581  $ 1 )
582  CALL ztrmv( 'Lower', 'Conjugate', 'Non-unit', nnb,
583  $ work( ppwo + 2*len*nnb ),
584  $ 2*nnb, work( pw + len ), 1 )
585  CALL zgemv( 'Conjugate', nnb, len, cone,
586  $ work( ppwo ), 2*nnb, a( jrow, j+1 ), 1,
587  $ cone, work( pw ), 1 )
588  CALL zgemv( 'Conjugate', len, nnb, cone,
589  $ work( ppwo + 2*len*nnb + nnb ), 2*nnb,
590  $ a( jrow+nnb, j+1 ), 1, cone,
591  $ work( pw+len ), 1 )
592  ppw = pw
593  DO i = jrow, jrow+len+nnb-1
594  a( i, j+1 ) = work( ppw )
595  ppw = ppw + 1
596  END DO
597  ppwo = ppwo + 4*nnb*nnb
598  END DO
599  END IF
600  END DO
601 *
602 * Apply accumulated unitary matrices to A.
603 *
604  cola = n - jcol - nnb + 1
605  j = ihi - nblst + 1
606  CALL zgemm( 'Conjugate', 'No Transpose', nblst,
607  $ cola, nblst, cone, work, nblst,
608  $ a( j, jcol+nnb ), lda, czero, work( pw ),
609  $ nblst )
610  CALL zlacpy( 'All', nblst, cola, work( pw ), nblst,
611  $ a( j, jcol+nnb ), lda )
612  ppwo = nblst*nblst + 1
613  j0 = j - nnb
614  DO j = j0, jcol+1, -nnb
615  IF ( blk22 ) THEN
616 *
617 * Exploit the structure of
618 *
619 * [ U11 U12 ]
620 * U = [ ]
621 * [ U21 U22 ],
622 *
623 * where all blocks are NNB-by-NNB, U21 is upper
624 * triangular and U12 is lower triangular.
625 *
626  CALL zunm22( 'Left', 'Conjugate', 2*nnb, cola, nnb,
627  $ nnb, work( ppwo ), 2*nnb,
628  $ a( j, jcol+nnb ), lda, work( pw ),
629  $ lwork-pw+1, ierr )
630  ELSE
631 *
632 * Ignore the structure of U.
633 *
634  CALL zgemm( 'Conjugate', 'No Transpose', 2*nnb,
635  $ cola, 2*nnb, cone, work( ppwo ), 2*nnb,
636  $ a( j, jcol+nnb ), lda, czero, work( pw ),
637  $ 2*nnb )
638  CALL zlacpy( 'All', 2*nnb, cola, work( pw ), 2*nnb,
639  $ a( j, jcol+nnb ), lda )
640  END IF
641  ppwo = ppwo + 4*nnb*nnb
642  END DO
643 *
644 * Apply accumulated unitary matrices to Q.
645 *
646  IF( wantq ) THEN
647  j = ihi - nblst + 1
648  IF ( initq ) THEN
649  topq = max( 2, j - jcol + 1 )
650  nh = ihi - topq + 1
651  ELSE
652  topq = 1
653  nh = n
654  END IF
655  CALL zgemm( 'No Transpose', 'No Transpose', nh,
656  $ nblst, nblst, cone, q( topq, j ), ldq,
657  $ work, nblst, czero, work( pw ), nh )
658  CALL zlacpy( 'All', nh, nblst, work( pw ), nh,
659  $ q( topq, j ), ldq )
660  ppwo = nblst*nblst + 1
661  j0 = j - nnb
662  DO j = j0, jcol+1, -nnb
663  IF ( initq ) THEN
664  topq = max( 2, j - jcol + 1 )
665  nh = ihi - topq + 1
666  END IF
667  IF ( blk22 ) THEN
668 *
669 * Exploit the structure of U.
670 *
671  CALL zunm22( 'Right', 'No Transpose', nh, 2*nnb,
672  $ nnb, nnb, work( ppwo ), 2*nnb,
673  $ q( topq, j ), ldq, work( pw ),
674  $ lwork-pw+1, ierr )
675  ELSE
676 *
677 * Ignore the structure of U.
678 *
679  CALL zgemm( 'No Transpose', 'No Transpose', nh,
680  $ 2*nnb, 2*nnb, cone, q( topq, j ), ldq,
681  $ work( ppwo ), 2*nnb, czero, work( pw ),
682  $ nh )
683  CALL zlacpy( 'All', nh, 2*nnb, work( pw ), nh,
684  $ q( topq, j ), ldq )
685  END IF
686  ppwo = ppwo + 4*nnb*nnb
687  END DO
688  END IF
689 *
690 * Accumulate right Givens rotations if required.
691 *
692  IF ( wantz .OR. top.GT.0 ) THEN
693 *
694 * Initialize small unitary factors that will hold the
695 * accumulated Givens rotations in workspace.
696 *
697  CALL zlaset( 'All', nblst, nblst, czero, cone, work,
698  $ nblst )
699  pw = nblst * nblst + 1
700  DO i = 1, n2nb
701  CALL zlaset( 'All', 2*nnb, 2*nnb, czero, cone,
702  $ work( pw ), 2*nnb )
703  pw = pw + 4*nnb*nnb
704  END DO
705 *
706 * Accumulate Givens rotations into workspace array.
707 *
708  DO j = jcol, jcol+nnb-1
709  ppw = ( nblst + 1 )*( nblst - 2 ) - j + jcol + 1
710  len = 2 + j - jcol
711  jrow = j + n2nb*nnb + 2
712  DO i = ihi, jrow, -1
713  ctemp = a( i, j )
714  a( i, j ) = czero
715  s = b( i, j )
716  b( i, j ) = czero
717  DO jj = ppw, ppw+len-1
718  temp = work( jj + nblst )
719  work( jj + nblst ) = ctemp*temp -
720  $ dconjg( s )*work( jj )
721  work( jj ) = s*temp + ctemp*work( jj )
722  END DO
723  len = len + 1
724  ppw = ppw - nblst - 1
725  END DO
726 *
727  ppwo = nblst*nblst + ( nnb+j-jcol-1 )*2*nnb + nnb
728  j0 = jrow - nnb
729  DO jrow = j0, j+2, -nnb
730  ppw = ppwo
731  len = 2 + j - jcol
732  DO i = jrow+nnb-1, jrow, -1
733  ctemp = a( i, j )
734  a( i, j ) = czero
735  s = b( i, j )
736  b( i, j ) = czero
737  DO jj = ppw, ppw+len-1
738  temp = work( jj + 2*nnb )
739  work( jj + 2*nnb ) = ctemp*temp -
740  $ dconjg( s )*work( jj )
741  work( jj ) = s*temp + ctemp*work( jj )
742  END DO
743  len = len + 1
744  ppw = ppw - 2*nnb - 1
745  END DO
746  ppwo = ppwo + 4*nnb*nnb
747  END DO
748  END DO
749  ELSE
750 *
751  CALL zlaset( 'Lower', ihi - jcol - 1, nnb, czero, czero,
752  $ a( jcol + 2, jcol ), lda )
753  CALL zlaset( 'Lower', ihi - jcol - 1, nnb, czero, czero,
754  $ b( jcol + 2, jcol ), ldb )
755  END IF
756 *
757 * Apply accumulated unitary matrices to A and B.
758 *
759  IF ( top.GT.0 ) THEN
760  j = ihi - nblst + 1
761  CALL zgemm( 'No Transpose', 'No Transpose', top,
762  $ nblst, nblst, cone, a( 1, j ), lda,
763  $ work, nblst, czero, work( pw ), top )
764  CALL zlacpy( 'All', top, nblst, work( pw ), top,
765  $ a( 1, j ), lda )
766  ppwo = nblst*nblst + 1
767  j0 = j - nnb
768  DO j = j0, jcol+1, -nnb
769  IF ( blk22 ) THEN
770 *
771 * Exploit the structure of U.
772 *
773  CALL zunm22( 'Right', 'No Transpose', top, 2*nnb,
774  $ nnb, nnb, work( ppwo ), 2*nnb,
775  $ a( 1, j ), lda, work( pw ),
776  $ lwork-pw+1, ierr )
777  ELSE
778 *
779 * Ignore the structure of U.
780 *
781  CALL zgemm( 'No Transpose', 'No Transpose', top,
782  $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
783  $ work( ppwo ), 2*nnb, czero,
784  $ work( pw ), top )
785  CALL zlacpy( 'All', top, 2*nnb, work( pw ), top,
786  $ a( 1, j ), lda )
787  END IF
788  ppwo = ppwo + 4*nnb*nnb
789  END DO
790 *
791  j = ihi - nblst + 1
792  CALL zgemm( 'No Transpose', 'No Transpose', top,
793  $ nblst, nblst, cone, b( 1, j ), ldb,
794  $ work, nblst, czero, work( pw ), top )
795  CALL zlacpy( 'All', top, nblst, work( pw ), top,
796  $ b( 1, j ), ldb )
797  ppwo = nblst*nblst + 1
798  j0 = j - nnb
799  DO j = j0, jcol+1, -nnb
800  IF ( blk22 ) THEN
801 *
802 * Exploit the structure of U.
803 *
804  CALL zunm22( 'Right', 'No Transpose', top, 2*nnb,
805  $ nnb, nnb, work( ppwo ), 2*nnb,
806  $ b( 1, j ), ldb, work( pw ),
807  $ lwork-pw+1, ierr )
808  ELSE
809 *
810 * Ignore the structure of U.
811 *
812  CALL zgemm( 'No Transpose', 'No Transpose', top,
813  $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
814  $ work( ppwo ), 2*nnb, czero,
815  $ work( pw ), top )
816  CALL zlacpy( 'All', top, 2*nnb, work( pw ), top,
817  $ b( 1, j ), ldb )
818  END IF
819  ppwo = ppwo + 4*nnb*nnb
820  END DO
821  END IF
822 *
823 * Apply accumulated unitary matrices to Z.
824 *
825  IF( wantz ) THEN
826  j = ihi - nblst + 1
827  IF ( initq ) THEN
828  topq = max( 2, j - jcol + 1 )
829  nh = ihi - topq + 1
830  ELSE
831  topq = 1
832  nh = n
833  END IF
834  CALL zgemm( 'No Transpose', 'No Transpose', nh,
835  $ nblst, nblst, cone, z( topq, j ), ldz,
836  $ work, nblst, czero, work( pw ), nh )
837  CALL zlacpy( 'All', nh, nblst, work( pw ), nh,
838  $ z( topq, j ), ldz )
839  ppwo = nblst*nblst + 1
840  j0 = j - nnb
841  DO j = j0, jcol+1, -nnb
842  IF ( initq ) THEN
843  topq = max( 2, j - jcol + 1 )
844  nh = ihi - topq + 1
845  END IF
846  IF ( blk22 ) THEN
847 *
848 * Exploit the structure of U.
849 *
850  CALL zunm22( 'Right', 'No Transpose', nh, 2*nnb,
851  $ nnb, nnb, work( ppwo ), 2*nnb,
852  $ z( topq, j ), ldz, work( pw ),
853  $ lwork-pw+1, ierr )
854  ELSE
855 *
856 * Ignore the structure of U.
857 *
858  CALL zgemm( 'No Transpose', 'No Transpose', nh,
859  $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
860  $ work( ppwo ), 2*nnb, czero, work( pw ),
861  $ nh )
862  CALL zlacpy( 'All', nh, 2*nnb, work( pw ), nh,
863  $ z( topq, j ), ldz )
864  END IF
865  ppwo = ppwo + 4*nnb*nnb
866  END DO
867  END IF
868  END DO
869  END IF
870 *
871 * Use unblocked code to reduce the rest of the matrix
872 * Avoid re-initialization of modified Q and Z.
873 *
874  compq2 = compq
875  compz2 = compz
876  IF ( jcol.NE.ilo ) THEN
877  IF ( wantq )
878  $ compq2 = 'V'
879  IF ( wantz )
880  $ compz2 = 'V'
881  END IF
882 *
883  IF ( jcol.LT.ihi )
884  $ CALL zgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
885  $ ldq, z, ldz, ierr )
886  work( 1 ) = dcmplx( lwkopt )
887 *
888  RETURN
889 *
890 * End of ZGGHD3
891 *
892  END
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:147
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: zrot.f:103
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
ZGGHD3
Definition: zgghd3.f:227
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:204
subroutine zunm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
ZUNM22 multiplies a general matrix by a banded unitary matrix.
Definition: zunm22.f:162