Cholesky.h

Go to the documentation of this file.

#ifndef REDUKT_CHOLESKY_H
#define REDUKT_CHOLESKY_H

#include "math.h" // for sqrt()

#include <boost/numeric/ublas/matrix.hpp>  // replacement for TNT::Array2D
#include <boost/numeric/ublas/vector.hpp>  // replacement for TNT::Array1D

namespace ublas = boost::numeric::ublas;

template<class ValueType>
class Cholesky
{
        typedef ublas::matrix<ValueType> Matrix;
        typedef ublas::vector<ValueType> Vector;
        
        Matrix L_;              // lower triangular factor
        int isspd;              // 1 if matrix to be factored was SPD

public:

        //~ Cholesky();
        Cholesky(const Matrix &A);

        Matrix getL() const { return L_; }
        
        
        Vector solve( const Vector &b )
        {
                int n = L_.dim1();
                if (b.dim1() != n)
                        return Vector();


                Vector x = b;


                  // Solve L*y = b;
                  for (int k = 0; k < n; k++) 
                  {
                         for (int i = 0; i < k; i++) 
                                   x[k] -= x[i]*L_(k,i);
                         x[k] /= L_(k,k);
                        
                  }

                  // Solve L'*X = Y;
                  for (int k = n-1; k >= 0; k--) 
                  {
                         for (int i = k+1; i < n; i++) 
                                   x[k] -= x[i]*L_(i,k);
                         x[k] /= L_(k,k);
                  }

                return x;       
        }

        
        Matrix solve(const Matrix &B)
        {
                int n = L_.dim1();
                if (B.dim1() != n)
                        return Matrix();


                Matrix X = B;
                int nx = B.size2();

                // Solve L*y = b;
                  for (int j=0; j< nx; j++)
                  {
                        for (int k = 0; k < n; k++) 
                        {
                                for (int i = 0; i < k; i++) 
                                   X(k,j) -= X(i,j)*L_(k,i);
                                X(k,j) /= L_(k,k);
                         }
                  }

                // Solve L'*X = Y;
                 for (int j=0; j<nx; j++)
                 {
                        for (int k = n-1; k >= 0; k--) 
                        {
                                for (int i = k+1; i < n; i++) 
                                   X(k,j) -= X(i,j)*L_(i,k);
                                X(k,j) /= L_(k,k);
                        }
                  }
                  
                return X;               
        }

        int is_spd() const { return isspd; }

};

//~ template <class ValueType>
//~ Cholesky<ValueType>::Cholesky() : L_(0,0), isspd(0) {}


template <class ValueType>
Cholesky<ValueType>::Cholesky( const ublas::matrix<ValueType> &A )
{
        int m = A.size1();
        int n = A.size2();
        
        isspd = (m == n);

        if (m != n)
        {
                L_ = Matrix(0,0);
                return;
        }

        L_ = Matrix(n,n);


      // Main loop.
     for (int j = 0; j < n; j++) 
         {
        ValueType d(0.0);
        for (int k = 0; k < j; k++) 
                {
            ValueType s(0.0);
            for (int i = 0; i < k; i++) 
                        {
               s += L_(k,i)*L_(j,i);
            }
            L_(j,k) = s = (A(j,k) - s)/L_(k,k);
            d = d + s*s;
            isspd = isspd && (A(k,j) == A(j,k)); 
         }
         d = A(j,j) - d;
         isspd = isspd && (d > 0.0);
         L_(j,j) = sqrt(d > 0.0 ? d : 0.0);
         for (int k = j+1; k < n; k++) 
                 {
            L_(j,k) = 0.0;
         }
        }
}

#endif