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