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