Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_tred2.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int tred2_(integer *nm, integer *n, doublereal *a,
00009 doublereal *d__, doublereal *e, doublereal *z__)
00010 {
00011
00012 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
00013 doublereal d__1;
00014
00015
00016 double d_sign(doublereal *, doublereal *);
00017
00018
00019 static doublereal f, g, h__;
00020 static integer i__, j, k, l;
00021 static doublereal scale, hh;
00022 static integer ii, jp1;
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 z_dim1 = *nm;
00068 z_offset = z_dim1 + 1;
00069 z__ -= z_offset;
00070 --e;
00071 --d__;
00072 a_dim1 = *nm;
00073 a_offset = a_dim1 + 1;
00074 a -= a_offset;
00075
00076
00077 i__1 = *n;
00078 for (i__ = 1; i__ <= i__1; ++i__) {
00079
00080 i__2 = *n;
00081 for (j = i__; j <= i__2; ++j) {
00082
00083 z__[j + i__ * z_dim1] = a[j + i__ * a_dim1];
00084 }
00085
00086 d__[i__] = a[*n + i__ * a_dim1];
00087
00088 }
00089
00090 if (*n == 1) {
00091 goto L510;
00092 }
00093
00094 i__1 = *n;
00095 for (ii = 2; ii <= i__1; ++ii) {
00096 i__ = *n + 2 - ii;
00097 l = i__ - 1;
00098 h__ = 0.;
00099 scale = 0.;
00100 if (l < 2) {
00101 goto L130;
00102 }
00103
00104 i__2 = l;
00105 for (k = 1; k <= i__2; ++k) {
00106
00107 scale += (d__1 = d__[k], abs(d__1));
00108 }
00109
00110 if (scale != 0.) {
00111 goto L140;
00112 }
00113 L130:
00114 e[i__] = d__[l];
00115
00116 i__2 = l;
00117 for (j = 1; j <= i__2; ++j) {
00118 d__[j] = z__[l + j * z_dim1];
00119 z__[i__ + j * z_dim1] = 0.;
00120 z__[j + i__ * z_dim1] = 0.;
00121
00122 }
00123
00124 goto L290;
00125
00126 L140:
00127 i__2 = l;
00128 for (k = 1; k <= i__2; ++k) {
00129 d__[k] /= scale;
00130 h__ += d__[k] * d__[k];
00131
00132 }
00133
00134 f = d__[l];
00135 d__1 = sqrt(h__);
00136 g = -d_sign(&d__1, &f);
00137 e[i__] = scale * g;
00138 h__ -= f * g;
00139 d__[l] = f - g;
00140
00141 i__2 = l;
00142 for (j = 1; j <= i__2; ++j) {
00143
00144 e[j] = 0.;
00145 }
00146
00147 i__2 = l;
00148 for (j = 1; j <= i__2; ++j) {
00149 f = d__[j];
00150 z__[j + i__ * z_dim1] = f;
00151 g = e[j] + z__[j + j * z_dim1] * f;
00152 jp1 = j + 1;
00153 if (l < jp1) {
00154 goto L220;
00155 }
00156
00157 i__3 = l;
00158 for (k = jp1; k <= i__3; ++k) {
00159 g += z__[k + j * z_dim1] * d__[k];
00160 e[k] += z__[k + j * z_dim1] * f;
00161
00162 }
00163
00164 L220:
00165 e[j] = g;
00166
00167 }
00168
00169 f = 0.;
00170
00171 i__2 = l;
00172 for (j = 1; j <= i__2; ++j) {
00173 e[j] /= h__;
00174 f += e[j] * d__[j];
00175
00176 }
00177
00178 hh = f / (h__ + h__);
00179
00180 i__2 = l;
00181 for (j = 1; j <= i__2; ++j) {
00182
00183 e[j] -= hh * d__[j];
00184 }
00185
00186 i__2 = l;
00187 for (j = 1; j <= i__2; ++j) {
00188 f = d__[j];
00189 g = e[j];
00190
00191 i__3 = l;
00192 for (k = j; k <= i__3; ++k) {
00193
00194 z__[k + j * z_dim1] = z__[k + j * z_dim1] - f * e[k] - g *
00195 d__[k];
00196 }
00197
00198 d__[j] = z__[l + j * z_dim1];
00199 z__[i__ + j * z_dim1] = 0.;
00200
00201 }
00202
00203 L290:
00204 d__[i__] = h__;
00205
00206 }
00207
00208 i__1 = *n;
00209 for (i__ = 2; i__ <= i__1; ++i__) {
00210 l = i__ - 1;
00211 z__[*n + l * z_dim1] = z__[l + l * z_dim1];
00212 z__[l + l * z_dim1] = 1.;
00213 h__ = d__[i__];
00214 if (h__ == 0.) {
00215 goto L380;
00216 }
00217
00218 i__2 = l;
00219 for (k = 1; k <= i__2; ++k) {
00220
00221 d__[k] = z__[k + i__ * z_dim1] / h__;
00222 }
00223
00224 i__2 = l;
00225 for (j = 1; j <= i__2; ++j) {
00226 g = 0.;
00227
00228 i__3 = l;
00229 for (k = 1; k <= i__3; ++k) {
00230
00231 g += z__[k + i__ * z_dim1] * z__[k + j * z_dim1];
00232 }
00233
00234 i__3 = l;
00235 for (k = 1; k <= i__3; ++k) {
00236 z__[k + j * z_dim1] -= g * d__[k];
00237
00238 }
00239 }
00240
00241 L380:
00242 i__3 = l;
00243 for (k = 1; k <= i__3; ++k) {
00244
00245 z__[k + i__ * z_dim1] = 0.;
00246 }
00247
00248
00249 }
00250
00251 L510:
00252 i__1 = *n;
00253 for (i__ = 1; i__ <= i__1; ++i__) {
00254 d__[i__] = z__[*n + i__ * z_dim1];
00255 z__[*n + i__ * z_dim1] = 0.;
00256
00257 }
00258
00259 z__[*n + *n * z_dim1] = 1.;
00260 e[1] = 0.;
00261 return 0;
00262 }
00263