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.GE. IF( LWORK6*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.LT..OR..GE. IF( NBNBMIN NBNH ) 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.EQ. BLK22 = KACC222
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.OR..GT. IF ( WANTZ TOP0 ) 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.GT. IF ( TOP0 ) 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