00001 SUBROUTINE DGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, 00002 $ LDQ, Z, LDZ, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER COMPQ, COMPZ 00011 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00015 $ Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DGGHRD reduces a pair of real matrices (A,B) to generalized upper 00022 * Hessenberg form using orthogonal transformations, where A is a 00023 * general matrix and B is upper triangular. The form of the 00024 * generalized eigenvalue problem is 00025 * A*x = lambda*B*x, 00026 * and B is typically made upper triangular by computing its QR 00027 * factorization and moving the orthogonal matrix Q to the left side 00028 * of the equation. 00029 * 00030 * This subroutine simultaneously reduces A to a Hessenberg matrix H: 00031 * Q**T*A*Z = H 00032 * and transforms B to another upper triangular matrix T: 00033 * Q**T*B*Z = T 00034 * in order to reduce the problem to its standard form 00035 * H*y = lambda*T*y 00036 * where y = Z**T*x. 00037 * 00038 * The orthogonal matrices Q and Z are determined as products of Givens 00039 * rotations. They may either be formed explicitly, or they may be 00040 * postmultiplied into input matrices Q1 and Z1, so that 00041 * 00042 * Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T 00043 * 00044 * Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T 00045 * 00046 * If Q1 is the orthogonal matrix from the QR factorization of B in the 00047 * original equation A*x = lambda*B*x, then DGGHRD reduces the original 00048 * problem to generalized Hessenberg form. 00049 * 00050 * Arguments 00051 * ========= 00052 * 00053 * COMPQ (input) CHARACTER*1 00054 * = 'N': do not compute Q; 00055 * = 'I': Q is initialized to the unit matrix, and the 00056 * orthogonal matrix Q is returned; 00057 * = 'V': Q must contain an orthogonal matrix Q1 on entry, 00058 * and the product Q1*Q is returned. 00059 * 00060 * COMPZ (input) CHARACTER*1 00061 * = 'N': do not compute Z; 00062 * = 'I': Z is initialized to the unit matrix, and the 00063 * orthogonal matrix Z is returned; 00064 * = 'V': Z must contain an orthogonal matrix Z1 on entry, 00065 * and the product Z1*Z is returned. 00066 * 00067 * N (input) INTEGER 00068 * The order of the matrices A and B. N >= 0. 00069 * 00070 * ILO (input) INTEGER 00071 * IHI (input) INTEGER 00072 * ILO and IHI mark the rows and columns of A which are to be 00073 * reduced. It is assumed that A is already upper triangular 00074 * in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are 00075 * normally set by a previous call to SGGBAL; otherwise they 00076 * should be set to 1 and N respectively. 00077 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 00078 * 00079 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00080 * On entry, the N-by-N general matrix to be reduced. 00081 * On exit, the upper triangle and the first subdiagonal of A 00082 * are overwritten with the upper Hessenberg matrix H, and the 00083 * rest is set to zero. 00084 * 00085 * LDA (input) INTEGER 00086 * The leading dimension of the array A. LDA >= max(1,N). 00087 * 00088 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00089 * On entry, the N-by-N upper triangular matrix B. 00090 * On exit, the upper triangular matrix T = Q**T B Z. The 00091 * elements below the diagonal are set to zero. 00092 * 00093 * LDB (input) INTEGER 00094 * The leading dimension of the array B. LDB >= max(1,N). 00095 * 00096 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) 00097 * On entry, if COMPQ = 'V', the orthogonal matrix Q1, 00098 * typically from the QR factorization of B. 00099 * On exit, if COMPQ='I', the orthogonal matrix Q, and if 00100 * COMPQ = 'V', the product Q1*Q. 00101 * Not referenced if COMPQ='N'. 00102 * 00103 * LDQ (input) INTEGER 00104 * The leading dimension of the array Q. 00105 * LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. 00106 * 00107 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 00108 * On entry, if COMPZ = 'V', the orthogonal matrix Z1. 00109 * On exit, if COMPZ='I', the orthogonal matrix Z, and if 00110 * COMPZ = 'V', the product Z1*Z. 00111 * Not referenced if COMPZ='N'. 00112 * 00113 * LDZ (input) INTEGER 00114 * The leading dimension of the array Z. 00115 * LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. 00116 * 00117 * INFO (output) INTEGER 00118 * = 0: successful exit. 00119 * < 0: if INFO = -i, the i-th argument had an illegal value. 00120 * 00121 * Further Details 00122 * =============== 00123 * 00124 * This routine reduces A to Hessenberg and B to triangular form by 00125 * an unblocked reduction, as described in _Matrix_Computations_, 00126 * by Golub and Van Loan (Johns Hopkins Press.) 00127 * 00128 * ===================================================================== 00129 * 00130 * .. Parameters .. 00131 DOUBLE PRECISION ONE, ZERO 00132 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00133 * .. 00134 * .. Local Scalars .. 00135 LOGICAL ILQ, ILZ 00136 INTEGER ICOMPQ, ICOMPZ, JCOL, JROW 00137 DOUBLE PRECISION C, S, TEMP 00138 * .. 00139 * .. External Functions .. 00140 LOGICAL LSAME 00141 EXTERNAL LSAME 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL DLARTG, DLASET, DROT, XERBLA 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC MAX 00148 * .. 00149 * .. Executable Statements .. 00150 * 00151 * Decode COMPQ 00152 * 00153 IF( LSAME( COMPQ, 'N' ) ) THEN 00154 ILQ = .FALSE. 00155 ICOMPQ = 1 00156 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 00157 ILQ = .TRUE. 00158 ICOMPQ = 2 00159 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 00160 ILQ = .TRUE. 00161 ICOMPQ = 3 00162 ELSE 00163 ICOMPQ = 0 00164 END IF 00165 * 00166 * Decode COMPZ 00167 * 00168 IF( LSAME( COMPZ, 'N' ) ) THEN 00169 ILZ = .FALSE. 00170 ICOMPZ = 1 00171 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00172 ILZ = .TRUE. 00173 ICOMPZ = 2 00174 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00175 ILZ = .TRUE. 00176 ICOMPZ = 3 00177 ELSE 00178 ICOMPZ = 0 00179 END IF 00180 * 00181 * Test the input parameters. 00182 * 00183 INFO = 0 00184 IF( ICOMPQ.LE.0 ) THEN 00185 INFO = -1 00186 ELSE IF( ICOMPZ.LE.0 ) THEN 00187 INFO = -2 00188 ELSE IF( N.LT.0 ) THEN 00189 INFO = -3 00190 ELSE IF( ILO.LT.1 ) THEN 00191 INFO = -4 00192 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 00193 INFO = -5 00194 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00195 INFO = -7 00196 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00197 INFO = -9 00198 ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN 00199 INFO = -11 00200 ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN 00201 INFO = -13 00202 END IF 00203 IF( INFO.NE.0 ) THEN 00204 CALL XERBLA( 'DGGHRD', -INFO ) 00205 RETURN 00206 END IF 00207 * 00208 * Initialize Q and Z if desired. 00209 * 00210 IF( ICOMPQ.EQ.3 ) 00211 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00212 IF( ICOMPZ.EQ.3 ) 00213 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 00214 * 00215 * Quick return if possible 00216 * 00217 IF( N.LE.1 ) 00218 $ RETURN 00219 * 00220 * Zero out lower triangle of B 00221 * 00222 DO 20 JCOL = 1, N - 1 00223 DO 10 JROW = JCOL + 1, N 00224 B( JROW, JCOL ) = ZERO 00225 10 CONTINUE 00226 20 CONTINUE 00227 * 00228 * Reduce A and B 00229 * 00230 DO 40 JCOL = ILO, IHI - 2 00231 * 00232 DO 30 JROW = IHI, JCOL + 2, -1 00233 * 00234 * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) 00235 * 00236 TEMP = A( JROW-1, JCOL ) 00237 CALL DLARTG( TEMP, A( JROW, JCOL ), C, S, 00238 $ A( JROW-1, JCOL ) ) 00239 A( JROW, JCOL ) = ZERO 00240 CALL DROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA, 00241 $ A( JROW, JCOL+1 ), LDA, C, S ) 00242 CALL DROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB, 00243 $ B( JROW, JROW-1 ), LDB, C, S ) 00244 IF( ILQ ) 00245 $ CALL DROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, S ) 00246 * 00247 * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) 00248 * 00249 TEMP = B( JROW, JROW ) 00250 CALL DLARTG( TEMP, B( JROW, JROW-1 ), C, S, 00251 $ B( JROW, JROW ) ) 00252 B( JROW, JROW-1 ) = ZERO 00253 CALL DROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S ) 00254 CALL DROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C, 00255 $ S ) 00256 IF( ILZ ) 00257 $ CALL DROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S ) 00258 30 CONTINUE 00259 40 CONTINUE 00260 * 00261 RETURN 00262 * 00263 * End of DGGHRD 00264 * 00265 END
1.7.2