OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cgghd3.f
Go to the documentation of this file.
1*> \brief \b CGGHD3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGGHD3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgghd3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgghd3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgghd3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
22* $ LDQ, Z, LDZ, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER COMPQ, COMPZ
26* INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
27* ..
28* .. Array Arguments ..
29* COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30* $ Z( LDZ, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*>
40*> CGGHD3 reduces a pair of complex matrices (A,B) to generalized upper
41*> Hessenberg form using unitary transformations, where A is a
42*> general matrix and B is upper triangular. The form of the
43*> generalized eigenvalue problem is
44*> A*x = lambda*B*x,
45*> and B is typically made upper triangular by computing its QR
46*> factorization and moving the unitary matrix Q to the left side
47*> of the equation.
48*>
49*> This subroutine simultaneously reduces A to a Hessenberg matrix H:
50*> Q**H*A*Z = H
51*> and transforms B to another upper triangular matrix T:
52*> Q**H*B*Z = T
53*> in order to reduce the problem to its standard form
54*> H*y = lambda*T*y
55*> where y = Z**H*x.
56*>
57*> The unitary matrices Q and Z are determined as products of Givens
58*> rotations. They may either be formed explicitly, or they may be
59*> postmultiplied into input matrices Q1 and Z1, so that
60*>
61*> Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
62*>
63*> Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
64*>
65*> If Q1 is the unitary matrix from the QR factorization of B in the
66*> original equation A*x = lambda*B*x, then CGGHD3 reduces the original
67*> problem to generalized Hessenberg form.
68*>
69*> This is a blocked variant of CGGHRD, using matrix-matrix
70*> multiplications for parts of the computation to enhance performance.
71*> \endverbatim
72*
73* Arguments:
74* ==========
75*
76*> \param[in] COMPQ
77*> \verbatim
78*> COMPQ is CHARACTER*1
79*> = 'N': do not compute Q;
80*> = 'I': Q is initialized to the unit matrix, and the
81*> unitary matrix Q is returned;
82*> = 'V': Q must contain a unitary matrix Q1 on entry,
83*> and the product Q1*Q is returned.
84*> \endverbatim
85*>
86*> \param[in] COMPZ
87*> \verbatim
88*> COMPZ is CHARACTER*1
89*> = 'N': do not compute Z;
90*> = 'I': Z is initialized to the unit matrix, and the
91*> unitary matrix Z is returned;
92*> = 'V': Z must contain a unitary matrix Z1 on entry,
93*> and the product Z1*Z is returned.
94*> \endverbatim
95*>
96*> \param[in] N
97*> \verbatim
98*> N is INTEGER
99*> The order of the matrices A and B. N >= 0.
100*> \endverbatim
101*>
102*> \param[in] ILO
103*> \verbatim
104*> ILO is INTEGER
105*> \endverbatim
106*>
107*> \param[in] IHI
108*> \verbatim
109*> IHI is INTEGER
110*>
111*> ILO and IHI mark the rows and columns of A which are to be
112*> reduced. It is assumed that A is already upper triangular
113*> in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
114*> normally set by a previous call to CGGBAL; otherwise they
115*> should be set to 1 and N respectively.
116*> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
117*> \endverbatim
118*>
119*> \param[in,out] A
120*> \verbatim
121*> A is COMPLEX array, dimension (LDA, N)
122*> On entry, the N-by-N general matrix to be reduced.
123*> On exit, the upper triangle and the first subdiagonal of A
124*> are overwritten with the upper Hessenberg matrix H, and the
125*> rest is set to zero.
126*> \endverbatim
127*>
128*> \param[in] LDA
129*> \verbatim
130*> LDA is INTEGER
131*> The leading dimension of the array A. LDA >= max(1,N).
132*> \endverbatim
133*>
134*> \param[in,out] B
135*> \verbatim
136*> B is COMPLEX array, dimension (LDB, N)
137*> On entry, the N-by-N upper triangular matrix B.
138*> On exit, the upper triangular matrix T = Q**H B Z. The
139*> elements below the diagonal are set to zero.
140*> \endverbatim
141*>
142*> \param[in] LDB
143*> \verbatim
144*> LDB is INTEGER
145*> The leading dimension of the array B. LDB >= max(1,N).
146*> \endverbatim
147*>
148*> \param[in,out] Q
149*> \verbatim
150*> Q is COMPLEX array, dimension (LDQ, N)
151*> On entry, if COMPQ = 'V', the unitary matrix Q1, typically
152*> from the QR factorization of B.
153*> On exit, if COMPQ='I', the unitary matrix Q, and if
154*> COMPQ = 'V', the product Q1*Q.
155*> Not referenced if COMPQ='N'.
156*> \endverbatim
157*>
158*> \param[in] LDQ
159*> \verbatim
160*> LDQ is INTEGER
161*> The leading dimension of the array Q.
162*> LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
163*> \endverbatim
164*>
165*> \param[in,out] Z
166*> \verbatim
167*> Z is COMPLEX array, dimension (LDZ, N)
168*> On entry, if COMPZ = 'V', the unitary matrix Z1.
169*> On exit, if COMPZ='I', the unitary matrix Z, and if
170*> COMPZ = 'V', the product Z1*Z.
171*> Not referenced if COMPZ='N'.
172*> \endverbatim
173*>
174*> \param[in] LDZ
175*> \verbatim
176*> LDZ is INTEGER
177*> The leading dimension of the array Z.
178*> LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
179*> \endverbatim
180*>
181*> \param[out] WORK
182*> \verbatim
183*> WORK is COMPLEX array, dimension (LWORK)
184*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
185*> \endverbatim
186*>
187*> \param[in] LWORK
188*> \verbatim
189*> LWORK is INTEGER
190*> The length of the array WORK. LWORK >= 1.
191*> For optimum performance LWORK >= 6*N*NB, where NB is the
192*> optimal blocksize.
193*>
194*> If LWORK = -1, then a workspace query is assumed; the routine
195*> only calculates the optimal size of the WORK array, returns
196*> this value as the first entry of the WORK array, and no error
197*> message related to LWORK is issued by XERBLA.
198*> \endverbatim
199*>
200*> \param[out] INFO
201*> \verbatim
202*> INFO is INTEGER
203*> = 0: successful exit.
204*> < 0: if INFO = -i, the i-th argument had an illegal value.
205*> \endverbatim
206*
207* Authors:
208* ========
209*
210*> \author Univ. of Tennessee
211*> \author Univ. of California Berkeley
212*> \author Univ. of Colorado Denver
213*> \author NAG Ltd.
214*
215*> \ingroup complexOTHERcomputational
216*
217*> \par Further Details:
218* =====================
219*>
220*> \verbatim
221*>
222*> This routine reduces A to Hessenberg form and maintains B in triangular form
223*> using a blocked variant of Moler and Stewart's original algorithm,
224*> as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
225*> (BIT 2008).
226*> \endverbatim
227*>
228* =====================================================================
229 SUBROUTINE cgghd3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
230 $ LDQ, Z, LDZ, WORK, LWORK, INFO )
231*
232* -- LAPACK computational routine --
233* -- LAPACK is a software package provided by Univ. of Tennessee, --
234* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235*
236*
237 IMPLICIT NONE
238*
239* .. Scalar Arguments ..
240 CHARACTER COMPQ, COMPZ
241 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
242* ..
243* .. Array Arguments ..
244 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
245 $ z( ldz, * ), work( * )
246* ..
247*
248* =====================================================================
249*
250* .. Parameters ..
251 COMPLEX CONE, CZERO
252 parameter( cone = ( 1.0e+0, 0.0e+0 ),
253 $ czero = ( 0.0e+0, 0.0e+0 ) )
254* ..
255* .. Local Scalars ..
256 LOGICAL BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
257 CHARACTER*1 COMPQ2, COMPZ2
258 INTEGER COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
259 $ kacc22, len, lwkopt, n2nb, nb, nblst, nbmin,
260 $ nh, nnb, nx, ppw, ppwo, pw, top, topq
261 REAL C
262 COMPLEX C1, C2, CTEMP, S, S1, S2, TEMP, TEMP1, TEMP2,
263 $ temp3
264* ..
265* .. External Functions ..
266 LOGICAL LSAME
267 INTEGER ILAENV
268 EXTERNAL ilaenv, lsame
269* ..
270* .. External Subroutines ..
271 EXTERNAL cgghrd, clartg, claset, cunm22, crot, cgemm,
272 $ cgemv, ctrmv, clacpy, xerbla
273* ..
274* .. Intrinsic Functions ..
275 INTRINSIC real, cmplx, conjg, max
276* ..
277* .. Executable Statements ..
278*
279* Decode and test the input parameters.
280*
281 info = 0
282 nb = ilaenv( 1, 'cgghd3', ' ', N, ILO, IHI, -1 )
283 LWKOPT = MAX( 6*N*NB, 1 )
284 WORK( 1 ) = CMPLX( LWKOPT )
285 INITQ = LSAME( COMPQ, 'i' )
286.OR. WANTQ = INITQ LSAME( COMPQ, 'v' )
287 INITZ = LSAME( COMPZ, 'i' )
288.OR. WANTZ = INITZ LSAME( COMPZ, 'v' )
289.EQ. LQUERY = ( LWORK-1 )
290*
291.NOT. IF( LSAME( COMPQ, 'n.AND..NOT.' ) WANTQ ) THEN
292 INFO = -1
293.NOT. ELSE IF( LSAME( COMPZ, 'n.AND..NOT.' ) WANTZ ) THEN
294 INFO = -2
295.LT. ELSE IF( N0 ) THEN
296 INFO = -3
297.LT. ELSE IF( ILO1 ) THEN
298 INFO = -4
299.GT..OR..LT. ELSE IF( IHIN IHIILO-1 ) THEN
300 INFO = -5
301.LT. ELSE IF( LDAMAX( 1, N ) ) THEN
302 INFO = -7
303.LT. ELSE IF( LDBMAX( 1, N ) ) THEN
304 INFO = -9
305.AND..LT..OR..LT. ELSE IF( ( WANTQ LDQN ) LDQ1 ) THEN
306 INFO = -11
307.AND..LT..OR..LT. ELSE IF( ( WANTZ LDZN ) LDZ1 ) THEN
308 INFO = -13
309.LT..AND..NOT. ELSE IF( LWORK1 LQUERY ) THEN
310 INFO = -15
311 END IF
312.NE. IF( INFO0 ) THEN
313 CALL XERBLA( 'cgghd3', -INFO )
314 RETURN
315 ELSE IF( LQUERY ) THEN
316 RETURN
317 END IF
318*
319* Initialize Q and Z if desired.
320*
321 IF( INITQ )
322 $ CALL CLASET( 'all', N, N, CZERO, CONE, Q, LDQ )
323 IF( INITZ )
324 $ CALL CLASET( 'all', N, N, CZERO, CONE, Z, LDZ )
325*
326* Zero out lower triangle of B.
327*
328.GT. IF( N1 )
329 $ CALL CLASET( 'lower', N-1, N-1, CZERO, CZERO, B(2, 1), LDB )
330*
331* Quick return if possible
332*
333 NH = IHI - ILO + 1
334.LE. IF( NH1 ) THEN
335 WORK( 1 ) = CONE
336 RETURN
337 END IF
338*
339* Determine the blocksize.
340*
341 NBMIN = ILAENV( 2, 'cgghd3', ' ', N, ILO, IHI, -1 )
342.GT..AND..LT. IF( NB1 NBNH ) THEN
343*
344* Determine when to use unblocked instead of blocked code.
345*
346 NX = MAX( NB, ILAENV( 3, 'cgghd3', ' ', N, ILO, IHI, -1 ) )
347.LT. IF( NXNH ) THEN
348*
349* Determine if workspace is large enough for blocked code.
350*
351.LT. IF( LWORKLWKOPT ) THEN
352*
353* Not enough workspace to use optimal NB: determine the
354* minimum value of NB, and reduce NB or force use of
355* unblocked code.
356*
357 NBMIN = MAX( 2, ILAENV( 2, 'cgghd3', ' ', N, ILO, IHI,
358 $ -1 ) )
359.GE. IF( LWORK6*N*NBMIN ) THEN
360 NB = LWORK / ( 6*N )
361 ELSE
362 NB = 1
363 END IF
364 END IF
365 END IF
366 END IF
367*
368.LT..OR..GE. IF( NBNBMIN NBNH ) THEN
369*
370* Use unblocked code below
371*
372 JCOL = ILO
373*
374 ELSE
375*
376* Use blocked code
377*
378 KACC22 = ILAENV( 16, 'cgghd3', ' ', N, ILO, IHI, -1 )
379.EQ. BLK22 = KACC222
380 DO JCOL = ILO, IHI-2, NB
381 NNB = MIN( NB, IHI-JCOL-1 )
382*
383* Initialize small unitary factors that will hold the
384* accumulated Givens rotations in workspace.
385* N2NB denotes the number of 2*NNB-by-2*NNB factors
386* NBLST denotes the (possibly smaller) order of the last
387* factor.
388*
389 N2NB = ( IHI-JCOL-1 ) / NNB - 1
390 NBLST = IHI - JCOL - N2NB*NNB
391 CALL CLASET( 'all', NBLST, NBLST, CZERO, CONE, WORK, NBLST )
392 PW = NBLST * NBLST + 1
393 DO I = 1, N2NB
394 CALL CLASET( 'all', 2*NNB, 2*NNB, CZERO, CONE,
395 $ WORK( PW ), 2*NNB )
396 PW = PW + 4*NNB*NNB
397 END DO
398*
399* Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
400*
401 DO J = JCOL, JCOL+NNB-1
402*
403* Reduce Jth column of A. Store cosines and sines in Jth
404* column of A and B, respectively.
405*
406 DO I = IHI, J+2, -1
407 TEMP = A( I-1, J )
408 CALL CLARTG( TEMP, A( I, J ), C, S, A( I-1, J ) )
409 A( I, J ) = CMPLX( C )
410 B( I, J ) = S
411 END DO
412*
413* Accumulate Givens rotations into workspace array.
414*
415 PPW = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
416 LEN = 2 + J - JCOL
417 JROW = J + N2NB*NNB + 2
418 DO I = IHI, JROW, -1
419 CTEMP = A( I, J )
420 S = B( I, J )
421 DO JJ = PPW, PPW+LEN-1
422 TEMP = WORK( JJ + NBLST )
423 WORK( JJ + NBLST ) = CTEMP*TEMP - S*WORK( JJ )
424 WORK( JJ ) = CONJG( S )*TEMP + CTEMP*WORK( JJ )
425 END DO
426 LEN = LEN + 1
427 PPW = PPW - NBLST - 1
428 END DO
429*
430 PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
431 J0 = JROW - NNB
432 DO JROW = J0, J+2, -NNB
433 PPW = PPWO
434 LEN = 2 + J - JCOL
435 DO I = JROW+NNB-1, JROW, -1
436 CTEMP = A( I, J )
437 S = B( I, J )
438 DO JJ = PPW, PPW+LEN-1
439 TEMP = WORK( JJ + 2*NNB )
440 WORK( JJ + 2*NNB ) = CTEMP*TEMP - S*WORK( JJ )
441 WORK( JJ ) = CONJG( S )*TEMP + CTEMP*WORK( JJ )
442 END DO
443 LEN = LEN + 1
444 PPW = PPW - 2*NNB - 1
445 END DO
446 PPWO = PPWO + 4*NNB*NNB
447 END DO
448*
449* TOP denotes the number of top rows in A and B that will
450* not be updated during the next steps.
451*
452.LE. IF( JCOL2 ) THEN
453 TOP = 0
454 ELSE
455 TOP = JCOL
456 END IF
457*
458* Propagate transformations through B and replace stored
459* left sines/cosines by right sines/cosines.
460*
461 DO JJ = N, J+1, -1
462*
463* Update JJth column of B.
464*
465 DO I = MIN( JJ+1, IHI ), J+2, -1
466 CTEMP = A( I, J )
467 S = B( I, J )
468 TEMP = B( I, JJ )
469 B( I, JJ ) = CTEMP*TEMP - CONJG( S )*B( I-1, JJ )
470 B( I-1, JJ ) = S*TEMP + CTEMP*B( I-1, JJ )
471 END DO
472*
473* Annihilate B( JJ+1, JJ ).
474*
475.LT. IF( JJIHI ) THEN
476 TEMP = B( JJ+1, JJ+1 )
477 CALL CLARTG( TEMP, B( JJ+1, JJ ), C, S,
478 $ B( JJ+1, JJ+1 ) )
479 B( JJ+1, JJ ) = CZERO
480 CALL CROT( JJ-TOP, B( TOP+1, JJ+1 ), 1,
481 $ B( TOP+1, JJ ), 1, C, S )
482 A( JJ+1, J ) = CMPLX( C )
483 B( JJ+1, J ) = -CONJG( S )
484 END IF
485 END DO
486*
487* Update A by transformations from right.
488*
489 JJ = MOD( IHI-J-1, 3 )
490 DO I = IHI-J-3, JJ+1, -3
491 CTEMP = A( J+1+I, J )
492 S = -B( J+1+I, J )
493 C1 = A( J+2+I, J )
494 S1 = -B( J+2+I, J )
495 C2 = A( J+3+I, J )
496 S2 = -B( J+3+I, J )
497*
498 DO K = TOP+1, IHI
499 TEMP = A( K, J+I )
500 TEMP1 = A( K, J+I+1 )
501 TEMP2 = A( K, J+I+2 )
502 TEMP3 = A( K, J+I+3 )
503 A( K, J+I+3 ) = C2*TEMP3 + CONJG( S2 )*TEMP2
504 TEMP2 = -S2*TEMP3 + C2*TEMP2
505 A( K, J+I+2 ) = C1*TEMP2 + CONJG( S1 )*TEMP1
506 TEMP1 = -S1*TEMP2 + C1*TEMP1
507 A( K, J+I+1 ) = CTEMP*TEMP1 + CONJG( S )*TEMP
508 A( K, J+I ) = -S*TEMP1 + CTEMP*TEMP
509 END DO
510 END DO
511*
512.GT. IF( JJ0 ) THEN
513 DO I = JJ, 1, -1
514 C = DBLE( A( J+1+I, J ) )
515 CALL CROT( IHI-TOP, A( TOP+1, J+I+1 ), 1,
516 $ A( TOP+1, J+I ), 1, C,
517 $ -CONJG( B( J+1+I, J ) ) )
518 END DO
519 END IF
520*
521* Update (J+1)th column of A by transformations from left.
522*
523.LT. IF ( J JCOL + NNB - 1 ) THEN
524 LEN = 1 + J - JCOL
525*
526* Multiply with the trailing accumulated unitary
527* matrix, which takes the form
528*
529* [ U11 U12 ]
530* U = [ ],
531* [ U21 U22 ]
532*
533* where U21 is a LEN-by-LEN matrix and U12 is lower
534* triangular.
535*
536 JROW = IHI - NBLST + 1
537 CALL CGEMV( 'conjugate', NBLST, LEN, CONE, WORK,
538 $ NBLST, A( JROW, J+1 ), 1, CZERO,
539 $ WORK( PW ), 1 )
540 PPW = PW + LEN
541 DO I = JROW, JROW+NBLST-LEN-1
542 WORK( PPW ) = A( I, J+1 )
543 PPW = PPW + 1
544 END DO
545 CALL CTRMV( 'lower', 'conjugate', 'non-unit',
546 $ NBLST-LEN, WORK( LEN*NBLST + 1 ), NBLST,
547 $ WORK( PW+LEN ), 1 )
548 CALL CGEMV( 'conjugate', LEN, NBLST-LEN, CONE,
549 $ WORK( (LEN+1)*NBLST - LEN + 1 ), NBLST,
550 $ A( JROW+NBLST-LEN, J+1 ), 1, CONE,
551 $ WORK( PW+LEN ), 1 )
552 PPW = PW
553 DO I = JROW, JROW+NBLST-1
554 A( I, J+1 ) = WORK( PPW )
555 PPW = PPW + 1
556 END DO
557*
558* Multiply with the other accumulated unitary
559* matrices, which take the form
560*
561* [ U11 U12 0 ]
562* [ ]
563* U = [ U21 U22 0 ],
564* [ ]
565* [ 0 0 I ]
566*
567* where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
568* matrix, U21 is a LEN-by-LEN upper triangular matrix
569* and U12 is an NNB-by-NNB lower triangular matrix.
570*
571 PPWO = 1 + NBLST*NBLST
572 J0 = JROW - NNB
573 DO JROW = J0, JCOL+1, -NNB
574 PPW = PW + LEN
575 DO I = JROW, JROW+NNB-1
576 WORK( PPW ) = A( I, J+1 )
577 PPW = PPW + 1
578 END DO
579 PPW = PW
580 DO I = JROW+NNB, JROW+NNB+LEN-1
581 WORK( PPW ) = A( I, J+1 )
582 PPW = PPW + 1
583 END DO
584 CALL CTRMV( 'upper', 'conjugate', 'non-unit', LEN,
585 $ WORK( PPWO + NNB ), 2*NNB, WORK( PW ),
586 $ 1 )
587 CALL CTRMV( 'lower', 'conjugate', 'non-unit', NNB,
588 $ WORK( PPWO + 2*LEN*NNB ),
589 $ 2*NNB, WORK( PW + LEN ), 1 )
590 CALL CGEMV( 'conjugate', NNB, LEN, CONE,
591 $ WORK( PPWO ), 2*NNB, A( JROW, J+1 ), 1,
592 $ CONE, WORK( PW ), 1 )
593 CALL CGEMV( 'conjugate', LEN, NNB, CONE,
594 $ WORK( PPWO + 2*LEN*NNB + NNB ), 2*NNB,
595 $ A( JROW+NNB, J+1 ), 1, CONE,
596 $ WORK( PW+LEN ), 1 )
597 PPW = PW
598 DO I = JROW, JROW+LEN+NNB-1
599 A( I, J+1 ) = WORK( PPW )
600 PPW = PPW + 1
601 END DO
602 PPWO = PPWO + 4*NNB*NNB
603 END DO
604 END IF
605 END DO
606*
607* Apply accumulated unitary matrices to A.
608*
609 COLA = N - JCOL - NNB + 1
610 J = IHI - NBLST + 1
611 CALL CGEMM( 'conjugate', 'no transpose', NBLST,
612 $ COLA, NBLST, CONE, WORK, NBLST,
613 $ A( J, JCOL+NNB ), LDA, CZERO, WORK( PW ),
614 $ NBLST )
615 CALL CLACPY( 'all', NBLST, COLA, WORK( PW ), NBLST,
616 $ A( J, JCOL+NNB ), LDA )
617 PPWO = NBLST*NBLST + 1
618 J0 = J - NNB
619 DO J = J0, JCOL+1, -NNB
620 IF ( BLK22 ) THEN
621*
622* Exploit the structure of
623*
624* [ U11 U12 ]
625* U = [ ]
626* [ U21 U22 ],
627*
628* where all blocks are NNB-by-NNB, U21 is upper
629* triangular and U12 is lower triangular.
630*
631 CALL CUNM22( 'left', 'conjugate', 2*NNB, COLA, NNB,
632 $ NNB, WORK( PPWO ), 2*NNB,
633 $ A( J, JCOL+NNB ), LDA, WORK( PW ),
634 $ LWORK-PW+1, IERR )
635 ELSE
636*
637* Ignore the structure of U.
638*
639 CALL CGEMM( 'conjugate', 'no transpose', 2*NNB,
640 $ COLA, 2*NNB, CONE, WORK( PPWO ), 2*NNB,
641 $ A( J, JCOL+NNB ), LDA, CZERO, WORK( PW ),
642 $ 2*NNB )
643 CALL CLACPY( 'all', 2*NNB, COLA, WORK( PW ), 2*NNB,
644 $ A( J, JCOL+NNB ), LDA )
645 END IF
646 PPWO = PPWO + 4*NNB*NNB
647 END DO
648*
649* Apply accumulated unitary matrices to Q.
650*
651 IF( WANTQ ) THEN
652 J = IHI - NBLST + 1
653 IF ( INITQ ) THEN
654 TOPQ = MAX( 2, J - JCOL + 1 )
655 NH = IHI - TOPQ + 1
656 ELSE
657 TOPQ = 1
658 NH = N
659 END IF
660 CALL CGEMM( 'no transpose', 'no transpose', NH,
661 $ NBLST, NBLST, CONE, Q( TOPQ, J ), LDQ,
662 $ WORK, NBLST, CZERO, WORK( PW ), NH )
663 CALL CLACPY( 'all', NH, NBLST, WORK( PW ), NH,
664 $ Q( TOPQ, J ), LDQ )
665 PPWO = NBLST*NBLST + 1
666 J0 = J - NNB
667 DO J = J0, JCOL+1, -NNB
668 IF ( INITQ ) THEN
669 TOPQ = MAX( 2, J - JCOL + 1 )
670 NH = IHI - TOPQ + 1
671 END IF
672 IF ( BLK22 ) THEN
673*
674* Exploit the structure of U.
675*
676 CALL CUNM22( 'right', 'no transpose', NH, 2*NNB,
677 $ NNB, NNB, WORK( PPWO ), 2*NNB,
678 $ Q( TOPQ, J ), LDQ, WORK( PW ),
679 $ LWORK-PW+1, IERR )
680 ELSE
681*
682* Ignore the structure of U.
683*
684 CALL CGEMM( 'no transpose', 'no transpose', NH,
685 $ 2*NNB, 2*NNB, CONE, Q( TOPQ, J ), LDQ,
686 $ WORK( PPWO ), 2*NNB, CZERO, WORK( PW ),
687 $ NH )
688 CALL CLACPY( 'all', NH, 2*NNB, WORK( PW ), NH,
689 $ Q( TOPQ, J ), LDQ )
690 END IF
691 PPWO = PPWO + 4*NNB*NNB
692 END DO
693 END IF
694*
695* Accumulate right Givens rotations if required.
696*
697.OR..GT. IF ( WANTZ TOP0 ) THEN
698*
699* Initialize small unitary factors that will hold the
700* accumulated Givens rotations in workspace.
701*
702 CALL CLASET( 'all', NBLST, NBLST, CZERO, CONE, WORK,
703 $ NBLST )
704 PW = NBLST * NBLST + 1
705 DO I = 1, N2NB
706 CALL CLASET( 'all', 2*NNB, 2*NNB, CZERO, CONE,
707 $ WORK( PW ), 2*NNB )
708 PW = PW + 4*NNB*NNB
709 END DO
710*
711* Accumulate Givens rotations into workspace array.
712*
713 DO J = JCOL, JCOL+NNB-1
714 PPW = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
715 LEN = 2 + J - JCOL
716 JROW = J + N2NB*NNB + 2
717 DO I = IHI, JROW, -1
718 CTEMP = A( I, J )
719 A( I, J ) = CZERO
720 S = B( I, J )
721 B( I, J ) = CZERO
722 DO JJ = PPW, PPW+LEN-1
723 TEMP = WORK( JJ + NBLST )
724 WORK( JJ + NBLST ) = CTEMP*TEMP -
725 $ CONJG( S )*WORK( JJ )
726 WORK( JJ ) = S*TEMP + CTEMP*WORK( JJ )
727 END DO
728 LEN = LEN + 1
729 PPW = PPW - NBLST - 1
730 END DO
731*
732 PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
733 J0 = JROW - NNB
734 DO JROW = J0, J+2, -NNB
735 PPW = PPWO
736 LEN = 2 + J - JCOL
737 DO I = JROW+NNB-1, JROW, -1
738 CTEMP = A( I, J )
739 A( I, J ) = CZERO
740 S = B( I, J )
741 B( I, J ) = CZERO
742 DO JJ = PPW, PPW+LEN-1
743 TEMP = WORK( JJ + 2*NNB )
744 WORK( JJ + 2*NNB ) = CTEMP*TEMP -
745 $ CONJG( S )*WORK( JJ )
746 WORK( JJ ) = S*TEMP + CTEMP*WORK( JJ )
747 END DO
748 LEN = LEN + 1
749 PPW = PPW - 2*NNB - 1
750 END DO
751 PPWO = PPWO + 4*NNB*NNB
752 END DO
753 END DO
754 ELSE
755*
756 CALL CLASET( 'lower', IHI - JCOL - 1, NNB, CZERO, CZERO,
757 $ A( JCOL + 2, JCOL ), LDA )
758 CALL CLASET( 'lower', IHI - JCOL - 1, NNB, CZERO, CZERO,
759 $ B( JCOL + 2, JCOL ), LDB )
760 END IF
761*
762* Apply accumulated unitary matrices to A and B.
763*
764.GT. IF ( TOP0 ) THEN
765 J = IHI - NBLST + 1
766 CALL CGEMM( 'no transpose', 'no transpose', TOP,
767 $ NBLST, NBLST, CONE, A( 1, J ), LDA,
768 $ WORK, NBLST, CZERO, WORK( PW ), TOP )
769 CALL CLACPY( 'all', TOP, NBLST, WORK( PW ), TOP,
770 $ A( 1, J ), LDA )
771 PPWO = NBLST*NBLST + 1
772 J0 = J - NNB
773 DO J = J0, JCOL+1, -NNB
774 IF ( BLK22 ) THEN
775*
776* Exploit the structure of U.
777*
778 CALL CUNM22( 'right', 'no transpose', TOP, 2*NNB,
779 $ NNB, NNB, WORK( PPWO ), 2*NNB,
780 $ A( 1, J ), LDA, WORK( PW ),
781 $ LWORK-PW+1, IERR )
782 ELSE
783*
784* Ignore the structure of U.
785*
786 CALL CGEMM( 'no transpose', 'No Transpose', top,
787 $ 2*nnb, 2*nnb, cone, a( 1, j ), lda,
788 $ work( ppwo ), 2*nnb, czero,
789 $ work( pw ), top )
790 CALL clacpy( 'All', top, 2*nnb, work( pw ), top,
791 $ a( 1, j ), lda )
792 END IF
793 ppwo = ppwo + 4*nnb*nnb
794 END DO
795*
796 j = ihi - nblst + 1
797 CALL cgemm( 'No Transpose', 'No Transpose', top,
798 $ nblst, nblst, cone, b( 1, j ), ldb,
799 $ work, nblst, czero, work( pw ), top )
800 CALL clacpy( 'All', top, nblst, work( pw ), top,
801 $ b( 1, j ), ldb )
802 ppwo = nblst*nblst + 1
803 j0 = j - nnb
804 DO j = j0, jcol+1, -nnb
805 IF ( blk22 ) THEN
806*
807* Exploit the structure of U.
808*
809 CALL cunm22( 'Right', 'No Transpose', top, 2*nnb,
810 $ nnb, nnb, work( ppwo ), 2*nnb,
811 $ b( 1, j ), ldb, work( pw ),
812 $ lwork-pw+1, ierr )
813 ELSE
814*
815* Ignore the structure of U.
816*
817 CALL cgemm( 'No Transpose', 'No Transpose', top,
818 $ 2*nnb, 2*nnb, cone, b( 1, j ), ldb,
819 $ work( ppwo ), 2*nnb, czero,
820 $ work( pw ), top )
821 CALL clacpy( 'All', top, 2*nnb, work( pw ), top,
822 $ b( 1, j ), ldb )
823 END IF
824 ppwo = ppwo + 4*nnb*nnb
825 END DO
826 END IF
827*
828* Apply accumulated unitary matrices to Z.
829*
830 IF( wantz ) THEN
831 j = ihi - nblst + 1
832 IF ( initq ) THEN
833 topq = max( 2, j - jcol + 1 )
834 nh = ihi - topq + 1
835 ELSE
836 topq = 1
837 nh = n
838 END IF
839 CALL cgemm( 'No Transpose', 'No Transpose', nh,
840 $ nblst, nblst, cone, z( topq, j ), ldz,
841 $ work, nblst, czero, work( pw ), nh )
842 CALL clacpy( 'All', nh, nblst, work( pw ), nh,
843 $ z( topq, j ), ldz )
844 ppwo = nblst*nblst + 1
845 j0 = j - nnb
846 DO j = j0, jcol+1, -nnb
847 IF ( initq ) THEN
848 topq = max( 2, j - jcol + 1 )
849 nh = ihi - topq + 1
850 END IF
851 IF ( blk22 ) THEN
852*
853* Exploit the structure of U.
854*
855 CALL cunm22( 'Right', 'No Transpose', nh, 2*nnb,
856 $ nnb, nnb, work( ppwo ), 2*nnb,
857 $ z( topq, j ), ldz, work( pw ),
858 $ lwork-pw+1, ierr )
859 ELSE
860*
861* Ignore the structure of U.
862*
863 CALL cgemm( 'No Transpose', 'No Transpose', nh,
864 $ 2*nnb, 2*nnb, cone, z( topq, j ), ldz,
865 $ work( ppwo ), 2*nnb, czero, work( pw ),
866 $ nh )
867 CALL clacpy( 'All', nh, 2*nnb, work( pw ), nh,
868 $ z( topq, j ), ldz )
869 END IF
870 ppwo = ppwo + 4*nnb*nnb
871 END DO
872 END IF
873 END DO
874 END IF
875*
876* Use unblocked code to reduce the rest of the matrix
877* Avoid re-initialization of modified Q and Z.
878*
879 compq2 = compq
880 compz2 = compz
881 IF ( jcol.NE.ilo ) THEN
882 IF ( wantq )
883 $ compq2 = 'V'
884 IF ( wantz )
885 $ compz2 = 'V'
886 END IF
887*
888 IF ( jcol.LT.ihi )
889 $ CALL cgghrd( compq2, compz2, n, jcol, ihi, a, lda, b, ldb, q,
890 $ ldq, z, ldz, ierr )
891 work( 1 ) = cmplx( lwkopt )
892*
893 RETURN
894*
895* End of CGGHD3
896*
897 END
cmplx
float cmplx[2]
Definition pblas.h:136
clartg
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:118
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
crot
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition crot.f:103
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
cunm22
subroutine cunm22(side, trans, m, n, n1, n2, q, ldq, c, ldc, work, lwork, info)
CUNM22 multiplies a general matrix by a banded unitary matrix.
Definition cunm22.f:162
cgghrd
subroutine cgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
CGGHRD
Definition cgghrd.f:204
cgghd3
subroutine cgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
CGGHD3
Definition cgghd3.f:231
ctrmv
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187
max
#define max(a, b)
Definition macros.h:21