Doxygen Source Code Documentation
eis_tinvit.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | tinvit_ (integer *nm, integer *n, doublereal *d__, doublereal *e, doublereal *e2, integer *m, doublereal *w, integer *ind, doublereal *z__, integer *ierr, doublereal *rv1, doublereal *rv2, doublereal *rv3, doublereal *rv4, doublereal *rv6) |
Function Documentation
|
Definition at line 8 of file eis_tinvit.c. References abs, epslon_(), ind, max, p, pythag_(), q, v, and x0. Referenced by rsm_().
00012 { 00013 /* System generated locals */ 00014 integer z_dim1, z_offset, i__1, i__2, i__3; 00015 doublereal d__1, d__2, d__3, d__4; 00016 00017 /* Builtin functions */ 00018 double sqrt(doublereal); 00019 00020 /* Local variables */ 00021 static doublereal norm; 00022 static integer i__, j, p, q, r__, s; 00023 static doublereal u, v, order; 00024 static integer group; 00025 static doublereal x0, x1; 00026 static integer ii, jj, ip; 00027 static doublereal uk, xu; 00028 extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal 00029 *); 00030 static integer tag, its; 00031 static doublereal eps2, eps3, eps4; 00032 00033 00034 00035 /* THIS SUBROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH- */ 00036 /* NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. */ 00037 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). */ 00038 00039 /* THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL */ 00040 /* SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, */ 00041 /* USING INVERSE ITERATION. */ 00042 00043 /* ON INPUT */ 00044 00045 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00046 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00047 /* DIMENSION STATEMENT. */ 00048 00049 /* N IS THE ORDER OF THE MATRIX. */ 00050 00051 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */ 00052 00053 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */ 00054 /* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */ 00055 00056 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E, */ 00057 /* WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E. */ 00058 /* E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN */ 00059 /* THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM */ 00060 /* OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN */ 00061 /* 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0D0 */ 00062 /* IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT, */ 00063 /* TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES, */ 00064 /* THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE. */ 00065 00066 /* M IS THE NUMBER OF SPECIFIED EIGENVALUES. */ 00067 00068 /* W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. 00069 */ 00070 00071 /* IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */ 00072 /* ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */ 00073 /* 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */ 00074 /* THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC. 00075 */ 00076 00077 /* ON OUTPUT */ 00078 00079 /* ALL INPUT ARRAYS ARE UNALTERED. */ 00080 00081 /* Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. */ 00082 /* ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. */ 00083 00084 /* IERR IS SET TO */ 00085 /* ZERO FOR NORMAL RETURN, */ 00086 /* -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */ 00087 /* EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS. */ 00088 00089 /* RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS. */ 00090 00091 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ 00092 00093 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00094 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00095 */ 00096 00097 /* THIS VERSION DATED AUGUST 1983. */ 00098 00099 /* ------------------------------------------------------------------ 00100 */ 00101 00102 /* Parameter adjustments */ 00103 --rv6; 00104 --rv4; 00105 --rv3; 00106 --rv2; 00107 --rv1; 00108 --e2; 00109 --e; 00110 --d__; 00111 z_dim1 = *nm; 00112 z_offset = z_dim1 + 1; 00113 z__ -= z_offset; 00114 --ind; 00115 --w; 00116 00117 /* Function Body */ 00118 *ierr = 0; 00119 if (*m == 0) { 00120 goto L1001; 00121 } 00122 tag = 0; 00123 order = 1. - e2[1]; 00124 q = 0; 00125 /* .......... ESTABLISH AND PROCESS NEXT SUBMATRIX .......... */ 00126 L100: 00127 p = q + 1; 00128 00129 i__1 = *n; 00130 for (q = p; q <= i__1; ++q) { 00131 if (q == *n) { 00132 goto L140; 00133 } 00134 if (e2[q + 1] == 0.) { 00135 goto L140; 00136 } 00137 /* L120: */ 00138 } 00139 /* .......... FIND VECTORS BY INVERSE ITERATION .......... */ 00140 L140: 00141 ++tag; 00142 s = 0; 00143 00144 i__1 = *m; 00145 for (r__ = 1; r__ <= i__1; ++r__) { 00146 if (ind[r__] != tag) { 00147 goto L920; 00148 } 00149 its = 1; 00150 x1 = w[r__]; 00151 if (s != 0) { 00152 goto L510; 00153 } 00154 /* .......... CHECK FOR ISOLATED ROOT .......... */ 00155 xu = 1.; 00156 if (p != q) { 00157 goto L490; 00158 } 00159 rv6[p] = 1.; 00160 goto L870; 00161 L490: 00162 norm = (d__1 = d__[p], abs(d__1)); 00163 ip = p + 1; 00164 00165 i__2 = q; 00166 for (i__ = ip; i__ <= i__2; ++i__) { 00167 /* L500: */ 00168 /* Computing MAX */ 00169 d__3 = norm, d__4 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[i__], 00170 abs(d__2)); 00171 norm = max(d__3,d__4); 00172 } 00173 /* .......... EPS2 IS THE CRITERION FOR GROUPING, */ 00174 /* EPS3 REPLACES ZERO PIVOTS AND EQUAL */ 00175 /* ROOTS ARE MODIFIED BY EPS3, */ 00176 /* EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ......... 00177 . */ 00178 eps2 = norm * .001; 00179 eps3 = epslon_(&norm); 00180 uk = (doublereal) (q - p + 1); 00181 eps4 = uk * eps3; 00182 uk = eps4 / sqrt(uk); 00183 s = p; 00184 L505: 00185 group = 0; 00186 goto L520; 00187 /* .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */ 00188 L510: 00189 if ((d__1 = x1 - x0, abs(d__1)) >= eps2) { 00190 goto L505; 00191 } 00192 ++group; 00193 if (order * (x1 - x0) <= 0.) { 00194 x1 = x0 + order * eps3; 00195 } 00196 /* .......... ELIMINATION WITH INTERCHANGES AND */ 00197 /* INITIALIZATION OF VECTOR .......... */ 00198 L520: 00199 v = 0.; 00200 00201 i__2 = q; 00202 for (i__ = p; i__ <= i__2; ++i__) { 00203 rv6[i__] = uk; 00204 if (i__ == p) { 00205 goto L560; 00206 } 00207 if ((d__1 = e[i__], abs(d__1)) < abs(u)) { 00208 goto L540; 00209 } 00210 /* .......... WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF */ 00211 /* E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY ...... 00212 .... */ 00213 xu = u / e[i__]; 00214 rv4[i__] = xu; 00215 rv1[i__ - 1] = e[i__]; 00216 rv2[i__ - 1] = d__[i__] - x1; 00217 rv3[i__ - 1] = 0.; 00218 if (i__ != q) { 00219 rv3[i__ - 1] = e[i__ + 1]; 00220 } 00221 u = v - xu * rv2[i__ - 1]; 00222 v = -xu * rv3[i__ - 1]; 00223 goto L580; 00224 L540: 00225 xu = e[i__] / u; 00226 rv4[i__] = xu; 00227 rv1[i__ - 1] = u; 00228 rv2[i__ - 1] = v; 00229 rv3[i__ - 1] = 0.; 00230 L560: 00231 u = d__[i__] - x1 - xu * v; 00232 if (i__ != q) { 00233 v = e[i__ + 1]; 00234 } 00235 L580: 00236 ; 00237 } 00238 00239 if (u == 0.) { 00240 u = eps3; 00241 } 00242 rv1[q] = u; 00243 rv2[q] = 0.; 00244 rv3[q] = 0.; 00245 /* .......... BACK SUBSTITUTION */ 00246 /* FOR I=Q STEP -1 UNTIL P DO -- .......... */ 00247 L600: 00248 i__2 = q; 00249 for (ii = p; ii <= i__2; ++ii) { 00250 i__ = p + q - ii; 00251 rv6[i__] = (rv6[i__] - u * rv2[i__] - v * rv3[i__]) / rv1[i__]; 00252 v = u; 00253 u = rv6[i__]; 00254 /* L620: */ 00255 } 00256 /* .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */ 00257 /* MEMBERS OF GROUP .......... */ 00258 if (group == 0) { 00259 goto L700; 00260 } 00261 j = r__; 00262 00263 i__2 = group; 00264 for (jj = 1; jj <= i__2; ++jj) { 00265 L630: 00266 --j; 00267 if (ind[j] != tag) { 00268 goto L630; 00269 } 00270 xu = 0.; 00271 00272 i__3 = q; 00273 for (i__ = p; i__ <= i__3; ++i__) { 00274 /* L640: */ 00275 xu += rv6[i__] * z__[i__ + j * z_dim1]; 00276 } 00277 00278 i__3 = q; 00279 for (i__ = p; i__ <= i__3; ++i__) { 00280 /* L660: */ 00281 rv6[i__] -= xu * z__[i__ + j * z_dim1]; 00282 } 00283 00284 /* L680: */ 00285 } 00286 00287 L700: 00288 norm = 0.; 00289 00290 i__2 = q; 00291 for (i__ = p; i__ <= i__2; ++i__) { 00292 /* L720: */ 00293 norm += (d__1 = rv6[i__], abs(d__1)); 00294 } 00295 00296 if (norm >= 1.) { 00297 goto L840; 00298 } 00299 /* .......... FORWARD SUBSTITUTION .......... */ 00300 if (its == 5) { 00301 goto L830; 00302 } 00303 if (norm != 0.) { 00304 goto L740; 00305 } 00306 rv6[s] = eps4; 00307 ++s; 00308 if (s > q) { 00309 s = p; 00310 } 00311 goto L780; 00312 L740: 00313 xu = eps4 / norm; 00314 00315 i__2 = q; 00316 for (i__ = p; i__ <= i__2; ++i__) { 00317 /* L760: */ 00318 rv6[i__] *= xu; 00319 } 00320 /* .......... ELIMINATION OPERATIONS ON NEXT VECTOR */ 00321 /* ITERATE .......... */ 00322 L780: 00323 i__2 = q; 00324 for (i__ = ip; i__ <= i__2; ++i__) { 00325 u = rv6[i__]; 00326 /* .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE */ 00327 /* WAS PERFORMED EARLIER IN THE */ 00328 /* TRIANGULARIZATION PROCESS .......... */ 00329 if (rv1[i__ - 1] != e[i__]) { 00330 goto L800; 00331 } 00332 u = rv6[i__ - 1]; 00333 rv6[i__ - 1] = rv6[i__]; 00334 L800: 00335 rv6[i__] = u - rv4[i__] * rv6[i__ - 1]; 00336 /* L820: */ 00337 } 00338 00339 ++its; 00340 goto L600; 00341 /* .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */ 00342 L830: 00343 *ierr = -r__; 00344 xu = 0.; 00345 goto L870; 00346 /* .......... NORMALIZE SO THAT SUM OF SQUARES IS */ 00347 /* 1 AND EXPAND TO FULL ORDER .......... */ 00348 L840: 00349 u = 0.; 00350 00351 i__2 = q; 00352 for (i__ = p; i__ <= i__2; ++i__) { 00353 /* L860: */ 00354 u = pythag_(&u, &rv6[i__]); 00355 } 00356 00357 xu = 1. / u; 00358 00359 L870: 00360 i__2 = *n; 00361 for (i__ = 1; i__ <= i__2; ++i__) { 00362 /* L880: */ 00363 z__[i__ + r__ * z_dim1] = 0.; 00364 } 00365 00366 i__2 = q; 00367 for (i__ = p; i__ <= i__2; ++i__) { 00368 /* L900: */ 00369 z__[i__ + r__ * z_dim1] = rv6[i__] * xu; 00370 } 00371 00372 x0 = x1; 00373 L920: 00374 ; 00375 } 00376 00377 if (q < *n) { 00378 goto L100; 00379 } 00380 L1001: 00381 return 0; 00382 } /* tinvit_ */ |