Doxygen Source Code Documentation
eis_htrid3.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | htrid3_ (integer *nm, integer *n, doublereal *a, doublereal *d__, doublereal *e, doublereal *e2, doublereal *tau) |
Function Documentation
|
Definition at line 8 of file eis_htrid3.c. References a, abs, l, pythag_(), and scale.
00010 { 00011 /* System generated locals */ 00012 integer a_dim1, a_offset, i__1, i__2, i__3; 00013 doublereal d__1, d__2; 00014 00015 /* Builtin functions */ 00016 double sqrt(doublereal); 00017 00018 /* Local variables */ 00019 static doublereal f, g, h__; 00020 static integer i__, j, k, l; 00021 static doublereal scale, fi, gi, hh; 00022 static integer ii; 00023 static doublereal si; 00024 extern doublereal pythag_(doublereal *, doublereal *); 00025 static integer jm1, jp1; 00026 00027 00028 00029 /* THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */ 00030 /* THE ALGOL PROCEDURE TRED3, NUM. MATH. 11, 181-195(1968) */ 00031 /* BY MARTIN, REINSCH, AND WILKINSON. */ 00032 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */ 00033 00034 /* THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX, STORED AS */ 00035 /* A SINGLE SQUARE ARRAY, TO A REAL SYMMETRIC TRIDIAGONAL MATRIX */ 00036 /* USING UNITARY SIMILARITY TRANSFORMATIONS. */ 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 MATRIX. */ 00045 00046 /* A CONTAINS THE LOWER TRIANGLE OF THE COMPLEX HERMITIAN INPUT */ 00047 /* MATRIX. THE REAL PARTS OF THE MATRIX ELEMENTS ARE STORED */ 00048 /* IN THE FULL LOWER TRIANGLE OF A, AND THE IMAGINARY PARTS */ 00049 /* ARE STORED IN THE TRANSPOSED POSITIONS OF THE STRICT UPPER */ 00050 /* TRIANGLE OF A. NO STORAGE IS REQUIRED FOR THE ZERO */ 00051 /* IMAGINARY PARTS OF THE DIAGONAL ELEMENTS. */ 00052 00053 /* ON OUTPUT */ 00054 00055 /* A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS */ 00056 /* USED IN THE REDUCTION. */ 00057 00058 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX. 00059 */ 00060 00061 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */ 00062 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */ 00063 00064 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */ 00065 /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */ 00066 00067 /* TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. */ 00068 00069 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ 00070 00071 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00072 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00073 */ 00074 00075 /* THIS VERSION DATED AUGUST 1983. */ 00076 00077 /* ------------------------------------------------------------------ 00078 */ 00079 00080 /* Parameter adjustments */ 00081 tau -= 3; 00082 --e2; 00083 --e; 00084 --d__; 00085 a_dim1 = *nm; 00086 a_offset = a_dim1 + 1; 00087 a -= a_offset; 00088 00089 /* Function Body */ 00090 tau[(*n << 1) + 1] = 1.; 00091 tau[(*n << 1) + 2] = 0.; 00092 /* .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */ 00093 i__1 = *n; 00094 for (ii = 1; ii <= i__1; ++ii) { 00095 i__ = *n + 1 - ii; 00096 l = i__ - 1; 00097 h__ = 0.; 00098 scale = 0.; 00099 if (l < 1) { 00100 goto L130; 00101 } 00102 /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */ 00103 i__2 = l; 00104 for (k = 1; k <= i__2; ++k) { 00105 /* L120: */ 00106 scale = scale + (d__1 = a[i__ + k * a_dim1], abs(d__1)) + (d__2 = 00107 a[k + i__ * a_dim1], abs(d__2)); 00108 } 00109 00110 if (scale != 0.) { 00111 goto L140; 00112 } 00113 tau[(l << 1) + 1] = 1.; 00114 tau[(l << 1) + 2] = 0.; 00115 L130: 00116 e[i__] = 0.; 00117 e2[i__] = 0.; 00118 goto L290; 00119 00120 L140: 00121 i__2 = l; 00122 for (k = 1; k <= i__2; ++k) { 00123 a[i__ + k * a_dim1] /= scale; 00124 a[k + i__ * a_dim1] /= scale; 00125 h__ = h__ + a[i__ + k * a_dim1] * a[i__ + k * a_dim1] + a[k + i__ 00126 * a_dim1] * a[k + i__ * a_dim1]; 00127 /* L150: */ 00128 } 00129 00130 e2[i__] = scale * scale * h__; 00131 g = sqrt(h__); 00132 e[i__] = scale * g; 00133 f = pythag_(&a[i__ + l * a_dim1], &a[l + i__ * a_dim1]); 00134 /* .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T .......... */ 00135 if (f == 0.) { 00136 goto L160; 00137 } 00138 tau[(l << 1) + 1] = (a[l + i__ * a_dim1] * tau[(i__ << 1) + 2] - a[ 00139 i__ + l * a_dim1] * tau[(i__ << 1) + 1]) / f; 00140 si = (a[i__ + l * a_dim1] * tau[(i__ << 1) + 2] + a[l + i__ * a_dim1] 00141 * tau[(i__ << 1) + 1]) / f; 00142 h__ += f * g; 00143 g = g / f + 1.; 00144 a[i__ + l * a_dim1] = g * a[i__ + l * a_dim1]; 00145 a[l + i__ * a_dim1] = g * a[l + i__ * a_dim1]; 00146 if (l == 1) { 00147 goto L270; 00148 } 00149 goto L170; 00150 L160: 00151 tau[(l << 1) + 1] = -tau[(i__ << 1) + 1]; 00152 si = tau[(i__ << 1) + 2]; 00153 a[i__ + l * a_dim1] = g; 00154 L170: 00155 f = 0.; 00156 00157 i__2 = l; 00158 for (j = 1; j <= i__2; ++j) { 00159 g = 0.; 00160 gi = 0.; 00161 if (j == 1) { 00162 goto L190; 00163 } 00164 jm1 = j - 1; 00165 /* .......... FORM ELEMENT OF A*U .......... */ 00166 i__3 = jm1; 00167 for (k = 1; k <= i__3; ++k) { 00168 g = g + a[j + k * a_dim1] * a[i__ + k * a_dim1] + a[k + j * 00169 a_dim1] * a[k + i__ * a_dim1]; 00170 gi = gi - a[j + k * a_dim1] * a[k + i__ * a_dim1] + a[k + j * 00171 a_dim1] * a[i__ + k * a_dim1]; 00172 /* L180: */ 00173 } 00174 00175 L190: 00176 g += a[j + j * a_dim1] * a[i__ + j * a_dim1]; 00177 gi -= a[j + j * a_dim1] * a[j + i__ * a_dim1]; 00178 jp1 = j + 1; 00179 if (l < jp1) { 00180 goto L220; 00181 } 00182 00183 i__3 = l; 00184 for (k = jp1; k <= i__3; ++k) { 00185 g = g + a[k + j * a_dim1] * a[i__ + k * a_dim1] - a[j + k * 00186 a_dim1] * a[k + i__ * a_dim1]; 00187 gi = gi - a[k + j * a_dim1] * a[k + i__ * a_dim1] - a[j + k * 00188 a_dim1] * a[i__ + k * a_dim1]; 00189 /* L200: */ 00190 } 00191 /* .......... FORM ELEMENT OF P .......... */ 00192 L220: 00193 e[j] = g / h__; 00194 tau[(j << 1) + 2] = gi / h__; 00195 f = f + e[j] * a[i__ + j * a_dim1] - tau[(j << 1) + 2] * a[j + 00196 i__ * a_dim1]; 00197 /* L240: */ 00198 } 00199 00200 hh = f / (h__ + h__); 00201 /* .......... FORM REDUCED A .......... */ 00202 i__2 = l; 00203 for (j = 1; j <= i__2; ++j) { 00204 f = a[i__ + j * a_dim1]; 00205 g = e[j] - hh * f; 00206 e[j] = g; 00207 fi = -a[j + i__ * a_dim1]; 00208 gi = tau[(j << 1) + 2] - hh * fi; 00209 tau[(j << 1) + 2] = -gi; 00210 a[j + j * a_dim1] -= (f * g + fi * gi) * 2.; 00211 if (j == 1) { 00212 goto L260; 00213 } 00214 jm1 = j - 1; 00215 00216 i__3 = jm1; 00217 for (k = 1; k <= i__3; ++k) { 00218 a[j + k * a_dim1] = a[j + k * a_dim1] - f * e[k] - g * a[i__ 00219 + k * a_dim1] + fi * tau[(k << 1) + 2] + gi * a[k + 00220 i__ * a_dim1]; 00221 a[k + j * a_dim1] = a[k + j * a_dim1] - f * tau[(k << 1) + 2] 00222 - g * a[k + i__ * a_dim1] - fi * e[k] - gi * a[i__ + 00223 k * a_dim1]; 00224 /* L250: */ 00225 } 00226 00227 L260: 00228 ; 00229 } 00230 00231 L270: 00232 i__2 = l; 00233 for (k = 1; k <= i__2; ++k) { 00234 a[i__ + k * a_dim1] = scale * a[i__ + k * a_dim1]; 00235 a[k + i__ * a_dim1] = scale * a[k + i__ * a_dim1]; 00236 /* L280: */ 00237 } 00238 00239 tau[(l << 1) + 2] = -si; 00240 L290: 00241 d__[i__] = a[i__ + i__ * a_dim1]; 00242 a[i__ + i__ * a_dim1] = scale * sqrt(h__); 00243 /* L300: */ 00244 } 00245 00246 return 0; 00247 } /* htrid3_ */ |