# Ublas :

From: Robbie Vogt (r.vogt_at_[hidden])
Date: 2005-07-25 20:19:59

I've written a cholesky factorisation in a similar fashion to the LU
factorisation already in uBLAS. At the moment I have only considered
dense matrix types, but possibly someone might have an implemtnation
appropriate for other types.

Robbie.

Vania Dos Santos Eleuterio wrote:

>Hi, i would like to know if there is some possibility to compute
>cholesky factorization using ublas.
>I'm very new in using C++ and also in using ublas.
>
>Thank you very much
>Vania
>
>_______________________________________________
>ublas mailing list
>ublas_at_[hidden]

>
>
>

// Cholesky factorizations in the spirit of lu.hpp

#ifndef BOOST_UBLAS_CHOLESKY_H
#define BOOST_UBLAS_CHOLESKY_H

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/triangular.hpp>

namespace boost { namespace numeric { namespace ublas {

// Cholesky factorization
template<class M>
bool cholesky_factorize (M &m) {
typedef M matrix_type;
typedef BOOST_UBLAS_TYPENAME M::size_type size_type;
typedef BOOST_UBLAS_TYPENAME M::value_type value_type;

BOOST_UBLAS_CHECK (m.size1() == m.size2(), external_logic("Cholesky decomposition is only valid for a square, positive definite matrix."));

size_type size = m.size1();
vector<value_type> d(size);

for (size_type i = 0; i < size; ++ i) {
matrix_row<M> mri (row (m, i));
for (size_type j = i; j < size; ++ j) {
matrix_row<M> mrj (row (m, j));

value_type elem = m(i,j) - inner_prod(project(mri,range(0,i)), project(mrj,range(0,i)));

if (i == j) {
if (elem <= 0.0) {
// matrix after rounding errors is not positive definite
return false;
}
else {
d(i) = type_traits<value_type>::sqrt(elem);
}
}
else {
m(j,i) = elem / d(i);
}
}
}
// put the diagonal back in
for (size_type i = 0; i < size; ++ i) {
m(i,i) = d(i);
}
// decomposition succeeded
return true;
}

// Cholesky substitution
template<class M, class E>
void cholesky_substitute (const M &m, vector_expression<E> &e) {
typedef const M const_matrix_type;
typedef vector<typename E::value_type> vector_type;

inplace_solve (m, e, lower_tag ());
inplace_solve (trans(m), e, upper_tag ());
}
template<class M, class E>
void cholesky_substitute (const M &m, matrix_expression<E> &e) {
typedef const M const_matrix_type;
typedef matrix<typename E::value_type> matrix_type;

inplace_solve (m, e, lower_tag ());
inplace_solve (trans(m), e, upper_tag ());
}
template<class E, class M>
void cholesky_substitute_left (vector_expression<E> &e, const M &m) {
typedef const M const_matrix_type;
typedef vector<typename E::value_type> vector_type;

inplace_solve (trans(m), e, upper_tag ());
inplace_solve (m, e, lower_tag ());
}
template<class E, class M>
void cholesky_substitute_left (matrix_expression<E> &e, const M &m) {
typedef const M const_matrix_type;
typedef matrix<typename E::value_type> matrix_type;

inplace_solve (trans(m), e, upper_tag ());
inplace_solve (m, e, lower_tag ());
}

// Cholesky matrix inversion
template<class M>
void cholesky_invert (M &m)
{
typedef typename M::size_type size_type;
typedef typename M::value_type value_type;

size_type size = m.size1();

// determine the inverse of the lower traingular matrix
for (size_type i = 0; i < size; ++ i) {
m(i,i) = 1 / m(i,i);

for (size_type j = i+1; j < size; ++ j) {
value_type elem(0);

for (size_type k = i; k < j; ++ k) {
elem -= m(j,k)*m(k,i);
}
m(j,i) = elem / m(j,j);
}
}

// multiply the upper and lower inverses together