Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_reduc.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int reduc_(integer *nm, integer *n, doublereal *a,
00009 doublereal *b, doublereal *dl, integer *ierr)
00010 {
00011
00012 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00013
00014
00015 double sqrt(doublereal);
00016
00017
00018 static integer i__, j, k;
00019 static doublereal x, y;
00020 static integer i1, j1, nn;
00021
00022
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
00070
00071
00072
00073
00074
00075
00076
00077 --dl;
00078 b_dim1 = *nm;
00079 b_offset = b_dim1 + 1;
00080 b -= b_offset;
00081 a_dim1 = *nm;
00082 a_offset = a_dim1 + 1;
00083 a -= a_offset;
00084
00085
00086 *ierr = 0;
00087 nn = abs(*n);
00088 if (*n < 0) {
00089 goto L100;
00090 }
00091
00092 i__1 = *n;
00093 for (i__ = 1; i__ <= i__1; ++i__) {
00094 i1 = i__ - 1;
00095
00096 i__2 = *n;
00097 for (j = i__; j <= i__2; ++j) {
00098 x = b[i__ + j * b_dim1];
00099 if (i__ == 1) {
00100 goto L40;
00101 }
00102
00103 i__3 = i1;
00104 for (k = 1; k <= i__3; ++k) {
00105
00106 x -= b[i__ + k * b_dim1] * b[j + k * b_dim1];
00107 }
00108
00109 L40:
00110 if (j != i__) {
00111 goto L60;
00112 }
00113 if (x <= 0.) {
00114 goto L1000;
00115 }
00116 y = sqrt(x);
00117 dl[i__] = y;
00118 goto L80;
00119 L60:
00120 b[j + i__ * b_dim1] = x / y;
00121 L80:
00122 ;
00123 }
00124 }
00125
00126
00127 L100:
00128 i__2 = nn;
00129 for (i__ = 1; i__ <= i__2; ++i__) {
00130 i1 = i__ - 1;
00131 y = dl[i__];
00132
00133 i__1 = nn;
00134 for (j = i__; j <= i__1; ++j) {
00135 x = a[i__ + j * a_dim1];
00136 if (i__ == 1) {
00137 goto L180;
00138 }
00139
00140 i__3 = i1;
00141 for (k = 1; k <= i__3; ++k) {
00142
00143 x -= b[i__ + k * b_dim1] * a[j + k * a_dim1];
00144 }
00145
00146 L180:
00147 a[j + i__ * a_dim1] = x / y;
00148
00149 }
00150 }
00151
00152 i__1 = nn;
00153 for (j = 1; j <= i__1; ++j) {
00154 j1 = j - 1;
00155
00156 i__2 = nn;
00157 for (i__ = j; i__ <= i__2; ++i__) {
00158 x = a[i__ + j * a_dim1];
00159 if (i__ == j) {
00160 goto L240;
00161 }
00162 i1 = i__ - 1;
00163
00164 i__3 = i1;
00165 for (k = j; k <= i__3; ++k) {
00166
00167 x -= a[k + j * a_dim1] * b[i__ + k * b_dim1];
00168 }
00169
00170 L240:
00171 if (j == 1) {
00172 goto L280;
00173 }
00174
00175 i__3 = j1;
00176 for (k = 1; k <= i__3; ++k) {
00177
00178 x -= a[j + k * a_dim1] * b[i__ + k * b_dim1];
00179 }
00180
00181 L280:
00182 a[i__ + j * a_dim1] = x / dl[i__];
00183
00184 }
00185 }
00186
00187 goto L1001;
00188
00189 L1000:
00190 *ierr = *n * 7 + 1;
00191 L1001:
00192 return 0;
00193 }
00194