Doxygen Source Code Documentation
eis_bandr.c File Reference
#include "f2c.h"Go to the source code of this file.
Functions | |
| int | bandr_ (integer *nm, integer *n, integer *mb, doublereal *a, doublereal *d__, doublereal *e, doublereal *e2, logical *matz, doublereal *z__) |
Function Documentation
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 8 of file eis_bandr.c. References a, i1, i2, l, m1, max, min, n2, and s2. Referenced by rsb_().
00011 {
00012 /* System generated locals */
00013 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5,
00014 i__6;
00015 doublereal d__1;
00016
00017 /* Builtin functions */
00018 double sqrt(doublereal);
00019
00020 /* Local variables */
00021 static doublereal dmin__;
00022 static integer maxl, maxr;
00023 static doublereal g;
00024 static integer j, k, l, r__;
00025 static doublereal u, b1, b2, c2, f1, f2;
00026 static integer i1, i2, j1, j2, m1, n2, r1;
00027 static doublereal s2;
00028 static integer kr, mr;
00029 static doublereal dminrt;
00030 static integer ugl;
00031
00032
00033
00034 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD, */
00035 /* NUM. MATH. 12, 231-241(1968) BY SCHWARZ. */
00036 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). */
00037
00038 /* THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX */
00039 /* TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY */
00040 /* ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
00041
00042 /* ON INPUT */
00043
00044 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00045 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00046 /* DIMENSION STATEMENT. */
00047
00048 /* N IS THE ORDER OF THE MATRIX. */
00049
00050 /* MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
00051 /* NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
00052 /* DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
00053 /* LOWER TRIANGLE OF THE MATRIX. */
00054
00055 /* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
00056 /* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */
00057 /* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
00058 /* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
00059 /* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
00060 /* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
00061 */
00062 /* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
00063
00064 /* MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS */
00065 /* TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE. */
00066
00067 /* ON OUTPUT */
00068
00069 /* A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH */
00070 /* CONTAIN A COPY OF THE TRIDIAGONAL MATRIX. */
00071
00072 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
00073
00074 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
00075 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */
00076
00077 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00078 /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
00079
00080 /* Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN */
00081 /* THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z */
00082 /* IS NOT REFERENCED. */
00083
00084 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00085 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
00086 */
00087
00088 /* THIS VERSION DATED AUGUST 1983. */
00089
00090 /* ------------------------------------------------------------------
00091 */
00092
00093 /* Parameter adjustments */
00094 z_dim1 = *nm;
00095 z_offset = z_dim1 + 1;
00096 z__ -= z_offset;
00097 --e2;
00098 --e;
00099 --d__;
00100 a_dim1 = *nm;
00101 a_offset = a_dim1 + 1;
00102 a -= a_offset;
00103
00104 /* Function Body */
00105 dmin__ = 5.4210108624275222e-20;
00106 dminrt = 2.3283064365386963e-10;
00107 /* .......... INITIALIZE DIAGONAL SCALING MATRIX .......... */
00108 i__1 = *n;
00109 for (j = 1; j <= i__1; ++j) {
00110 /* L30: */
00111 d__[j] = 1.;
00112 }
00113
00114 if (! (*matz)) {
00115 goto L60;
00116 }
00117
00118 i__1 = *n;
00119 for (j = 1; j <= i__1; ++j) {
00120
00121 i__2 = *n;
00122 for (k = 1; k <= i__2; ++k) {
00123 /* L40: */
00124 z__[j + k * z_dim1] = 0.;
00125 }
00126
00127 z__[j + j * z_dim1] = 1.;
00128 /* L50: */
00129 }
00130
00131 L60:
00132 m1 = *mb - 1;
00133 if ((i__1 = m1 - 1) < 0) {
00134 goto L900;
00135 } else if (i__1 == 0) {
00136 goto L800;
00137 } else {
00138 goto L70;
00139 }
00140 L70:
00141 n2 = *n - 2;
00142
00143 i__1 = n2;
00144 for (k = 1; k <= i__1; ++k) {
00145 /* Computing MIN */
00146 i__2 = m1, i__3 = *n - k;
00147 maxr = min(i__2,i__3);
00148 /* .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- .......... */
00149 i__2 = maxr;
00150 for (r1 = 2; r1 <= i__2; ++r1) {
00151 r__ = maxr + 2 - r1;
00152 kr = k + r__;
00153 mr = *mb - r__;
00154 g = a[kr + mr * a_dim1];
00155 a[kr - 1 + a_dim1] = a[kr - 1 + (mr + 1) * a_dim1];
00156 ugl = k;
00157
00158 i__3 = *n;
00159 i__4 = m1;
00160 for (j = kr; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
00161 j1 = j - 1;
00162 j2 = j1 - 1;
00163 if (g == 0.) {
00164 goto L600;
00165 }
00166 b1 = a[j1 + a_dim1] / g;
00167 b2 = b1 * d__[j1] / d__[j];
00168 s2 = 1. / (b1 * b2 + 1.);
00169 if (s2 >= .5) {
00170 goto L450;
00171 }
00172 b1 = g / a[j1 + a_dim1];
00173 b2 = b1 * d__[j] / d__[j1];
00174 c2 = 1. - s2;
00175 d__[j1] = c2 * d__[j1];
00176 d__[j] = c2 * d__[j];
00177 f1 = a[j + m1 * a_dim1] * 2.;
00178 f2 = b1 * a[j1 + *mb * a_dim1];
00179 a[j + m1 * a_dim1] = -b2 * (b1 * a[j + m1 * a_dim1] - a[j + *
00180 mb * a_dim1]) - f2 + a[j + m1 * a_dim1];
00181 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j + *mb * a_dim1] + f1) +
00182 a[j1 + *mb * a_dim1];
00183 a[j + *mb * a_dim1] = b1 * (f2 - f1) + a[j + *mb * a_dim1];
00184
00185 i__5 = j2;
00186 for (l = ugl; l <= i__5; ++l) {
00187 i2 = *mb - j + l;
00188 u = a[j1 + (i2 + 1) * a_dim1] + b2 * a[j + i2 * a_dim1];
00189 a[j + i2 * a_dim1] = -b1 * a[j1 + (i2 + 1) * a_dim1] + a[
00190 j + i2 * a_dim1];
00191 a[j1 + (i2 + 1) * a_dim1] = u;
00192 /* L200: */
00193 }
00194
00195 ugl = j;
00196 a[j1 + a_dim1] += b2 * g;
00197 if (j == *n) {
00198 goto L350;
00199 }
00200 /* Computing MIN */
00201 i__5 = m1, i__6 = *n - j1;
00202 maxl = min(i__5,i__6);
00203
00204 i__5 = maxl;
00205 for (l = 2; l <= i__5; ++l) {
00206 i1 = j1 + l;
00207 i2 = *mb - l;
00208 u = a[i1 + i2 * a_dim1] + b2 * a[i1 + (i2 + 1) * a_dim1];
00209 a[i1 + (i2 + 1) * a_dim1] = -b1 * a[i1 + i2 * a_dim1] + a[
00210 i1 + (i2 + 1) * a_dim1];
00211 a[i1 + i2 * a_dim1] = u;
00212 /* L300: */
00213 }
00214
00215 i1 = j + m1;
00216 if (i1 > *n) {
00217 goto L350;
00218 }
00219 g = b2 * a[i1 + a_dim1];
00220 L350:
00221 if (! (*matz)) {
00222 goto L500;
00223 }
00224
00225 i__5 = *n;
00226 for (l = 1; l <= i__5; ++l) {
00227 u = z__[l + j1 * z_dim1] + b2 * z__[l + j * z_dim1];
00228 z__[l + j * z_dim1] = -b1 * z__[l + j1 * z_dim1] + z__[l
00229 + j * z_dim1];
00230 z__[l + j1 * z_dim1] = u;
00231 /* L400: */
00232 }
00233
00234 goto L500;
00235
00236 L450:
00237 u = d__[j1];
00238 d__[j1] = s2 * d__[j];
00239 d__[j] = s2 * u;
00240 f1 = a[j + m1 * a_dim1] * 2.;
00241 f2 = b1 * a[j + *mb * a_dim1];
00242 u = b1 * (f2 - f1) + a[j1 + *mb * a_dim1];
00243 a[j + m1 * a_dim1] = b2 * (b1 * a[j + m1 * a_dim1] - a[j1 + *
00244 mb * a_dim1]) + f2 - a[j + m1 * a_dim1];
00245 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j1 + *mb * a_dim1] + f1)
00246 + a[j + *mb * a_dim1];
00247 a[j + *mb * a_dim1] = u;
00248
00249 i__5 = j2;
00250 for (l = ugl; l <= i__5; ++l) {
00251 i2 = *mb - j + l;
00252 u = b2 * a[j1 + (i2 + 1) * a_dim1] + a[j + i2 * a_dim1];
00253 a[j + i2 * a_dim1] = -a[j1 + (i2 + 1) * a_dim1] + b1 * a[
00254 j + i2 * a_dim1];
00255 a[j1 + (i2 + 1) * a_dim1] = u;
00256 /* L460: */
00257 }
00258
00259 ugl = j;
00260 a[j1 + a_dim1] = b2 * a[j1 + a_dim1] + g;
00261 if (j == *n) {
00262 goto L480;
00263 }
00264 /* Computing MIN */
00265 i__5 = m1, i__6 = *n - j1;
00266 maxl = min(i__5,i__6);
00267
00268 i__5 = maxl;
00269 for (l = 2; l <= i__5; ++l) {
00270 i1 = j1 + l;
00271 i2 = *mb - l;
00272 u = b2 * a[i1 + i2 * a_dim1] + a[i1 + (i2 + 1) * a_dim1];
00273 a[i1 + (i2 + 1) * a_dim1] = -a[i1 + i2 * a_dim1] + b1 * a[
00274 i1 + (i2 + 1) * a_dim1];
00275 a[i1 + i2 * a_dim1] = u;
00276 /* L470: */
00277 }
00278
00279 i1 = j + m1;
00280 if (i1 > *n) {
00281 goto L480;
00282 }
00283 g = a[i1 + a_dim1];
00284 a[i1 + a_dim1] = b1 * a[i1 + a_dim1];
00285 L480:
00286 if (! (*matz)) {
00287 goto L500;
00288 }
00289
00290 i__5 = *n;
00291 for (l = 1; l <= i__5; ++l) {
00292 u = b2 * z__[l + j1 * z_dim1] + z__[l + j * z_dim1];
00293 z__[l + j * z_dim1] = -z__[l + j1 * z_dim1] + b1 * z__[l
00294 + j * z_dim1];
00295 z__[l + j1 * z_dim1] = u;
00296 /* L490: */
00297 }
00298
00299 L500:
00300 ;
00301 }
00302
00303 L600:
00304 ;
00305 }
00306
00307 if (k % 64 != 0) {
00308 goto L700;
00309 }
00310 /* .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW .......... */
00311 i__2 = *n;
00312 for (j = k; j <= i__2; ++j) {
00313 if (d__[j] >= dmin__) {
00314 goto L650;
00315 }
00316 /* Computing MAX */
00317 i__4 = 1, i__3 = *mb + 1 - j;
00318 maxl = max(i__4,i__3);
00319
00320 i__4 = m1;
00321 for (l = maxl; l <= i__4; ++l) {
00322 /* L610: */
00323 a[j + l * a_dim1] = dminrt * a[j + l * a_dim1];
00324 }
00325
00326 if (j == *n) {
00327 goto L630;
00328 }
00329 /* Computing MIN */
00330 i__4 = m1, i__3 = *n - j;
00331 maxl = min(i__4,i__3);
00332
00333 i__4 = maxl;
00334 for (l = 1; l <= i__4; ++l) {
00335 i1 = j + l;
00336 i2 = *mb - l;
00337 a[i1 + i2 * a_dim1] = dminrt * a[i1 + i2 * a_dim1];
00338 /* L620: */
00339 }
00340
00341 L630:
00342 if (! (*matz)) {
00343 goto L645;
00344 }
00345
00346 i__4 = *n;
00347 for (l = 1; l <= i__4; ++l) {
00348 /* L640: */
00349 z__[l + j * z_dim1] = dminrt * z__[l + j * z_dim1];
00350 }
00351
00352 L645:
00353 a[j + *mb * a_dim1] = dmin__ * a[j + *mb * a_dim1];
00354 d__[j] /= dmin__;
00355 L650:
00356 ;
00357 }
00358
00359 L700:
00360 ;
00361 }
00362 /* .......... FORM SQUARE ROOT OF SCALING MATRIX .......... */
00363 L800:
00364 i__1 = *n;
00365 for (j = 2; j <= i__1; ++j) {
00366 /* L810: */
00367 e[j] = sqrt(d__[j]);
00368 }
00369
00370 if (! (*matz)) {
00371 goto L840;
00372 }
00373
00374 i__1 = *n;
00375 for (j = 1; j <= i__1; ++j) {
00376
00377 i__2 = *n;
00378 for (k = 2; k <= i__2; ++k) {
00379 /* L820: */
00380 z__[j + k * z_dim1] = e[k] * z__[j + k * z_dim1];
00381 }
00382
00383 /* L830: */
00384 }
00385
00386 L840:
00387 u = 1.;
00388
00389 i__1 = *n;
00390 for (j = 2; j <= i__1; ++j) {
00391 a[j + m1 * a_dim1] = u * e[j] * a[j + m1 * a_dim1];
00392 u = e[j];
00393 /* Computing 2nd power */
00394 d__1 = a[j + m1 * a_dim1];
00395 e2[j] = d__1 * d__1;
00396 a[j + *mb * a_dim1] = d__[j] * a[j + *mb * a_dim1];
00397 d__[j] = a[j + *mb * a_dim1];
00398 e[j] = a[j + m1 * a_dim1];
00399 /* L850: */
00400 }
00401
00402 d__[1] = a[*mb * a_dim1 + 1];
00403 e[1] = 0.;
00404 e2[1] = 0.;
00405 goto L1001;
00406
00407 L900:
00408 i__1 = *n;
00409 for (j = 1; j <= i__1; ++j) {
00410 d__[j] = a[j + *mb * a_dim1];
00411 e[j] = 0.;
00412 e2[j] = 0.;
00413 /* L950: */
00414 }
00415
00416 L1001:
00417 return 0;
00418 } /* bandr_ */
|
