24 #define GRT_DLL_EXPORTS
30 debugLog.setProceedingText(
"[DEBUG LUdcmp]");
31 errorLog.setProceedingText(
"[ERROR LUdcmp]");
32 warningLog.setProceedingText(
"[WARNING LUdcmp]");
37 Cholesky::Cholesky(
MatrixFloat &a) : N(a.getNumRows()), el(a) {
39 debugLog.setProceedingText(
"[DEBUG LUdcmp]");
40 errorLog.setProceedingText(
"[ERROR LUdcmp]");
41 warningLog.setProceedingText(
"[WARNING LUdcmp]");
47 if( el.getNumCols() != N ){
48 errorLog <<
"The input matrix is not square!" << std::endl;
55 for (sum=el[i][j],k=i-1;k>=0;k--) sum -= el[i][k]*el[j][k];
58 errorLog <<
"Sum is <=0.0" << std::endl;
62 }
else el[j][i]=sum/el[i][i];
77 if (b.size() != N || x.size() != N){
78 errorLog <<
":solve(vector<Float> &b, vector<Float> &x) - The input vectors are not the same size!" << std::endl;
82 for(sum=b[i],k=i-1;k>=0;k--) sum -= el[i][k]*x[k];
85 for (i=n-1; i>=0; i--) {
86 for (sum=x[i],k=i+1;k<n;k++) sum -= el[k][i]*x[k];
94 if (b.size() != N || y.size() != N){
95 errorLog <<
"elmult(vector<Float> &y vector<Float> &b) - The input vectors are not the same size!" << std::endl;
100 for (j=0; j<=i; j++) b[i] += el[i][j]*y[j];
109 if (b.size() != N || y.size() != N){
110 errorLog <<
"elsolve(vector<Float> &b vector<Float> &y) - The input vectors are not the same size!" << std::endl;
113 for (i=0; i<N; i++) {
114 for (sum=b[i],j=0; j<i; j++) sum -= el[i][j]*y[j];
122 const int n = int(N);
126 for(i=0; i<n; i++)
for(j=0; j<=i; j++){
127 sum = (i==j? 1. : 0.);
128 for(k=i-1; k>=j; k--) sum -= el[i][k]*ainv[j][k];
129 ainv[j][i]= sum/el[i][i];
131 for(i=n-1; i>=0; i--)
for(j=0; j<=i; j++){
132 sum = (i<j? 0. : ainv[j][i]);
133 for(k=i+1; k<n; k++) sum -= el[k][i]*ainv[j][k];
134 ainv[i][j] = ainv[j][i] = sum/el[i][i];
139 Float Cholesky::logdet(){
141 for(
unsigned int i=0; i<N; i++) sum += log(el[i][i]);
This code is based on the LU Decomposition code from Numerical Recipes (3rd Edition) ...
virtual bool resize(const unsigned int r, const unsigned int c)