24 #include "LUDecomposition.h"
28 LUDecomposition::LUDecomposition(
const MatrixDouble &a) : sing(false){
30 debugLog.setProceedingText(
"[DEBUG LUDecomposition]");
31 errorLog.setProceedingText(
"[ERROR LUDecomposition]");
32 warningLog.setProceedingText(
"[WARNING LUDecomposition]");
39 const double TINY=1.0e-40;
40 unsigned int i,imax,j,k;
48 if ((temp=fabs( lu[i][j] )) > big) big=temp;
51 errorLog <<
"Error in LUDecomposition constructor, big == 0.0" << endl;
59 temp=vv[i]*fabs(lu[i][k]);
68 lu[imax][j] = lu[k][j];
75 if (lu[k][k] == 0.0) lu[k][k] = TINY;
76 for (i=k+1; i<N; i++) {
77 temp = lu[i][k] /= lu[k][k];
79 lu[i][j] -= temp * lu[k][j];
85 LUDecomposition::~LUDecomposition(){
89 bool LUDecomposition::solve_vector(
const VectorDouble &b,VectorDouble &x)
91 int i=0,ii=0,ip=0,j=0;
95 if (b.size() != N || x.size() != N){
96 errorLog <<
"solve_vector(const VectorDouble &b,VectorDouble &x) - the size of the two vectors does not match!" << endl;
99 for (i=0;i<n;i++) x[i] = b[i];
105 for (j=ii-1;j<i;j++) sum -= lu[i][j] * x[j];
110 for (i=n-1;i>=0;i--) {
112 for (j=i+1;j<n;j++) sum -= lu[i][j] * x[j];
113 x[i] = sum / lu[i][i];
119 bool LUDecomposition::solve(
const MatrixDouble &b,MatrixDouble &x)
122 if (b.getNumRows() != N || x.getNumRows() != N || b.getNumCols() != x.getNumCols() ){
123 errorLog <<
"solve(const MatrixDouble &b,MatrixDouble &x) - the size of the two matrices does not match!" << endl;
127 for (
unsigned int j=0; j<m; j++) {
128 for(
unsigned int i=0; i<N; i++) xx[i] = b[i][j];
130 for(
unsigned int i=0; i<N; i++) x[i][j] = xx[i];
135 bool LUDecomposition::inverse(MatrixDouble &ainv)
140 for (j=0;j<N;j++) ainv[i][j] = 0.0;
143 return solve(ainv,ainv);
146 double LUDecomposition::det()
149 for (
unsigned int i=0;i<N;i++) dd *= lu[i][i];
153 bool LUDecomposition::mprove(
const VectorDouble &b,VectorDouble &x)
158 long double sdp = -b[i];
160 sdp += (
long double) aref[i][j] * (
long double)x[j];
163 if( !solve_vector(r,r) ){
166 for (i=0;i<N;i++) x[i] -= r[i];
170 bool LUDecomposition::getIsSingular(){
174 MatrixDouble LUDecomposition::getLU(){
unsigned int getNumCols() const