00001 /* dgghrd.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static doublereal c_b10 = 0.; 00019 static doublereal c_b11 = 1.; 00020 static integer c__1 = 1; 00021 00022 /* Subroutine */ int dgghrd_(char *compq, char *compz, integer *n, integer * 00023 ilo, integer *ihi, doublereal *a, integer *lda, doublereal *b, 00024 integer *ldb, doublereal *q, integer *ldq, doublereal *z__, integer * 00025 ldz, integer *info) 00026 { 00027 /* System generated locals */ 00028 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 00029 z_offset, i__1, i__2, i__3; 00030 00031 /* Local variables */ 00032 doublereal c__, s; 00033 logical ilq, ilz; 00034 integer jcol; 00035 doublereal temp; 00036 extern /* Subroutine */ int drot_(integer *, doublereal *, integer *, 00037 doublereal *, integer *, doublereal *, doublereal *); 00038 integer jrow; 00039 extern logical lsame_(char *, char *); 00040 extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 00041 doublereal *, doublereal *, doublereal *, integer *), 00042 dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, 00043 doublereal *), xerbla_(char *, integer *); 00044 integer icompq, icompz; 00045 00046 00047 /* -- LAPACK routine (version 3.2) -- */ 00048 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00049 /* November 2006 */ 00050 00051 /* .. Scalar Arguments .. */ 00052 /* .. */ 00053 /* .. Array Arguments .. */ 00054 /* .. */ 00055 00056 /* Purpose */ 00057 /* ======= */ 00058 00059 /* DGGHRD reduces a pair of real matrices (A,B) to generalized upper */ 00060 /* Hessenberg form using orthogonal transformations, where A is a */ 00061 /* general matrix and B is upper triangular. The form of the */ 00062 /* generalized eigenvalue problem is */ 00063 /* A*x = lambda*B*x, */ 00064 /* and B is typically made upper triangular by computing its QR */ 00065 /* factorization and moving the orthogonal matrix Q to the left side */ 00066 /* of the equation. */ 00067 00068 /* This subroutine simultaneously reduces A to a Hessenberg matrix H: */ 00069 /* Q**T*A*Z = H */ 00070 /* and transforms B to another upper triangular matrix T: */ 00071 /* Q**T*B*Z = T */ 00072 /* in order to reduce the problem to its standard form */ 00073 /* H*y = lambda*T*y */ 00074 /* where y = Z**T*x. */ 00075 00076 /* The orthogonal matrices Q and Z are determined as products of Givens */ 00077 /* rotations. They may either be formed explicitly, or they may be */ 00078 /* postmultiplied into input matrices Q1 and Z1, so that */ 00079 00080 /* Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T */ 00081 00082 /* Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T */ 00083 00084 /* If Q1 is the orthogonal matrix from the QR factorization of B in the */ 00085 /* original equation A*x = lambda*B*x, then DGGHRD reduces the original */ 00086 /* problem to generalized Hessenberg form. */ 00087 00088 /* Arguments */ 00089 /* ========= */ 00090 00091 /* COMPQ (input) CHARACTER*1 */ 00092 /* = 'N': do not compute Q; */ 00093 /* = 'I': Q is initialized to the unit matrix, and the */ 00094 /* orthogonal matrix Q is returned; */ 00095 /* = 'V': Q must contain an orthogonal matrix Q1 on entry, */ 00096 /* and the product Q1*Q is returned. */ 00097 00098 /* COMPZ (input) CHARACTER*1 */ 00099 /* = 'N': do not compute Z; */ 00100 /* = 'I': Z is initialized to the unit matrix, and the */ 00101 /* orthogonal matrix Z is returned; */ 00102 /* = 'V': Z must contain an orthogonal matrix Z1 on entry, */ 00103 /* and the product Z1*Z is returned. */ 00104 00105 /* N (input) INTEGER */ 00106 /* The order of the matrices A and B. N >= 0. */ 00107 00108 /* ILO (input) INTEGER */ 00109 /* IHI (input) INTEGER */ 00110 /* ILO and IHI mark the rows and columns of A which are to be */ 00111 /* reduced. It is assumed that A is already upper triangular */ 00112 /* in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are */ 00113 /* normally set by a previous call to SGGBAL; otherwise they */ 00114 /* should be set to 1 and N respectively. */ 00115 /* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */ 00116 00117 /* A (input/output) 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 00123 /* LDA (input) INTEGER */ 00124 /* The leading dimension of the array A. LDA >= max(1,N). */ 00125 00126 /* B (input/output) DOUBLE PRECISION array, dimension (LDB, N) */ 00127 /* On entry, the N-by-N upper triangular matrix B. */ 00128 /* On exit, the upper triangular matrix T = Q**T B Z. The */ 00129 /* elements below the diagonal are set to zero. */ 00130 00131 /* LDB (input) INTEGER */ 00132 /* The leading dimension of the array B. LDB >= max(1,N). */ 00133 00134 /* Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) */ 00135 /* On entry, if COMPQ = 'V', the orthogonal matrix Q1, */ 00136 /* typically from the QR factorization of B. */ 00137 /* On exit, if COMPQ='I', the orthogonal matrix Q, and if */ 00138 /* COMPQ = 'V', the product Q1*Q. */ 00139 /* Not referenced if COMPQ='N'. */ 00140 00141 /* LDQ (input) INTEGER */ 00142 /* The leading dimension of the array Q. */ 00143 /* LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. */ 00144 00145 /* Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) */ 00146 /* On entry, if COMPZ = 'V', the orthogonal matrix Z1. */ 00147 /* On exit, if COMPZ='I', the orthogonal matrix Z, and if */ 00148 /* COMPZ = 'V', the product Z1*Z. */ 00149 /* Not referenced if COMPZ='N'. */ 00150 00151 /* LDZ (input) INTEGER */ 00152 /* The leading dimension of the array Z. */ 00153 /* LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. */ 00154 00155 /* INFO (output) INTEGER */ 00156 /* = 0: successful exit. */ 00157 /* < 0: if INFO = -i, the i-th argument had an illegal value. */ 00158 00159 /* Further Details */ 00160 /* =============== */ 00161 00162 /* This routine reduces A to Hessenberg and B to triangular form by */ 00163 /* an unblocked reduction, as described in _Matrix_Computations_, */ 00164 /* by Golub and Van Loan (Johns Hopkins Press.) */ 00165 00166 /* ===================================================================== */ 00167 00168 /* .. Parameters .. */ 00169 /* .. */ 00170 /* .. Local Scalars .. */ 00171 /* .. */ 00172 /* .. External Functions .. */ 00173 /* .. */ 00174 /* .. External Subroutines .. */ 00175 /* .. */ 00176 /* .. Intrinsic Functions .. */ 00177 /* .. */ 00178 /* .. Executable Statements .. */ 00179 00180 /* Decode COMPQ */ 00181 00182 /* Parameter adjustments */ 00183 a_dim1 = *lda; 00184 a_offset = 1 + a_dim1; 00185 a -= a_offset; 00186 b_dim1 = *ldb; 00187 b_offset = 1 + b_dim1; 00188 b -= b_offset; 00189 q_dim1 = *ldq; 00190 q_offset = 1 + q_dim1; 00191 q -= q_offset; 00192 z_dim1 = *ldz; 00193 z_offset = 1 + z_dim1; 00194 z__ -= z_offset; 00195 00196 /* Function Body */ 00197 if (lsame_(compq, "N")) { 00198 ilq = FALSE_; 00199 icompq = 1; 00200 } else if (lsame_(compq, "V")) { 00201 ilq = TRUE_; 00202 icompq = 2; 00203 } else if (lsame_(compq, "I")) { 00204 ilq = TRUE_; 00205 icompq = 3; 00206 } else { 00207 icompq = 0; 00208 } 00209 00210 /* Decode COMPZ */ 00211 00212 if (lsame_(compz, "N")) { 00213 ilz = FALSE_; 00214 icompz = 1; 00215 } else if (lsame_(compz, "V")) { 00216 ilz = TRUE_; 00217 icompz = 2; 00218 } else if (lsame_(compz, "I")) { 00219 ilz = TRUE_; 00220 icompz = 3; 00221 } else { 00222 icompz = 0; 00223 } 00224 00225 /* Test the input parameters. */ 00226 00227 *info = 0; 00228 if (icompq <= 0) { 00229 *info = -1; 00230 } else if (icompz <= 0) { 00231 *info = -2; 00232 } else if (*n < 0) { 00233 *info = -3; 00234 } else if (*ilo < 1) { 00235 *info = -4; 00236 } else if (*ihi > *n || *ihi < *ilo - 1) { 00237 *info = -5; 00238 } else if (*lda < max(1,*n)) { 00239 *info = -7; 00240 } else if (*ldb < max(1,*n)) { 00241 *info = -9; 00242 } else if (ilq && *ldq < *n || *ldq < 1) { 00243 *info = -11; 00244 } else if (ilz && *ldz < *n || *ldz < 1) { 00245 *info = -13; 00246 } 00247 if (*info != 0) { 00248 i__1 = -(*info); 00249 xerbla_("DGGHRD", &i__1); 00250 return 0; 00251 } 00252 00253 /* Initialize Q and Z if desired. */ 00254 00255 if (icompq == 3) { 00256 dlaset_("Full", n, n, &c_b10, &c_b11, &q[q_offset], ldq); 00257 } 00258 if (icompz == 3) { 00259 dlaset_("Full", n, n, &c_b10, &c_b11, &z__[z_offset], ldz); 00260 } 00261 00262 /* Quick return if possible */ 00263 00264 if (*n <= 1) { 00265 return 0; 00266 } 00267 00268 /* Zero out lower triangle of B */ 00269 00270 i__1 = *n - 1; 00271 for (jcol = 1; jcol <= i__1; ++jcol) { 00272 i__2 = *n; 00273 for (jrow = jcol + 1; jrow <= i__2; ++jrow) { 00274 b[jrow + jcol * b_dim1] = 0.; 00275 /* L10: */ 00276 } 00277 /* L20: */ 00278 } 00279 00280 /* Reduce A and B */ 00281 00282 i__1 = *ihi - 2; 00283 for (jcol = *ilo; jcol <= i__1; ++jcol) { 00284 00285 i__2 = jcol + 2; 00286 for (jrow = *ihi; jrow >= i__2; --jrow) { 00287 00288 /* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) */ 00289 00290 temp = a[jrow - 1 + jcol * a_dim1]; 00291 dlartg_(&temp, &a[jrow + jcol * a_dim1], &c__, &s, &a[jrow - 1 + 00292 jcol * a_dim1]); 00293 a[jrow + jcol * a_dim1] = 0.; 00294 i__3 = *n - jcol; 00295 drot_(&i__3, &a[jrow - 1 + (jcol + 1) * a_dim1], lda, &a[jrow + ( 00296 jcol + 1) * a_dim1], lda, &c__, &s); 00297 i__3 = *n + 2 - jrow; 00298 drot_(&i__3, &b[jrow - 1 + (jrow - 1) * b_dim1], ldb, &b[jrow + ( 00299 jrow - 1) * b_dim1], ldb, &c__, &s); 00300 if (ilq) { 00301 drot_(n, &q[(jrow - 1) * q_dim1 + 1], &c__1, &q[jrow * q_dim1 00302 + 1], &c__1, &c__, &s); 00303 } 00304 00305 /* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) */ 00306 00307 temp = b[jrow + jrow * b_dim1]; 00308 dlartg_(&temp, &b[jrow + (jrow - 1) * b_dim1], &c__, &s, &b[jrow 00309 + jrow * b_dim1]); 00310 b[jrow + (jrow - 1) * b_dim1] = 0.; 00311 drot_(ihi, &a[jrow * a_dim1 + 1], &c__1, &a[(jrow - 1) * a_dim1 + 00312 1], &c__1, &c__, &s); 00313 i__3 = jrow - 1; 00314 drot_(&i__3, &b[jrow * b_dim1 + 1], &c__1, &b[(jrow - 1) * b_dim1 00315 + 1], &c__1, &c__, &s); 00316 if (ilz) { 00317 drot_(n, &z__[jrow * z_dim1 + 1], &c__1, &z__[(jrow - 1) * 00318 z_dim1 + 1], &c__1, &c__, &s); 00319 } 00320 /* L30: */ 00321 } 00322 /* L40: */ 00323 } 00324 00325 return 0; 00326 00327 /* End of DGGHRD */ 00328 00329 } /* dgghrd_ */