LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgghrd.f
Go to the documentation of this file.
00001 *> \brief \b DGGHRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DGGHRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgghrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgghrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgghrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
00022 *                          LDQ, Z, LDZ, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          COMPQ, COMPZ
00026 *       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00030 *      $                   Z( LDZ, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DGGHRD reduces a pair of real matrices (A,B) to generalized upper
00040 *> Hessenberg form using orthogonal transformations, where A is a
00041 *> general matrix and B is upper triangular.  The form of the
00042 *> generalized eigenvalue problem is
00043 *>    A*x = lambda*B*x,
00044 *> and B is typically made upper triangular by computing its QR
00045 *> factorization and moving the orthogonal matrix Q to the left side
00046 *> of the equation.
00047 *>
00048 *> This subroutine simultaneously reduces A to a Hessenberg matrix H:
00049 *>    Q**T*A*Z = H
00050 *> and transforms B to another upper triangular matrix T:
00051 *>    Q**T*B*Z = T
00052 *> in order to reduce the problem to its standard form
00053 *>    H*y = lambda*T*y
00054 *> where y = Z**T*x.
00055 *>
00056 *> The orthogonal matrices Q and Z are determined as products of Givens
00057 *> rotations.  They may either be formed explicitly, or they may be
00058 *> postmultiplied into input matrices Q1 and Z1, so that
00059 *>
00060 *>      Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
00061 *>
00062 *>      Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
00063 *>
00064 *> If Q1 is the orthogonal matrix from the QR factorization of B in the
00065 *> original equation A*x = lambda*B*x, then DGGHRD reduces the original
00066 *> problem to generalized Hessenberg form.
00067 *> \endverbatim
00068 *
00069 *  Arguments:
00070 *  ==========
00071 *
00072 *> \param[in] COMPQ
00073 *> \verbatim
00074 *>          COMPQ is CHARACTER*1
00075 *>          = 'N': do not compute Q;
00076 *>          = 'I': Q is initialized to the unit matrix, and the
00077 *>                 orthogonal matrix Q is returned;
00078 *>          = 'V': Q must contain an orthogonal matrix Q1 on entry,
00079 *>                 and the product Q1*Q is returned.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] COMPZ
00083 *> \verbatim
00084 *>          COMPZ is CHARACTER*1
00085 *>          = 'N': do not compute Z;
00086 *>          = 'I': Z is initialized to the unit matrix, and the
00087 *>                 orthogonal matrix Z is returned;
00088 *>          = 'V': Z must contain an orthogonal matrix Z1 on entry,
00089 *>                 and the product Z1*Z is returned.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] N
00093 *> \verbatim
00094 *>          N is INTEGER
00095 *>          The order of the matrices A and B.  N >= 0.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] ILO
00099 *> \verbatim
00100 *>          ILO is INTEGER
00101 *> \endverbatim
00102 *>
00103 *> \param[in] IHI
00104 *> \verbatim
00105 *>          IHI is INTEGER
00106 *>
00107 *>          ILO and IHI mark the rows and columns of A which are to be
00108 *>          reduced.  It is assumed that A is already upper triangular
00109 *>          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
00110 *>          normally set by a previous call to DGGBAL; otherwise they
00111 *>          should be set to 1 and N respectively.
00112 *>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] A
00116 *> \verbatim
00117 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
00118 *>          On entry, the N-by-N general matrix to be reduced.
00119 *>          On exit, the upper triangle and the first subdiagonal of A
00120 *>          are overwritten with the upper Hessenberg matrix H, and the
00121 *>          rest is set to zero.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDA
00125 *> \verbatim
00126 *>          LDA is INTEGER
00127 *>          The leading dimension of the array A.  LDA >= max(1,N).
00128 *> \endverbatim
00129 *>
00130 *> \param[in,out] B
00131 *> \verbatim
00132 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
00133 *>          On entry, the N-by-N upper triangular matrix B.
00134 *>          On exit, the upper triangular matrix T = Q**T B Z.  The
00135 *>          elements below the diagonal are set to zero.
00136 *> \endverbatim
00137 *>
00138 *> \param[in] LDB
00139 *> \verbatim
00140 *>          LDB is INTEGER
00141 *>          The leading dimension of the array B.  LDB >= max(1,N).
00142 *> \endverbatim
00143 *>
00144 *> \param[in,out] Q
00145 *> \verbatim
00146 *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
00147 *>          On entry, if COMPQ = 'V', the orthogonal matrix Q1,
00148 *>          typically from the QR factorization of B.
00149 *>          On exit, if COMPQ='I', the orthogonal matrix Q, and if
00150 *>          COMPQ = 'V', the product Q1*Q.
00151 *>          Not referenced if COMPQ='N'.
00152 *> \endverbatim
00153 *>
00154 *> \param[in] LDQ
00155 *> \verbatim
00156 *>          LDQ is INTEGER
00157 *>          The leading dimension of the array Q.
00158 *>          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
00159 *> \endverbatim
00160 *>
00161 *> \param[in,out] Z
00162 *> \verbatim
00163 *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
00164 *>          On entry, if COMPZ = 'V', the orthogonal matrix Z1.
00165 *>          On exit, if COMPZ='I', the orthogonal matrix Z, and if
00166 *>          COMPZ = 'V', the product Z1*Z.
00167 *>          Not referenced if COMPZ='N'.
00168 *> \endverbatim
00169 *>
00170 *> \param[in] LDZ
00171 *> \verbatim
00172 *>          LDZ is INTEGER
00173 *>          The leading dimension of the array Z.
00174 *>          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
00175 *> \endverbatim
00176 *>
00177 *> \param[out] INFO
00178 *> \verbatim
00179 *>          INFO is INTEGER
00180 *>          = 0:  successful exit.
00181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00182 *> \endverbatim
00183 *
00184 *  Authors:
00185 *  ========
00186 *
00187 *> \author Univ. of Tennessee 
00188 *> \author Univ. of California Berkeley 
00189 *> \author Univ. of Colorado Denver 
00190 *> \author NAG Ltd. 
00191 *
00192 *> \date November 2011
00193 *
00194 *> \ingroup doubleOTHERcomputational
00195 *
00196 *> \par Further Details:
00197 *  =====================
00198 *>
00199 *> \verbatim
00200 *>
00201 *>  This routine reduces A to Hessenberg and B to triangular form by
00202 *>  an unblocked reduction, as described in _Matrix_Computations_,
00203 *>  by Golub and Van Loan (Johns Hopkins Press.)
00204 *> \endverbatim
00205 *>
00206 *  =====================================================================
00207       SUBROUTINE DGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
00208      $                   LDQ, Z, LDZ, INFO )
00209 *
00210 *  -- LAPACK computational routine (version 3.4.0) --
00211 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00212 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00213 *     November 2011
00214 *
00215 *     .. Scalar Arguments ..
00216       CHARACTER          COMPQ, COMPZ
00217       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
00218 *     ..
00219 *     .. Array Arguments ..
00220       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00221      $                   Z( LDZ, * )
00222 *     ..
00223 *
00224 *  =====================================================================
00225 *
00226 *     .. Parameters ..
00227       DOUBLE PRECISION   ONE, ZERO
00228       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00229 *     ..
00230 *     .. Local Scalars ..
00231       LOGICAL            ILQ, ILZ
00232       INTEGER            ICOMPQ, ICOMPZ, JCOL, JROW
00233       DOUBLE PRECISION   C, S, TEMP
00234 *     ..
00235 *     .. External Functions ..
00236       LOGICAL            LSAME
00237       EXTERNAL           LSAME
00238 *     ..
00239 *     .. External Subroutines ..
00240       EXTERNAL           DLARTG, DLASET, DROT, XERBLA
00241 *     ..
00242 *     .. Intrinsic Functions ..
00243       INTRINSIC          MAX
00244 *     ..
00245 *     .. Executable Statements ..
00246 *
00247 *     Decode COMPQ
00248 *
00249       IF( LSAME( COMPQ, 'N' ) ) THEN
00250          ILQ = .FALSE.
00251          ICOMPQ = 1
00252       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
00253          ILQ = .TRUE.
00254          ICOMPQ = 2
00255       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
00256          ILQ = .TRUE.
00257          ICOMPQ = 3
00258       ELSE
00259          ICOMPQ = 0
00260       END IF
00261 *
00262 *     Decode COMPZ
00263 *
00264       IF( LSAME( COMPZ, 'N' ) ) THEN
00265          ILZ = .FALSE.
00266          ICOMPZ = 1
00267       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00268          ILZ = .TRUE.
00269          ICOMPZ = 2
00270       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00271          ILZ = .TRUE.
00272          ICOMPZ = 3
00273       ELSE
00274          ICOMPZ = 0
00275       END IF
00276 *
00277 *     Test the input parameters.
00278 *
00279       INFO = 0
00280       IF( ICOMPQ.LE.0 ) THEN
00281          INFO = -1
00282       ELSE IF( ICOMPZ.LE.0 ) THEN
00283          INFO = -2
00284       ELSE IF( N.LT.0 ) THEN
00285          INFO = -3
00286       ELSE IF( ILO.LT.1 ) THEN
00287          INFO = -4
00288       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
00289          INFO = -5
00290       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00291          INFO = -7
00292       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00293          INFO = -9
00294       ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
00295          INFO = -11
00296       ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
00297          INFO = -13
00298       END IF
00299       IF( INFO.NE.0 ) THEN
00300          CALL XERBLA( 'DGGHRD', -INFO )
00301          RETURN
00302       END IF
00303 *
00304 *     Initialize Q and Z if desired.
00305 *
00306       IF( ICOMPQ.EQ.3 )
00307      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00308       IF( ICOMPZ.EQ.3 )
00309      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
00310 *
00311 *     Quick return if possible
00312 *
00313       IF( N.LE.1 )
00314      $   RETURN
00315 *
00316 *     Zero out lower triangle of B
00317 *
00318       DO 20 JCOL = 1, N - 1
00319          DO 10 JROW = JCOL + 1, N
00320             B( JROW, JCOL ) = ZERO
00321    10    CONTINUE
00322    20 CONTINUE
00323 *
00324 *     Reduce A and B
00325 *
00326       DO 40 JCOL = ILO, IHI - 2
00327 *
00328          DO 30 JROW = IHI, JCOL + 2, -1
00329 *
00330 *           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
00331 *
00332             TEMP = A( JROW-1, JCOL )
00333             CALL DLARTG( TEMP, A( JROW, JCOL ), C, S,
00334      $                   A( JROW-1, JCOL ) )
00335             A( JROW, JCOL ) = ZERO
00336             CALL DROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
00337      $                 A( JROW, JCOL+1 ), LDA, C, S )
00338             CALL DROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
00339      $                 B( JROW, JROW-1 ), LDB, C, S )
00340             IF( ILQ )
00341      $         CALL DROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, S )
00342 *
00343 *           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
00344 *
00345             TEMP = B( JROW, JROW )
00346             CALL DLARTG( TEMP, B( JROW, JROW-1 ), C, S,
00347      $                   B( JROW, JROW ) )
00348             B( JROW, JROW-1 ) = ZERO
00349             CALL DROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
00350             CALL DROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
00351      $                 S )
00352             IF( ILZ )
00353      $         CALL DROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
00354    30    CONTINUE
00355    40 CONTINUE
00356 *
00357       RETURN
00358 *
00359 *     End of DGGHRD
00360 *
00361       END
 All Files Functions