---

jama_cholesky_old.h

Go to the documentation of this file.

#ifndef JAMA_CHOLESKY_H
#define JAMA_CHOLESKY_H

#include "math.h"
        /* needed for sqrt() below. */


namespace JAMA
{


template <class Real>
class Cholesky
{
        TNT::Array2D<Real> L;           // lower triangular factor
        int isspd;                              // 1 if matrix to be factored was SPD


public:

        Cholesky();
        Cholesky(const TNT::Array2D<Real> &A);
        TNT::Array2D<Real> getL() const;
        int is_spd() const;

};

template <class Real>
Cholesky<Real>::Cholesky() : L(0,0), isspd(0) {}

template <class Real>
int Cholesky<Real>::is_spd() const
{
        return isspd;
}

template <class Real>
TNT::Array2D<Real> Cholesky<Real>::getL() const
{
        return L;
}

template <class Real>
Cholesky<Real>::Cholesky(const TNT::Array2D<Real> &A)
{


        int m = A.dim1();
        int n = A.dim2();
        
        isspd = (m == n);

        if (m != n)
        {
                L = TNT::Array2D<Real>(0,0);
                return;
        }

        L = TNT::Array2D<Real>(n,n);


      // Main loop.
     for (int j = 0; j < n; j++) 
         {
        double d = 0.0;
        for (int k = 0; k < j; k++) 
                {
            Real 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;
         }
        }
}

}
// namespace JAMA

#endif
// JAMA_CHOLESKY_H

---