Go to the documentation of this file.00001
00006 #ifndef LU_solver_h
00007 #define LU_solver_h
00008
00009 #include <vector>
00010 #define TINY 1.0e-20
00011
00029 namespace LU_Solver {
00030
00038 int LUdecmpDoolittle(double *A, int n)
00039 {
00040 int i, j, k, l;
00041 double *pK, *pRow, *pCol;
00042
00043 for (k = 0, pK = A; k < n; pK += n, k++) {
00044 for (j = k; j < n; j++) {
00045 for (l = 0, pCol = A; l < k; pCol += n, l++)
00046 * (pK + j) -= *(pK + l) * *(pCol + j);
00047 }
00048 if ( *(pK + k) == 0.0 ) return -1;
00049 for (i = k + 1, pRow = pK + n; i < n; pRow += n, i++) {
00050 for (l = 0, pCol = A; l < k; pCol += n, l++)
00051 * (pRow + k) -= *(pRow + l) * *(pCol + k);
00052 *(pRow + k) /= *(pK + k);
00053 }
00054 }
00055 return 0;
00056 }
00057
00068 int LUsolveDoolittle(double *LU, double *X , int nRows, int nCols)
00069 {
00070 int i, j, k;
00071 int iRows, iCols;
00072
00073 for (i = 0; i < nRows; i++) {
00074 for (j = 0; j < nCols; j++) {
00075 iRows = i * nRows;
00076 iCols = i * nCols;
00077 if (*(LU + iCols + j) == 0.0) return -1;
00078 for (k = 0; k < i; k++)
00079 *(X + iCols + j) -= *(LU + iRows + k) * *(X + k * nCols + j);
00080 }
00081 }
00082 for (i = nRows - 1; i >= 0; i--) {
00083 for (j = nCols - 1; j >= 0; j--) {
00084 iRows = i * nRows;
00085 iCols = i * nCols;
00086 for (k = nRows - 1; k > i; k--)
00087 *(X + iCols + j) -= *(LU + iRows + k)* *(X + k * nCols + j);
00088 *(X + iCols + j) /= *(LU + iRows + i);
00089 }
00090 }
00091 return 0;
00092 }
00093
00103 int LUdecmpCrout(double * A, int n, int * indx)
00104 {
00105 int i;
00106 int imax = 0;
00107 int j, k;
00108 double big = 0;
00109 double dum, sum, temp;
00110 std::vector<double> vv(n);
00111
00112 for (i = 0; i < n; i++)
00113 *(indx + i) = i;
00114
00115 for (i = 0; i < n; i++) {
00116 for (j = 0; j < n; j++)
00117 if ((temp = fabs(*(A + i * n + j))) > big) big = temp;
00118 if (big == 0.0) return -1;
00119 vv[i] = 1.0 / big;
00120 }
00121 for (j = 0; j < n; j++) {
00122 for (i = 0; i < j; i++) {
00123 sum = *(A + i * n + j);
00124 for (k = 0; k < i; k++) sum -= *(A + i * n + k)* *(A + k * n + j);
00125 *(A + i * n + j) = sum;
00126 }
00127 big = 0.0;
00128 for (i = j; i < n; i++) {
00129 sum = *(A + i * n + j);
00130 for (k = 0; k < j; k++)
00131 sum -= *(A + i * n + k) * *(A + k * n + j);
00132 *(A + i * n + j) = sum;
00133 if ( (dum = vv[i] * fabs(sum)) >= big) {
00134 big = dum;
00135 imax = i;
00136 }
00137 }
00138 if (j != imax) {
00139 for (k = 0; k < n; k++) {
00140 dum = *(A + imax * n + k);
00141 *(A + imax * n + k) = *(A + j * n + k);
00142 *(A + j * n + k) = dum;
00143 }
00144
00145 temp = *(indx + j);
00146 *(indx + j) = *(indx + imax);
00147 *(indx + imax) = temp;
00148
00149 vv[imax] = vv[j];
00150 }
00151 if (*(A + j * n + j) == 0.0) *(A + j * n + j) = TINY;
00152
00153 if (j != n - 1) {
00154 dum = 1.0 / *(A + j * n + j);
00155 for (i = j + 1; i < n; i++) *(A + i * n + j) *= dum;
00156 }
00157 }
00158 return 0;
00159 }
00160
00172 int LUsolveCrout(double *LU, double *B, double *X, int nRows, int nCols, int *order)
00173 {
00174 int i, j, k;
00175 int iRows, iCols;
00176
00177 for (i = 0; i < nRows; i++) {
00178 for (j = 0; j < nCols; j++) {
00179 iRows = i * nRows;
00180 iCols = i * nCols;
00181 if (*(LU + iRows + i) == 0.0)
00182 return -1;
00183 *(X + iCols + j) = *(B + * (order + i) * nCols + j);
00184
00185 for (k = 0; k < i; k++)
00186 *(X + iCols + j) -= *(LU + iRows + k) * *(X + k * nCols + j);
00187 }
00188 }
00189
00190 for (i = nRows - 1; i >= 0; i--) {
00191 for (j = nCols - 1; j >= 0; j--) {
00192 iRows = i * nRows;
00193 iCols = i * nCols;
00194 for (k = nRows - 1; k > i; k--)
00195 *(X + iCols + j) -= *(LU + iRows + k)* *(X + k * nCols + j);
00196 *(X + iCols + j) /= *(LU + iRows + i);
00197 }
00198 }
00199 return 0;
00200 }
00201
00202 }
00203 #endif