Doxygen Source Code Documentation
eis_tql1.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | tql1_ (integer *n, doublereal *d__, doublereal *e, integer *ierr) |
Variables | |
doublereal | c_b10 = 1. |
Function Documentation
|
Definition at line 12 of file eis_tql1.c. References abs, c_b10, d_sign(), l, p, pythag_(), and s2.
00014 { 00015 /* System generated locals */ 00016 integer i__1, i__2; 00017 doublereal d__1, d__2; 00018 00019 /* Builtin functions */ 00020 double d_sign(doublereal *, doublereal *); 00021 00022 /* Local variables */ 00023 static doublereal c__, f, g, h__; 00024 static integer i__, j, l, m; 00025 static doublereal p, r__, s, c2, c3; 00026 static integer l1, l2; 00027 static doublereal s2; 00028 static integer ii; 00029 extern doublereal pythag_(doublereal *, doublereal *); 00030 static doublereal dl1, el1; 00031 static integer mml; 00032 static doublereal tst1, tst2; 00033 00034 00035 00036 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL1, */ 00037 /* NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */ 00038 /* WILKINSON. */ 00039 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */ 00040 00041 /* THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC */ 00042 /* TRIDIAGONAL MATRIX BY THE QL METHOD. */ 00043 00044 /* ON INPUT */ 00045 00046 /* N IS THE ORDER OF THE MATRIX. */ 00047 00048 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */ 00049 00050 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */ 00051 /* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */ 00052 00053 /* ON OUTPUT */ 00054 00055 /* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */ 00056 /* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */ 00057 /* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */ 00058 /* THE SMALLEST EIGENVALUES. */ 00059 00060 /* E HAS BEEN DESTROYED. */ 00061 00062 /* IERR IS SET TO */ 00063 /* ZERO FOR NORMAL RETURN, */ 00064 /* J IF THE J-TH EIGENVALUE HAS NOT BEEN */ 00065 /* DETERMINED AFTER 30 ITERATIONS. */ 00066 00067 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ 00068 00069 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00070 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00071 */ 00072 00073 /* THIS VERSION DATED AUGUST 1983. */ 00074 00075 /* ------------------------------------------------------------------ 00076 */ 00077 00078 /* Parameter adjustments */ 00079 --e; 00080 --d__; 00081 00082 /* Function Body */ 00083 *ierr = 0; 00084 if (*n == 1) { 00085 goto L1001; 00086 } 00087 00088 i__1 = *n; 00089 for (i__ = 2; i__ <= i__1; ++i__) { 00090 /* L100: */ 00091 e[i__ - 1] = e[i__]; 00092 } 00093 00094 f = 0.; 00095 tst1 = 0.; 00096 e[*n] = 0.; 00097 00098 i__1 = *n; 00099 for (l = 1; l <= i__1; ++l) { 00100 j = 0; 00101 h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2)); 00102 if (tst1 < h__) { 00103 tst1 = h__; 00104 } 00105 /* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */ 00106 i__2 = *n; 00107 for (m = l; m <= i__2; ++m) { 00108 tst2 = tst1 + (d__1 = e[m], abs(d__1)); 00109 if (tst2 == tst1) { 00110 goto L120; 00111 } 00112 /* .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT */ 00113 /* THROUGH THE BOTTOM OF THE LOOP .......... */ 00114 /* L110: */ 00115 } 00116 00117 L120: 00118 if (m == l) { 00119 goto L210; 00120 } 00121 L130: 00122 if (j == 30) { 00123 goto L1000; 00124 } 00125 ++j; 00126 /* .......... FORM SHIFT .......... */ 00127 l1 = l + 1; 00128 l2 = l1 + 1; 00129 g = d__[l]; 00130 p = (d__[l1] - g) / (e[l] * 2.); 00131 r__ = pythag_(&p, &c_b10); 00132 d__[l] = e[l] / (p + d_sign(&r__, &p)); 00133 d__[l1] = e[l] * (p + d_sign(&r__, &p)); 00134 dl1 = d__[l1]; 00135 h__ = g - d__[l]; 00136 if (l2 > *n) { 00137 goto L145; 00138 } 00139 00140 i__2 = *n; 00141 for (i__ = l2; i__ <= i__2; ++i__) { 00142 /* L140: */ 00143 d__[i__] -= h__; 00144 } 00145 00146 L145: 00147 f += h__; 00148 /* .......... QL TRANSFORMATION .......... */ 00149 p = d__[m]; 00150 c__ = 1.; 00151 c2 = c__; 00152 el1 = e[l1]; 00153 s = 0.; 00154 mml = m - l; 00155 /* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */ 00156 i__2 = mml; 00157 for (ii = 1; ii <= i__2; ++ii) { 00158 c3 = c2; 00159 c2 = c__; 00160 s2 = s; 00161 i__ = m - ii; 00162 g = c__ * e[i__]; 00163 h__ = c__ * p; 00164 r__ = pythag_(&p, &e[i__]); 00165 e[i__ + 1] = s * r__; 00166 s = e[i__] / r__; 00167 c__ = p / r__; 00168 p = c__ * d__[i__] - s * g; 00169 d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]); 00170 /* L200: */ 00171 } 00172 00173 p = -s * s2 * c3 * el1 * e[l] / dl1; 00174 e[l] = s * p; 00175 d__[l] = c__ * p; 00176 tst2 = tst1 + (d__1 = e[l], abs(d__1)); 00177 if (tst2 > tst1) { 00178 goto L130; 00179 } 00180 L210: 00181 p = d__[l] + f; 00182 /* .......... ORDER EIGENVALUES .......... */ 00183 if (l == 1) { 00184 goto L250; 00185 } 00186 /* .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */ 00187 i__2 = l; 00188 for (ii = 2; ii <= i__2; ++ii) { 00189 i__ = l + 2 - ii; 00190 if (p >= d__[i__ - 1]) { 00191 goto L270; 00192 } 00193 d__[i__] = d__[i__ - 1]; 00194 /* L230: */ 00195 } 00196 00197 L250: 00198 i__ = 1; 00199 L270: 00200 d__[i__] = p; 00201 /* L290: */ 00202 } 00203 00204 goto L1001; 00205 /* .......... SET ERROR -- NO CONVERGENCE TO AN */ 00206 /* EIGENVALUE AFTER 30 ITERATIONS .......... */ 00207 L1000: 00208 *ierr = l; 00209 L1001: 00210 return 0; 00211 } /* tql1_ */ |
Variable Documentation
|
Definition at line 10 of file eis_tql1.c. Referenced by tql1_(). |