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