Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_orthes.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int orthes_(integer *nm, integer *n, integer *low, integer *
00009 igh, doublereal *a, doublereal *ort)
00010 {
00011
00012 integer a_dim1, a_offset, i__1, i__2, i__3;
00013 doublereal d__1;
00014
00015
00016 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00017
00018
00019 static doublereal f, g, h__;
00020 static integer i__, j, m;
00021 static doublereal scale;
00022 static integer la, ii, jj, mp, kp1;
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
00068
00069 a_dim1 = *nm;
00070 a_offset = a_dim1 + 1;
00071 a -= a_offset;
00072 --ort;
00073
00074
00075 la = *igh - 1;
00076 kp1 = *low + 1;
00077 if (la < kp1) {
00078 goto L200;
00079 }
00080
00081 i__1 = la;
00082 for (m = kp1; m <= i__1; ++m) {
00083 h__ = 0.;
00084 ort[m] = 0.;
00085 scale = 0.;
00086
00087
00088 i__2 = *igh;
00089 for (i__ = m; i__ <= i__2; ++i__) {
00090
00091 scale += (d__1 = a[i__ + (m - 1) * a_dim1], abs(d__1));
00092 }
00093
00094 if (scale == 0.) {
00095 goto L180;
00096 }
00097 mp = m + *igh;
00098
00099 i__2 = *igh;
00100 for (ii = m; ii <= i__2; ++ii) {
00101 i__ = mp - ii;
00102 ort[i__] = a[i__ + (m - 1) * a_dim1] / scale;
00103 h__ += ort[i__] * ort[i__];
00104
00105 }
00106
00107 d__1 = sqrt(h__);
00108 g = -d_sign(&d__1, &ort[m]);
00109 h__ -= ort[m] * g;
00110 ort[m] -= g;
00111
00112 i__2 = *n;
00113 for (j = m; j <= i__2; ++j) {
00114 f = 0.;
00115
00116 i__3 = *igh;
00117 for (ii = m; ii <= i__3; ++ii) {
00118 i__ = mp - ii;
00119 f += ort[i__] * a[i__ + j * a_dim1];
00120
00121 }
00122
00123 f /= h__;
00124
00125 i__3 = *igh;
00126 for (i__ = m; i__ <= i__3; ++i__) {
00127
00128 a[i__ + j * a_dim1] -= f * ort[i__];
00129 }
00130
00131
00132 }
00133
00134 i__2 = *igh;
00135 for (i__ = 1; i__ <= i__2; ++i__) {
00136 f = 0.;
00137
00138 i__3 = *igh;
00139 for (jj = m; jj <= i__3; ++jj) {
00140 j = mp - jj;
00141 f += ort[j] * a[i__ + j * a_dim1];
00142
00143 }
00144
00145 f /= h__;
00146
00147 i__3 = *igh;
00148 for (j = m; j <= i__3; ++j) {
00149
00150 a[i__ + j * a_dim1] -= f * ort[j];
00151 }
00152
00153
00154 }
00155
00156 ort[m] = scale * ort[m];
00157 a[m + (m - 1) * a_dim1] = scale * g;
00158 L180:
00159 ;
00160 }
00161
00162 L200:
00163 return 0;
00164 }
00165