Doxygen Source Code Documentation
eis_tred3.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | tred3_ (integer *n, integer *nv, doublereal *a, doublereal *d__, doublereal *e, doublereal *e2) |
Function Documentation
|
Definition at line 8 of file eis_tred3.c. References a, abs, d_sign(), l, and scale. Referenced by rsp_().
00010 { 00011 /* System generated locals */ 00012 integer i__1, i__2, i__3; 00013 doublereal d__1; 00014 00015 /* Builtin functions */ 00016 double sqrt(doublereal), d_sign(doublereal *, doublereal *); 00017 00018 /* Local variables */ 00019 static doublereal f, g, h__; 00020 static integer i__, j, k, l; 00021 static doublereal scale, hh; 00022 static integer ii, jk, iz, jm1; 00023 00024 00025 00026 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3, */ 00027 /* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. */ 00028 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */ 00029 00030 /* THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS */ 00031 /* A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX */ 00032 /* USING ORTHOGONAL SIMILARITY TRANSFORMATIONS. */ 00033 00034 /* ON INPUT */ 00035 00036 /* N IS THE ORDER OF THE MATRIX. */ 00037 00038 /* NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A */ 00039 /* AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */ 00040 00041 /* A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC */ 00042 /* INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL */ 00043 /* ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS. */ 00044 00045 /* ON OUTPUT */ 00046 00047 /* A CONTAINS INFORMATION ABOUT THE ORTHOGONAL */ 00048 /* TRANSFORMATIONS USED IN THE REDUCTION. */ 00049 00050 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */ 00051 00052 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */ 00053 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */ 00054 00055 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */ 00056 /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */ 00057 00058 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00059 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00060 */ 00061 00062 /* THIS VERSION DATED AUGUST 1983. */ 00063 00064 /* ------------------------------------------------------------------ 00065 */ 00066 00067 /* .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */ 00068 /* Parameter adjustments */ 00069 --e2; 00070 --e; 00071 --d__; 00072 --a; 00073 00074 /* Function Body */ 00075 i__1 = *n; 00076 for (ii = 1; ii <= i__1; ++ii) { 00077 i__ = *n + 1 - ii; 00078 l = i__ - 1; 00079 iz = i__ * l / 2; 00080 h__ = 0.; 00081 scale = 0.; 00082 if (l < 1) { 00083 goto L130; 00084 } 00085 /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */ 00086 i__2 = l; 00087 for (k = 1; k <= i__2; ++k) { 00088 ++iz; 00089 d__[k] = a[iz]; 00090 scale += (d__1 = d__[k], abs(d__1)); 00091 /* L120: */ 00092 } 00093 00094 if (scale != 0.) { 00095 goto L140; 00096 } 00097 L130: 00098 e[i__] = 0.; 00099 e2[i__] = 0.; 00100 goto L290; 00101 00102 L140: 00103 i__2 = l; 00104 for (k = 1; k <= i__2; ++k) { 00105 d__[k] /= scale; 00106 h__ += d__[k] * d__[k]; 00107 /* L150: */ 00108 } 00109 00110 e2[i__] = scale * scale * h__; 00111 f = d__[l]; 00112 d__1 = sqrt(h__); 00113 g = -d_sign(&d__1, &f); 00114 e[i__] = scale * g; 00115 h__ -= f * g; 00116 d__[l] = f - g; 00117 a[iz] = scale * d__[l]; 00118 if (l == 1) { 00119 goto L290; 00120 } 00121 jk = 1; 00122 00123 i__2 = l; 00124 for (j = 1; j <= i__2; ++j) { 00125 f = d__[j]; 00126 g = 0.; 00127 jm1 = j - 1; 00128 if (jm1 < 1) { 00129 goto L220; 00130 } 00131 00132 i__3 = jm1; 00133 for (k = 1; k <= i__3; ++k) { 00134 g += a[jk] * d__[k]; 00135 e[k] += a[jk] * f; 00136 ++jk; 00137 /* L200: */ 00138 } 00139 00140 L220: 00141 e[j] = g + a[jk] * f; 00142 ++jk; 00143 /* L240: */ 00144 } 00145 /* .......... FORM P .......... */ 00146 f = 0.; 00147 00148 i__2 = l; 00149 for (j = 1; j <= i__2; ++j) { 00150 e[j] /= h__; 00151 f += e[j] * d__[j]; 00152 /* L245: */ 00153 } 00154 00155 hh = f / (h__ + h__); 00156 /* .......... FORM Q .......... */ 00157 i__2 = l; 00158 for (j = 1; j <= i__2; ++j) { 00159 /* L250: */ 00160 e[j] -= hh * d__[j]; 00161 } 00162 00163 jk = 1; 00164 /* .......... FORM REDUCED A .......... */ 00165 i__2 = l; 00166 for (j = 1; j <= i__2; ++j) { 00167 f = d__[j]; 00168 g = e[j]; 00169 00170 i__3 = j; 00171 for (k = 1; k <= i__3; ++k) { 00172 a[jk] = a[jk] - f * e[k] - g * d__[k]; 00173 ++jk; 00174 /* L260: */ 00175 } 00176 00177 /* L280: */ 00178 } 00179 00180 L290: 00181 d__[i__] = a[iz + 1]; 00182 a[iz + 1] = scale * sqrt(h__); 00183 /* L300: */ 00184 } 00185 00186 return 0; 00187 } /* tred3_ */ |