nojhan/cholesky
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

beautifying comments, missing headers

Browse source
This commit is contained in:
nojhan 2011-12-15 23:25:42 +01:00
commit 541b8bdc6b
2 changed files with 8 additions and 4 deletions
Download patch file Download diff file Expand all files Collapse all files

10
cholesky.h
View file

@ -21,7 +21,9 @@ Authors:
Caner Candan <caner.candan@thalesgroup.com> Caner Candan <caner.candan@thalesgroup.com>
*/ */
#include <string>
#include <stdexcept> #include <stdexcept>
#include <cmath>
#include <boost/numeric/ublas/symmetric.hpp> #include <boost/numeric/ublas/symmetric.hpp>
using namespace boost::numeric; using namespace boost::numeric;
@ -277,11 +279,12 @@ public:
} }
/** Compute the final symetric matrix: _L = L D^1/2
* remember that V = ( L D^1/2) ( L D^1/2)^T
* the factorization is thus L*D^1/2
*/
inline typename Cholesky<T>::FactorMat root( typename Cholesky<T>::FactorMat& L, typename Cholesky<T>::FactorMat& D ) inline typename Cholesky<T>::FactorMat root( typename Cholesky<T>::FactorMat& L, typename Cholesky<T>::FactorMat& D )
{ {
// now compute the final symetric matrix: this->_L = L D^1/2
// remember that V = ( L D^1/2) ( L D^1/2)^T
// fortunately, the square root of a diagonal matrix is the square // fortunately, the square root of a diagonal matrix is the square
// root of all its elements // root of all its elements
typename Cholesky<T>::FactorMat sqrt_D = D; typename Cholesky<T>::FactorMat sqrt_D = D;
@ -289,7 +292,6 @@ public:
sqrt_D(i,i) = sqrt(D(i,i)); sqrt_D(i,i) = sqrt(D(i,i));
} }
// the factorization is thus L*D^1/2
return ublas::prod( L, sqrt_D ); return ublas::prod( L, sqrt_D );
} }
}; };

2
test.cpp
View file

@ -161,6 +161,8 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig
} }
// a variance-covariance matrix of size N*N // a variance-covariance matrix of size N*N
// Note: a covariance matrix is necessarily semi-positive definite
// thus, any failure in the Cholesky factorization is due to round-off errors
CovarMat V = ublas::prod( ublas::trans(S), S ); CovarMat V = ublas::prod( ublas::trans(S), S );
assert( V.size1() == N && V.size2() == N ); assert( V.size1() == N && V.size2() == N );
#ifndef NDEBUG #ifndef NDEBUG