SRC/zgghrd.f

   1 *> \brief \b ZGGHRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGGHRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgghrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgghrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgghrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
  22 *                          LDQ, Z, LDZ, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          COMPQ, COMPZ
  26 *       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  30 *      $                   Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> ZGGHRD 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 ZGGHRD reduces the original
  63 *> problem to generalized Hessenberg form.
  64 *> \endverbatim
  65 *
  66 *  Arguments:
  67 *  ==========
  68 *
  69 *> \param[in] COMPQ
  70 *> \verbatim
  71 *>          COMPQ is CHARACTER*1
  72 *>          = 'N': do not compute Q;
  73 *>          = 'I': Q is initialized to the unit matrix, and the
  74 *>                 unitary matrix Q is returned;
  75 *>          = 'V': Q must contain a unitary matrix Q1 on entry,
  76 *>                 and the product Q1*Q is returned.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] COMPZ
  80 *> \verbatim
  81 *>          COMPZ is CHARACTER*1
  82 *>          = 'N': do not compute Z;
  83 *>          = 'I': Z is initialized to the unit matrix, and the
  84 *>                 unitary matrix Z is returned;
  85 *>          = 'V': Z must contain a unitary matrix Z1 on entry,
  86 *>                 and the product Z1*Z is returned.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] N
  90 *> \verbatim
  91 *>          N is INTEGER
  92 *>          The order of the matrices A and B.  N >= 0.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] ILO
  96 *> \verbatim
  97 *>          ILO is INTEGER
  98 *> \endverbatim
  99 *>
 100 *> \param[in] IHI
 101 *> \verbatim
 102 *>          IHI is INTEGER
 103 *>
 104 *>          ILO and IHI mark the rows and columns of A which are to be
 105 *>          reduced.  It is assumed that A is already upper triangular
 106 *>          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
 107 *>          normally set by a previous call to ZGGBAL; otherwise they
 108 *>          should be set to 1 and N respectively.
 109 *>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
 110 *> \endverbatim
 111 *>
 112 *> \param[in,out] A
 113 *> \verbatim
 114 *>          A is COMPLEX*16 array, dimension (LDA, N)
 115 *>          On entry, the N-by-N general matrix to be reduced.
 116 *>          On exit, the upper triangle and the first subdiagonal of A
 117 *>          are overwritten with the upper Hessenberg matrix H, and the
 118 *>          rest is set to zero.
 119 *> \endverbatim
 120 *>
 121 *> \param[in] LDA
 122 *> \verbatim
 123 *>          LDA is INTEGER
 124 *>          The leading dimension of the array A.  LDA >= max(1,N).
 125 *> \endverbatim
 126 *>
 127 *> \param[in,out] B
 128 *> \verbatim
 129 *>          B is COMPLEX*16 array, dimension (LDB, N)
 130 *>          On entry, the N-by-N upper triangular matrix B.
 131 *>          On exit, the upper triangular matrix T = Q**H B Z.  The
 132 *>          elements below the diagonal are set to zero.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] LDB
 136 *> \verbatim
 137 *>          LDB is INTEGER
 138 *>          The leading dimension of the array B.  LDB >= max(1,N).
 139 *> \endverbatim
 140 *>
 141 *> \param[in,out] Q
 142 *> \verbatim
 143 *>          Q is COMPLEX*16 array, dimension (LDQ, N)
 144 *>          On entry, if COMPQ = 'V', the unitary matrix Q1, typically
 145 *>          from the QR factorization of B.
 146 *>          On exit, if COMPQ='I', the unitary matrix Q, and if
 147 *>          COMPQ = 'V', the product Q1*Q.
 148 *>          Not referenced if COMPQ='N'.
 149 *> \endverbatim
 150 *>
 151 *> \param[in] LDQ
 152 *> \verbatim
 153 *>          LDQ is INTEGER
 154 *>          The leading dimension of the array Q.
 155 *>          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
 156 *> \endverbatim
 157 *>
 158 *> \param[in,out] Z
 159 *> \verbatim
 160 *>          Z is COMPLEX*16 array, dimension (LDZ, N)
 161 *>          On entry, if COMPZ = 'V', the unitary matrix Z1.
 162 *>          On exit, if COMPZ='I', the unitary matrix Z, and if
 163 *>          COMPZ = 'V', the product Z1*Z.
 164 *>          Not referenced if COMPZ='N'.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] LDZ
 168 *> \verbatim
 169 *>          LDZ is INTEGER
 170 *>          The leading dimension of the array Z.
 171 *>          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
 172 *> \endverbatim
 173 *>
 174 *> \param[out] INFO
 175 *> \verbatim
 176 *>          INFO is INTEGER
 177 *>          = 0:  successful exit.
 178 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 179 *> \endverbatim
 180 *
 181 *  Authors:
 182 *  ========
 183 *
 184 *> \author Univ. of Tennessee
 185 *> \author Univ. of California Berkeley
 186 *> \author Univ. of Colorado Denver
 187 *> \author NAG Ltd.
 188 *
 189 *> \date November 2015
 190 *
 191 *> \ingroup complex16OTHERcomputational
 192 *
 193 *> \par Further Details:
 194 *  =====================
 195 *>
 196 *> \verbatim
 197 *>
 198 *>  This routine reduces A to Hessenberg and B to triangular form by
 199 *>  an unblocked reduction, as described in _Matrix_Computations_,
 200 *>  by Golub and van Loan (Johns Hopkins Press).
 201 *> \endverbatim
 202 *>
 203 *  =====================================================================
 204       SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
 205      $                   LDQ, Z, LDZ, INFO )
 206 *
 207 *  -- LAPACK computational routine (version 3.6.0) --
 208 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 209 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 210 *     November 2015
 211 *
 212 *     .. Scalar Arguments ..
 213       CHARACTER          COMPQ, COMPZ
 214       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
 215 *     ..
 216 *     .. Array Arguments ..
 217       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 218      $                   Z( LDZ, * )
 219 *     ..
 220 *
 221 *  =====================================================================
 222 *
 223 *     .. Parameters ..
 224       COMPLEX*16         CONE, CZERO
 225       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
 226      $                   CZERO = ( 0.0D+0, 0.0D+0 ) )
 227 *     ..
 228 *     .. Local Scalars ..
 229       LOGICAL            ILQ, ILZ
 230       INTEGER            ICOMPQ, ICOMPZ, JCOL, JROW
 231       DOUBLE PRECISION   C
 232       COMPLEX*16         CTEMP, S
 233 *     ..
 234 *     .. External Functions ..
 235       LOGICAL            LSAME
 236       EXTERNAL           LSAME
 237 *     ..
 238 *     .. External Subroutines ..
 239       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT
 240 *     ..
 241 *     .. Intrinsic Functions ..
 242       INTRINSIC          DCONJG, MAX
 243 *     ..
 244 *     .. Executable Statements ..
 245 *
 246 *     Decode COMPQ
 247 *
 248       IF( LSAME( COMPQ, 'N' ) ) THEN
 249          ILQ = .FALSE.
 250          ICOMPQ = 1
 251       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
 252          ILQ = .TRUE.
 253          ICOMPQ = 2
 254       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
 255          ILQ = .TRUE.
 256          ICOMPQ = 3
 257       ELSE
 258          ICOMPQ = 0
 259       END IF
 260 *
 261 *     Decode COMPZ
 262 *
 263       IF( LSAME( COMPZ, 'N' ) ) THEN
 264          ILZ = .FALSE.
 265          ICOMPZ = 1
 266       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 267          ILZ = .TRUE.
 268          ICOMPZ = 2
 269       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 270          ILZ = .TRUE.
 271          ICOMPZ = 3
 272       ELSE
 273          ICOMPZ = 0
 274       END IF
 275 *
 276 *     Test the input parameters.
 277 *
 278       INFO = 0
 279       IF( ICOMPQ.LE.0 ) THEN
 280          INFO = -1
 281       ELSE IF( ICOMPZ.LE.0 ) THEN
 282          INFO = -2
 283       ELSE IF( N.LT.0 ) THEN
 284          INFO = -3
 285       ELSE IF( ILO.LT.1 ) THEN
 286          INFO = -4
 287       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
 288          INFO = -5
 289       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 290          INFO = -7
 291       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 292          INFO = -9
 293       ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
 294          INFO = -11
 295       ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
 296          INFO = -13
 297       END IF
 298       IF( INFO.NE.0 ) THEN
 299          CALL XERBLA( 'ZGGHRD', -INFO )
 300          RETURN
 301       END IF
 302 *
 303 *     Initialize Q and Z if desired.
 304 *
 305       IF( ICOMPQ.EQ.3 )
 306      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
 307       IF( ICOMPZ.EQ.3 )
 308      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
 309 *
 310 *     Quick return if possible
 311 *
 312       IF( N.LE.1 )
 313      $   RETURN
 314 *
 315 *     Zero out lower triangle of B
 316 *
 317       DO 20 JCOL = 1, N - 1
 318          DO 10 JROW = JCOL + 1, N
 319             B( JROW, JCOL ) = CZERO
 320    10    CONTINUE
 321    20 CONTINUE
 322 *
 323 *     Reduce A and B
 324 *
 325       DO 40 JCOL = ILO, IHI - 2
 326 *
 327          DO 30 JROW = IHI, JCOL + 2, -1
 328 *
 329 *           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
 330 *
 331             CTEMP = A( JROW-1, JCOL )
 332             CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
 333      $                   A( JROW-1, JCOL ) )
 334             A( JROW, JCOL ) = CZERO
 335             CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
 336      $                 A( JROW, JCOL+1 ), LDA, C, S )
 337             CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
 338      $                 B( JROW, JROW-1 ), LDB, C, S )
 339             IF( ILQ )
 340      $         CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
 341      $                    DCONJG( S ) )
 342 *
 343 *           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
 344 *
 345             CTEMP = B( JROW, JROW )
 346             CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
 347      $                   B( JROW, JROW ) )
 348             B( JROW, JROW-1 ) = CZERO
 349             CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
 350             CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
 351      $                 S )
 352             IF( ILZ )
 353      $         CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
 354    30    CONTINUE
 355    40 CONTINUE
 356 *
 357       RETURN
 358 *
 359 *     End of ZGGHRD
 360 *
 361       END