25 #include <boost/numeric/ublas/vector.hpp>
26 #include <boost/numeric/ublas/vector_proxy.hpp>
28 #include <boost/numeric/ublas/matrix.hpp>
29 #include <boost/numeric/ublas/matrix_proxy.hpp>
32 namespace ublas = boost::numeric::ublas;
44 template <
class MATRIX,
class TRIA >
47 using namespace ublas;
49 typedef typename MATRIX::value_type T;
51 assert( A.size1() == A.size2() );
52 assert( A.size1() == L.size1() );
53 assert( A.size2() == L.size2() );
55 const size_t n = A.size1();
57 for (
size_t k=0 ; k < n; k++ )
60 double qL_kk = A( k,k ) - inner_prod(
project( row( L, k ), range( 0, k ) ),
61 project( row( L, k ), range( 0, k ) ) );
70 double L_kk = sqrt( qL_kk );
73 matrix_column<TRIA> cLk( L, k );
75 = (
project( column( A, k ), range( k+1, n ) )
76 - prod(
project( L, range( k+1, n ), range( 0, k ) ),
77 project( row( L, k ), range( 0, k ) ) ) ) / L_kk;