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