SRC/cgghd3.f

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