Doxygen Source Code Documentation
eis_tred1.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | tred1_ (integer *nm, integer *n, doublereal *a, doublereal *d__, doublereal *e, doublereal *e2) |
Function Documentation
|
Definition at line 8 of file eis_tred1.c. References a, abs, d_sign(), l, and scale. Referenced by rs_(), rsg_(), rsgab_(), rsgba_(), and rsm_().
00010 { 00011 /* System generated locals */ 00012 integer a_dim1, a_offset, i__1, i__2, i__3; 00013 doublereal d__1; 00014 00015 /* Builtin functions */ 00016 double 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; 00022 static integer ii, jp1; 00023 00024 00025 00026 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, */ 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 */ 00031 /* TO A SYMMETRIC TRIDIAGONAL MATRIX USING */ 00032 /* ORTHOGONAL SIMILARITY TRANSFORMATIONS. */ 00033 00034 /* ON INPUT */ 00035 00036 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00037 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00038 /* DIMENSION STATEMENT. */ 00039 00040 /* N IS THE ORDER OF THE MATRIX. */ 00041 00042 /* A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE */ 00043 /* LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */ 00044 00045 /* ON OUTPUT */ 00046 00047 /* A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- */ 00048 /* FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER */ 00049 /* TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED. */ 00050 00051 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */ 00052 00053 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */ 00054 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */ 00055 00056 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */ 00057 /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */ 00058 00059 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00060 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00061 */ 00062 00063 /* THIS VERSION DATED AUGUST 1983. */ 00064 00065 /* ------------------------------------------------------------------ 00066 */ 00067 00068 /* Parameter adjustments */ 00069 --e2; 00070 --e; 00071 --d__; 00072 a_dim1 = *nm; 00073 a_offset = a_dim1 + 1; 00074 a -= a_offset; 00075 00076 /* Function Body */ 00077 i__1 = *n; 00078 for (i__ = 1; i__ <= i__1; ++i__) { 00079 d__[i__] = a[*n + i__ * a_dim1]; 00080 a[*n + i__ * a_dim1] = a[i__ + i__ * a_dim1]; 00081 /* L100: */ 00082 } 00083 /* .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */ 00084 i__1 = *n; 00085 for (ii = 1; ii <= i__1; ++ii) { 00086 i__ = *n + 1 - ii; 00087 l = i__ - 1; 00088 h__ = 0.; 00089 scale = 0.; 00090 if (l < 1) { 00091 goto L130; 00092 } 00093 /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */ 00094 i__2 = l; 00095 for (k = 1; k <= i__2; ++k) { 00096 /* L120: */ 00097 scale += (d__1 = d__[k], abs(d__1)); 00098 } 00099 00100 if (scale != 0.) { 00101 goto L140; 00102 } 00103 00104 i__2 = l; 00105 for (j = 1; j <= i__2; ++j) { 00106 d__[j] = a[l + j * a_dim1]; 00107 a[l + j * a_dim1] = a[i__ + j * a_dim1]; 00108 a[i__ + j * a_dim1] = 0.; 00109 /* L125: */ 00110 } 00111 00112 L130: 00113 e[i__] = 0.; 00114 e2[i__] = 0.; 00115 goto L300; 00116 00117 L140: 00118 i__2 = l; 00119 for (k = 1; k <= i__2; ++k) { 00120 d__[k] /= scale; 00121 h__ += d__[k] * d__[k]; 00122 /* L150: */ 00123 } 00124 00125 e2[i__] = scale * scale * h__; 00126 f = d__[l]; 00127 d__1 = sqrt(h__); 00128 g = -d_sign(&d__1, &f); 00129 e[i__] = scale * g; 00130 h__ -= f * g; 00131 d__[l] = f - g; 00132 if (l == 1) { 00133 goto L285; 00134 } 00135 /* .......... FORM A*U .......... */ 00136 i__2 = l; 00137 for (j = 1; j <= i__2; ++j) { 00138 /* L170: */ 00139 e[j] = 0.; 00140 } 00141 00142 i__2 = l; 00143 for (j = 1; j <= i__2; ++j) { 00144 f = d__[j]; 00145 g = e[j] + a[j + j * a_dim1] * f; 00146 jp1 = j + 1; 00147 if (l < jp1) { 00148 goto L220; 00149 } 00150 00151 i__3 = l; 00152 for (k = jp1; k <= i__3; ++k) { 00153 g += a[k + j * a_dim1] * d__[k]; 00154 e[k] += a[k + j * a_dim1] * f; 00155 /* L200: */ 00156 } 00157 00158 L220: 00159 e[j] = g; 00160 /* L240: */ 00161 } 00162 /* .......... FORM P .......... */ 00163 f = 0.; 00164 00165 i__2 = l; 00166 for (j = 1; j <= i__2; ++j) { 00167 e[j] /= h__; 00168 f += e[j] * d__[j]; 00169 /* L245: */ 00170 } 00171 00172 h__ = f / (h__ + h__); 00173 /* .......... FORM Q .......... */ 00174 i__2 = l; 00175 for (j = 1; j <= i__2; ++j) { 00176 /* L250: */ 00177 e[j] -= h__ * d__[j]; 00178 } 00179 /* .......... FORM REDUCED A .......... */ 00180 i__2 = l; 00181 for (j = 1; j <= i__2; ++j) { 00182 f = d__[j]; 00183 g = e[j]; 00184 00185 i__3 = l; 00186 for (k = j; k <= i__3; ++k) { 00187 /* L260: */ 00188 a[k + j * a_dim1] = a[k + j * a_dim1] - f * e[k] - g * d__[k]; 00189 } 00190 00191 /* L280: */ 00192 } 00193 00194 L285: 00195 i__2 = l; 00196 for (j = 1; j <= i__2; ++j) { 00197 f = d__[j]; 00198 d__[j] = a[l + j * a_dim1]; 00199 a[l + j * a_dim1] = a[i__ + j * a_dim1]; 00200 a[i__ + j * a_dim1] = f * scale; 00201 /* L290: */ 00202 } 00203 00204 L300: 00205 ; 00206 } 00207 00208 return 0; 00209 } /* tred1_ */ |