1 /* tred1.f -- translated by f2c (version 19961017).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
tred1_(integer * nm,integer * n,doublereal * a,doublereal * d__,doublereal * e,doublereal * e2)8 /* Subroutine */ int tred1_(integer *nm, integer *n, doublereal *a,
9 doublereal *d__, doublereal *e, doublereal *e2)
10 {
11 /* System generated locals */
12 integer a_dim1, a_offset, i__1, i__2, i__3;
13 doublereal d__1;
14
15 /* Builtin functions */
16 double d_sign(doublereal *, doublereal *);
17
18 /* Local variables */
19 doublereal f, g, h__;
20 integer i__, j, k, l;
21 doublereal scale;
22 integer ii, jp1;
23
24
25
26 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, */
27 /* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. */
28 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */
29
30 /* THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX */
31 /* TO A SYMMETRIC TRIDIAGONAL MATRIX USING */
32 /* ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
33
34 /* ON INPUT */
35
36 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
37 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
38 /* DIMENSION STATEMENT. */
39
40 /* N IS THE ORDER OF THE MATRIX. */
41
42 /* A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE */
43 /* LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */
44
45 /* ON OUTPUT */
46
47 /* A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- */
48 /* FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER */
49 /* TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED. */
50
51 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
52
53 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
54 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */
55
56 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
57 /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
58
59 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
60 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
61 */
62
63 /* THIS VERSION DATED AUGUST 1983. */
64
65 /* ------------------------------------------------------------------
66 */
67
68 /* Parameter adjustments */
69 --e2;
70 --e;
71 --d__;
72 a_dim1 = *nm;
73 a_offset = a_dim1 + 1;
74 a -= a_offset;
75
76 /* Function Body */
77 i__1 = *n;
78 for (i__ = 1; i__ <= i__1; ++i__) {
79 d__[i__] = a[*n + i__ * a_dim1];
80 a[*n + i__ * a_dim1] = a[i__ + i__ * a_dim1];
81 /* L100: */
82 }
83 /* .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */
84 i__1 = *n;
85 for (ii = 1; ii <= i__1; ++ii) {
86 i__ = *n + 1 - ii;
87 l = i__ - 1;
88 h__ = 0.;
89 scale = 0.;
90 if (l < 1) {
91 goto L130;
92 }
93 /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
94 i__2 = l;
95 for (k = 1; k <= i__2; ++k) {
96 /* L120: */
97 scale += (d__1 = d__[k], abs(d__1));
98 }
99
100 if (scale != 0.) {
101 goto L140;
102 }
103
104 i__2 = l;
105 for (j = 1; j <= i__2; ++j) {
106 d__[j] = a[l + j * a_dim1];
107 a[l + j * a_dim1] = a[i__ + j * a_dim1];
108 a[i__ + j * a_dim1] = 0.;
109 /* L125: */
110 }
111
112 L130:
113 e[i__] = 0.;
114 e2[i__] = 0.;
115 goto L300;
116
117 L140:
118 i__2 = l;
119 for (k = 1; k <= i__2; ++k) {
120 d__[k] /= scale;
121 h__ += d__[k] * d__[k];
122 /* L150: */
123 }
124
125 e2[i__] = scale * scale * h__;
126 f = d__[l];
127 d__1 = sqrt(h__);
128 g = -d_sign(&d__1, &f);
129 e[i__] = scale * g;
130 h__ -= f * g;
131 d__[l] = f - g;
132 if (l == 1) {
133 goto L285;
134 }
135 /* .......... FORM A*U .......... */
136 i__2 = l;
137 for (j = 1; j <= i__2; ++j) {
138 /* L170: */
139 e[j] = 0.;
140 }
141
142 i__2 = l;
143 for (j = 1; j <= i__2; ++j) {
144 f = d__[j];
145 g = e[j] + a[j + j * a_dim1] * f;
146 jp1 = j + 1;
147 if (l < jp1) {
148 goto L220;
149 }
150
151 i__3 = l;
152 for (k = jp1; k <= i__3; ++k) {
153 g += a[k + j * a_dim1] * d__[k];
154 e[k] += a[k + j * a_dim1] * f;
155 /* L200: */
156 }
157
158 L220:
159 e[j] = g;
160 /* L240: */
161 }
162 /* .......... FORM P .......... */
163 f = 0.;
164
165 i__2 = l;
166 for (j = 1; j <= i__2; ++j) {
167 e[j] /= h__;
168 f += e[j] * d__[j];
169 /* L245: */
170 }
171
172 h__ = f / (h__ + h__);
173 /* .......... FORM Q .......... */
174 i__2 = l;
175 for (j = 1; j <= i__2; ++j) {
176 /* L250: */
177 e[j] -= h__ * d__[j];
178 }
179 /* .......... FORM REDUCED A .......... */
180 i__2 = l;
181 for (j = 1; j <= i__2; ++j) {
182 f = d__[j];
183 g = e[j];
184
185 i__3 = l;
186 for (k = j; k <= i__3; ++k) {
187 /* L260: */
188 a[k + j * a_dim1] = a[k + j * a_dim1] - f * e[k] - g * d__[k];
189 }
190
191 /* L280: */
192 }
193
194 L285:
195 i__2 = l;
196 for (j = 1; j <= i__2; ++j) {
197 f = d__[j];
198 d__[j] = a[l + j * a_dim1];
199 a[l + j * a_dim1] = a[i__ + j * a_dim1];
200 a[i__ + j * a_dim1] = f * scale;
201 /* L290: */
202 }
203
204 L300:
205 ;
206 }
207
208 return 0;
209 } /* tred1_ */
210
211