Doxygen Source Code Documentation
eis_qzhes.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | qzhes_ (integer *nm, integer *n, doublereal *a, doublereal *b, logical *matz, doublereal *z__) |
Function Documentation
|
Definition at line 8 of file eis_qzhes.c. References a, abs, d_sign(), l, and v1. Referenced by rgg_().
00010 { 00011 /* System generated locals */ 00012 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2, 00013 i__3; 00014 doublereal d__1, d__2; 00015 00016 /* Builtin functions */ 00017 double sqrt(doublereal), d_sign(doublereal *, doublereal *); 00018 00019 /* Local variables */ 00020 static integer i__, j, k, l; 00021 static doublereal r__, s, t; 00022 static integer l1; 00023 static doublereal u1, u2, v1, v2; 00024 static integer lb, nk1, nm1, nm2; 00025 static doublereal rho; 00026 00027 00028 00029 /* THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM */ 00030 /* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */ 00031 /* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */ 00032 00033 /* THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND */ 00034 /* REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER */ 00035 /* TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. */ 00036 /* IT IS USUALLY FOLLOWED BY QZIT, QZVAL AND, POSSIBLY, QZVEC. */ 00037 00038 /* ON INPUT */ 00039 00040 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00041 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00042 /* DIMENSION STATEMENT. */ 00043 00044 /* N IS THE ORDER OF THE MATRICES. */ 00045 00046 /* A CONTAINS A REAL GENERAL MATRIX. */ 00047 00048 /* B CONTAINS A REAL GENERAL MATRIX. */ 00049 00050 /* MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 00051 */ 00052 /* ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */ 00053 /* EIGENVECTORS, AND TO .FALSE. OTHERWISE. */ 00054 00055 /* ON OUTPUT */ 00056 00057 /* A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS */ 00058 /* BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO. */ 00059 00060 /* B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS */ 00061 /* BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO. */ 00062 00063 /* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF */ 00064 /* MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z IS NOT REFERENCED. 00065 */ 00066 00067 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00068 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00069 */ 00070 00071 /* THIS VERSION DATED AUGUST 1983. */ 00072 00073 /* ------------------------------------------------------------------ 00074 */ 00075 00076 /* .......... INITIALIZE Z .......... */ 00077 /* Parameter adjustments */ 00078 z_dim1 = *nm; 00079 z_offset = z_dim1 + 1; 00080 z__ -= z_offset; 00081 b_dim1 = *nm; 00082 b_offset = b_dim1 + 1; 00083 b -= b_offset; 00084 a_dim1 = *nm; 00085 a_offset = a_dim1 + 1; 00086 a -= a_offset; 00087 00088 /* Function Body */ 00089 if (! (*matz)) { 00090 goto L10; 00091 } 00092 00093 i__1 = *n; 00094 for (j = 1; j <= i__1; ++j) { 00095 00096 i__2 = *n; 00097 for (i__ = 1; i__ <= i__2; ++i__) { 00098 z__[i__ + j * z_dim1] = 0.; 00099 /* L2: */ 00100 } 00101 00102 z__[j + j * z_dim1] = 1.; 00103 /* L3: */ 00104 } 00105 /* .......... REDUCE B TO UPPER TRIANGULAR FORM .......... */ 00106 L10: 00107 if (*n <= 1) { 00108 goto L170; 00109 } 00110 nm1 = *n - 1; 00111 00112 i__1 = nm1; 00113 for (l = 1; l <= i__1; ++l) { 00114 l1 = l + 1; 00115 s = 0.; 00116 00117 i__2 = *n; 00118 for (i__ = l1; i__ <= i__2; ++i__) { 00119 s += (d__1 = b[i__ + l * b_dim1], abs(d__1)); 00120 /* L20: */ 00121 } 00122 00123 if (s == 0.) { 00124 goto L100; 00125 } 00126 s += (d__1 = b[l + l * b_dim1], abs(d__1)); 00127 r__ = 0.; 00128 00129 i__2 = *n; 00130 for (i__ = l; i__ <= i__2; ++i__) { 00131 b[i__ + l * b_dim1] /= s; 00132 /* Computing 2nd power */ 00133 d__1 = b[i__ + l * b_dim1]; 00134 r__ += d__1 * d__1; 00135 /* L25: */ 00136 } 00137 00138 d__1 = sqrt(r__); 00139 r__ = d_sign(&d__1, &b[l + l * b_dim1]); 00140 b[l + l * b_dim1] += r__; 00141 rho = r__ * b[l + l * b_dim1]; 00142 00143 i__2 = *n; 00144 for (j = l1; j <= i__2; ++j) { 00145 t = 0.; 00146 00147 i__3 = *n; 00148 for (i__ = l; i__ <= i__3; ++i__) { 00149 t += b[i__ + l * b_dim1] * b[i__ + j * b_dim1]; 00150 /* L30: */ 00151 } 00152 00153 t = -t / rho; 00154 00155 i__3 = *n; 00156 for (i__ = l; i__ <= i__3; ++i__) { 00157 b[i__ + j * b_dim1] += t * b[i__ + l * b_dim1]; 00158 /* L40: */ 00159 } 00160 00161 /* L50: */ 00162 } 00163 00164 i__2 = *n; 00165 for (j = 1; j <= i__2; ++j) { 00166 t = 0.; 00167 00168 i__3 = *n; 00169 for (i__ = l; i__ <= i__3; ++i__) { 00170 t += b[i__ + l * b_dim1] * a[i__ + j * a_dim1]; 00171 /* L60: */ 00172 } 00173 00174 t = -t / rho; 00175 00176 i__3 = *n; 00177 for (i__ = l; i__ <= i__3; ++i__) { 00178 a[i__ + j * a_dim1] += t * b[i__ + l * b_dim1]; 00179 /* L70: */ 00180 } 00181 00182 /* L80: */ 00183 } 00184 00185 b[l + l * b_dim1] = -s * r__; 00186 00187 i__2 = *n; 00188 for (i__ = l1; i__ <= i__2; ++i__) { 00189 b[i__ + l * b_dim1] = 0.; 00190 /* L90: */ 00191 } 00192 00193 L100: 00194 ; 00195 } 00196 /* .......... REDUCE A TO UPPER HESSENBERG FORM, WHILE */ 00197 /* KEEPING B TRIANGULAR .......... */ 00198 if (*n == 2) { 00199 goto L170; 00200 } 00201 nm2 = *n - 2; 00202 00203 i__1 = nm2; 00204 for (k = 1; k <= i__1; ++k) { 00205 nk1 = nm1 - k; 00206 /* .......... FOR L=N-1 STEP -1 UNTIL K+1 DO -- .......... */ 00207 i__2 = nk1; 00208 for (lb = 1; lb <= i__2; ++lb) { 00209 l = *n - lb; 00210 l1 = l + 1; 00211 /* .......... ZERO A(L+1,K) .......... */ 00212 s = (d__1 = a[l + k * a_dim1], abs(d__1)) + (d__2 = a[l1 + k * 00213 a_dim1], abs(d__2)); 00214 if (s == 0.) { 00215 goto L150; 00216 } 00217 u1 = a[l + k * a_dim1] / s; 00218 u2 = a[l1 + k * a_dim1] / s; 00219 d__1 = sqrt(u1 * u1 + u2 * u2); 00220 r__ = d_sign(&d__1, &u1); 00221 v1 = -(u1 + r__) / r__; 00222 v2 = -u2 / r__; 00223 u2 = v2 / v1; 00224 00225 i__3 = *n; 00226 for (j = k; j <= i__3; ++j) { 00227 t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1]; 00228 a[l + j * a_dim1] += t * v1; 00229 a[l1 + j * a_dim1] += t * v2; 00230 /* L110: */ 00231 } 00232 00233 a[l1 + k * a_dim1] = 0.; 00234 00235 i__3 = *n; 00236 for (j = l; j <= i__3; ++j) { 00237 t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1]; 00238 b[l + j * b_dim1] += t * v1; 00239 b[l1 + j * b_dim1] += t * v2; 00240 /* L120: */ 00241 } 00242 /* .......... ZERO B(L+1,L) .......... */ 00243 s = (d__1 = b[l1 + l1 * b_dim1], abs(d__1)) + (d__2 = b[l1 + l * 00244 b_dim1], abs(d__2)); 00245 if (s == 0.) { 00246 goto L150; 00247 } 00248 u1 = b[l1 + l1 * b_dim1] / s; 00249 u2 = b[l1 + l * b_dim1] / s; 00250 d__1 = sqrt(u1 * u1 + u2 * u2); 00251 r__ = d_sign(&d__1, &u1); 00252 v1 = -(u1 + r__) / r__; 00253 v2 = -u2 / r__; 00254 u2 = v2 / v1; 00255 00256 i__3 = l1; 00257 for (i__ = 1; i__ <= i__3; ++i__) { 00258 t = b[i__ + l1 * b_dim1] + u2 * b[i__ + l * b_dim1]; 00259 b[i__ + l1 * b_dim1] += t * v1; 00260 b[i__ + l * b_dim1] += t * v2; 00261 /* L130: */ 00262 } 00263 00264 b[l1 + l * b_dim1] = 0.; 00265 00266 i__3 = *n; 00267 for (i__ = 1; i__ <= i__3; ++i__) { 00268 t = a[i__ + l1 * a_dim1] + u2 * a[i__ + l * a_dim1]; 00269 a[i__ + l1 * a_dim1] += t * v1; 00270 a[i__ + l * a_dim1] += t * v2; 00271 /* L140: */ 00272 } 00273 00274 if (! (*matz)) { 00275 goto L150; 00276 } 00277 00278 i__3 = *n; 00279 for (i__ = 1; i__ <= i__3; ++i__) { 00280 t = z__[i__ + l1 * z_dim1] + u2 * z__[i__ + l * z_dim1]; 00281 z__[i__ + l1 * z_dim1] += t * v1; 00282 z__[i__ + l * z_dim1] += t * v2; 00283 /* L145: */ 00284 } 00285 00286 L150: 00287 ; 00288 } 00289 00290 /* L160: */ 00291 } 00292 00293 L170: 00294 return 0; 00295 } /* qzhes_ */ |