OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:118
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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:103
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:103
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:106
zgghd3
subroutine zgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
ZGGHD3
Definition zgghd3.f:227
zgghrd
subroutine zgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
ZGGHRD
Definition zgghrd.f:204
ztrmv
subroutine ztrmv(uplo, trans, diag, n, a, lda, x, incx)
ZTRMV
Definition ztrmv.f:147
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:158
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
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:162
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21