SRC/sgghrd.f

   1 *> \brief \b SGGHRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGGHRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgghrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgghrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgghrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGGHRD( 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 *       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  30 *      $                   Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SGGHRD reduces a pair of real matrices (A,B) to generalized upper
  40 *> Hessenberg form using orthogonal 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 orthogonal matrix Q to the left side
  46 *> of the equation.
  47 *>
  48 *> This subroutine simultaneously reduces A to a Hessenberg matrix H:
  49 *>    Q**T*A*Z = H
  50 *> and transforms B to another upper triangular matrix T:
  51 *>    Q**T*B*Z = T
  52 *> in order to reduce the problem to its standard form
  53 *>    H*y = lambda*T*y
  54 *> where y = Z**T*x.
  55 *>
  56 *> The orthogonal 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 *>
  60 *>      Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
  61 *>
  62 *>      Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
  63 *>
  64 *> If Q1 is the orthogonal matrix from the QR factorization of B in the
  65 *> original equation A*x = lambda*B*x, then SGGHRD reduces the original
  66 *> problem to generalized Hessenberg form.
  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 *>                 orthogonal matrix Q is returned;
  78 *>          = 'V': Q must contain an orthogonal 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 *>                 orthogonal matrix Z is returned;
  88 *>          = 'V': Z must contain an orthogonal 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 SGGBAL; 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 REAL 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 REAL array, dimension (LDB, N)
 133 *>          On entry, the N-by-N upper triangular matrix B.
 134 *>          On exit, the upper triangular matrix T = Q**T 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 REAL array, dimension (LDQ, N)
 147 *>          On entry, if COMPQ = 'V', the orthogonal matrix Q1,
 148 *>          typically from the QR factorization of B.
 149 *>          On exit, if COMPQ='I', the orthogonal 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 REAL array, dimension (LDZ, N)
 164 *>          On entry, if COMPZ = 'V', the orthogonal matrix Z1.
 165 *>          On exit, if COMPZ='I', the orthogonal 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] INFO
 178 *> \verbatim
 179 *>          INFO is INTEGER
 180 *>          = 0:  successful exit.
 181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 182 *> \endverbatim
 183 *
 184 *  Authors:
 185 *  ========
 186 *
 187 *> \author Univ. of Tennessee
 188 *> \author Univ. of California Berkeley
 189 *> \author Univ. of Colorado Denver
 190 *> \author NAG Ltd.
 191 *
 192 *> \date November 2011
 193 *
 194 *> \ingroup realOTHERcomputational
 195 *
 196 *> \par Further Details:
 197 *  =====================
 198 *>
 199 *> \verbatim
 200 *>
 201 *>  This routine reduces A to Hessenberg and B to triangular form by
 202 *>  an unblocked reduction, as described in _Matrix_Computations_,
 203 *>  by Golub and Van Loan (Johns Hopkins Press.)
 204 *> \endverbatim
 205 *>
 206 *  =====================================================================
 207       SUBROUTINE SGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
 208      $                   LDQ, Z, LDZ, INFO )
 209 *
 210 *  -- LAPACK computational routine (version 3.4.0) --
 211 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 212 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 213 *     November 2011
 214 *
 215 *     .. Scalar Arguments ..
 216       CHARACTER          COMPQ, COMPZ
 217       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
 218 *     ..
 219 *     .. Array Arguments ..
 220       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 221      $                   Z( LDZ, * )
 222 *     ..
 223 *
 224 *  =====================================================================
 225 *
 226 *     .. Parameters ..
 227       REAL               ONE, ZERO
 228       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 229 *     ..
 230 *     .. Local Scalars ..
 231       LOGICAL            ILQ, ILZ
 232       INTEGER            ICOMPQ, ICOMPZ, JCOL, JROW
 233       REAL               C, S, TEMP
 234 *     ..
 235 *     .. External Functions ..
 236       LOGICAL            LSAME
 237       EXTERNAL           LSAME
 238 *     ..
 239 *     .. External Subroutines ..
 240       EXTERNAL           SLARTG, SLASET, SROT, XERBLA
 241 *     ..
 242 *     .. Intrinsic Functions ..
 243       INTRINSIC          MAX
 244 *     ..
 245 *     .. Executable Statements ..
 246 *
 247 *     Decode COMPQ
 248 *
 249       IF( LSAME( COMPQ, 'N' ) ) THEN
 250          ILQ = .FALSE.
 251          ICOMPQ = 1
 252       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
 253          ILQ = .TRUE.
 254          ICOMPQ = 2
 255       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
 256          ILQ = .TRUE.
 257          ICOMPQ = 3
 258       ELSE
 259          ICOMPQ = 0
 260       END IF
 261 *
 262 *     Decode COMPZ
 263 *
 264       IF( LSAME( COMPZ, 'N' ) ) THEN
 265          ILZ = .FALSE.
 266          ICOMPZ = 1
 267       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 268          ILZ = .TRUE.
 269          ICOMPZ = 2
 270       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 271          ILZ = .TRUE.
 272          ICOMPZ = 3
 273       ELSE
 274          ICOMPZ = 0
 275       END IF
 276 *
 277 *     Test the input parameters.
 278 *
 279       INFO = 0
 280       IF( ICOMPQ.LE.0 ) THEN
 281          INFO = -1
 282       ELSE IF( ICOMPZ.LE.0 ) THEN
 283          INFO = -2
 284       ELSE IF( N.LT.0 ) THEN
 285          INFO = -3
 286       ELSE IF( ILO.LT.1 ) THEN
 287          INFO = -4
 288       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
 289          INFO = -5
 290       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 291          INFO = -7
 292       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 293          INFO = -9
 294       ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
 295          INFO = -11
 296       ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
 297          INFO = -13
 298       END IF
 299       IF( INFO.NE.0 ) THEN
 300          CALL XERBLA( 'SGGHRD', -INFO )
 301          RETURN
 302       END IF
 303 *
 304 *     Initialize Q and Z if desired.
 305 *
 306       IF( ICOMPQ.EQ.3 )
 307      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
 308       IF( ICOMPZ.EQ.3 )
 309      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
 310 *
 311 *     Quick return if possible
 312 *
 313       IF( N.LE.1 )
 314      $   RETURN
 315 *
 316 *     Zero out lower triangle of B
 317 *
 318       DO 20 JCOL = 1, N - 1
 319          DO 10 JROW = JCOL + 1, N
 320             B( JROW, JCOL ) = ZERO
 321    10    CONTINUE
 322    20 CONTINUE
 323 *
 324 *     Reduce A and B
 325 *
 326       DO 40 JCOL = ILO, IHI - 2
 327 *
 328          DO 30 JROW = IHI, JCOL + 2, -1
 329 *
 330 *           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
 331 *
 332             TEMP = A( JROW-1, JCOL )
 333             CALL SLARTG( TEMP, A( JROW, JCOL ), C, S,
 334      $                   A( JROW-1, JCOL ) )
 335             A( JROW, JCOL ) = ZERO
 336             CALL SROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
 337      $                 A( JROW, JCOL+1 ), LDA, C, S )
 338             CALL SROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
 339      $                 B( JROW, JROW-1 ), LDB, C, S )
 340             IF( ILQ )
 341      $         CALL SROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, S )
 342 *
 343 *           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
 344 *
 345             TEMP = B( JROW, JROW )
 346             CALL SLARTG( TEMP, B( JROW, JROW-1 ), C, S,
 347      $                   B( JROW, JROW ) )
 348             B( JROW, JROW-1 ) = ZERO
 349             CALL SROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
 350             CALL SROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
 351      $                 S )
 352             IF( ILZ )
 353      $         CALL SROT( 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 SGGHRD
 360 *
 361       END