Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_corth.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int corth_(integer *nm, integer *n, integer *low, integer *
00009 igh, doublereal *ar, doublereal *ai, doublereal *ortr, doublereal *
00010 orti)
00011 {
00012
00013 integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3;
00014 doublereal d__1, d__2;
00015
00016
00017 double sqrt(doublereal);
00018
00019
00020 static doublereal f, g, h__;
00021 static integer i__, j, m;
00022 static doublereal scale;
00023 static integer la;
00024 static doublereal fi;
00025 static integer ii, jj;
00026 static doublereal fr;
00027 static integer mp;
00028 extern doublereal pythag_(doublereal *, doublereal *);
00029 static integer kp1;
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
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 ai_dim1 = *nm;
00082 ai_offset = ai_dim1 + 1;
00083 ai -= ai_offset;
00084 ar_dim1 = *nm;
00085 ar_offset = ar_dim1 + 1;
00086 ar -= ar_offset;
00087 --orti;
00088 --ortr;
00089
00090
00091 la = *igh - 1;
00092 kp1 = *low + 1;
00093 if (la < kp1) {
00094 goto L200;
00095 }
00096
00097 i__1 = la;
00098 for (m = kp1; m <= i__1; ++m) {
00099 h__ = 0.;
00100 ortr[m] = 0.;
00101 orti[m] = 0.;
00102 scale = 0.;
00103
00104
00105 i__2 = *igh;
00106 for (i__ = m; i__ <= i__2; ++i__) {
00107
00108 scale = scale + (d__1 = ar[i__ + (m - 1) * ar_dim1], abs(d__1)) +
00109 (d__2 = ai[i__ + (m - 1) * ai_dim1], abs(d__2));
00110 }
00111
00112 if (scale == 0.) {
00113 goto L180;
00114 }
00115 mp = m + *igh;
00116
00117 i__2 = *igh;
00118 for (ii = m; ii <= i__2; ++ii) {
00119 i__ = mp - ii;
00120 ortr[i__] = ar[i__ + (m - 1) * ar_dim1] / scale;
00121 orti[i__] = ai[i__ + (m - 1) * ai_dim1] / scale;
00122 h__ = h__ + ortr[i__] * ortr[i__] + orti[i__] * orti[i__];
00123
00124 }
00125
00126 g = sqrt(h__);
00127 f = pythag_(&ortr[m], &orti[m]);
00128 if (f == 0.) {
00129 goto L103;
00130 }
00131 h__ += f * g;
00132 g /= f;
00133 ortr[m] = (g + 1.) * ortr[m];
00134 orti[m] = (g + 1.) * orti[m];
00135 goto L105;
00136
00137 L103:
00138 ortr[m] = g;
00139 ar[m + (m - 1) * ar_dim1] = scale;
00140
00141 L105:
00142 i__2 = *n;
00143 for (j = m; j <= i__2; ++j) {
00144 fr = 0.;
00145 fi = 0.;
00146
00147 i__3 = *igh;
00148 for (ii = m; ii <= i__3; ++ii) {
00149 i__ = mp - ii;
00150 fr = fr + ortr[i__] * ar[i__ + j * ar_dim1] + orti[i__] * ai[
00151 i__ + j * ai_dim1];
00152 fi = fi + ortr[i__] * ai[i__ + j * ai_dim1] - orti[i__] * ar[
00153 i__ + j * ar_dim1];
00154
00155 }
00156
00157 fr /= h__;
00158 fi /= h__;
00159
00160 i__3 = *igh;
00161 for (i__ = m; i__ <= i__3; ++i__) {
00162 ar[i__ + j * ar_dim1] = ar[i__ + j * ar_dim1] - fr * ortr[i__]
00163 + fi * orti[i__];
00164 ai[i__ + j * ai_dim1] = ai[i__ + j * ai_dim1] - fr * orti[i__]
00165 - fi * ortr[i__];
00166
00167 }
00168
00169
00170 }
00171
00172 i__2 = *igh;
00173 for (i__ = 1; i__ <= i__2; ++i__) {
00174 fr = 0.;
00175 fi = 0.;
00176
00177 i__3 = *igh;
00178 for (jj = m; jj <= i__3; ++jj) {
00179 j = mp - jj;
00180 fr = fr + ortr[j] * ar[i__ + j * ar_dim1] - orti[j] * ai[i__
00181 + j * ai_dim1];
00182 fi = fi + ortr[j] * ai[i__ + j * ai_dim1] + orti[j] * ar[i__
00183 + j * ar_dim1];
00184
00185 }
00186
00187 fr /= h__;
00188 fi /= h__;
00189
00190 i__3 = *igh;
00191 for (j = m; j <= i__3; ++j) {
00192 ar[i__ + j * ar_dim1] = ar[i__ + j * ar_dim1] - fr * ortr[j]
00193 - fi * orti[j];
00194 ai[i__ + j * ai_dim1] = ai[i__ + j * ai_dim1] + fr * orti[j]
00195 - fi * ortr[j];
00196
00197 }
00198
00199
00200 }
00201
00202 ortr[m] = scale * ortr[m];
00203 orti[m] = scale * orti[m];
00204 ar[m + (m - 1) * ar_dim1] = -g * ar[m + (m - 1) * ar_dim1];
00205 ai[m + (m - 1) * ai_dim1] = -g * ai[m + (m - 1) * ai_dim1];
00206 L180:
00207 ;
00208 }
00209
00210 L200:
00211 return 0;
00212 }
00213