shells/utilities/CholeskyDecomp.h

Go to the documentation of this file.

// -*- C++ -*- 
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//                                   Fehmi Cirak
//                        California Institute of Technology
//                           (C) 2004 All Rights Reserved
//
// <LicenseText>
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
#ifndef CHOLESKYDECOMP_H
#define CHOLESKYDECOMP_H

#include <iostream>
#include <cassert>
#include <cmath>


namespace utilities {
    template<int ND>
    void choleskyDecomposition(double a[ND][ND], double p[ND]);

    template <int ND>
    void choleskySolve(double a[ND][ND], double p[ND], double rhs[ND], double x[ND]);
}


// implementation
namespace utilities {

template <int ND>
void choleskyDecomposition(double a[ND][ND], double p[ND]) 
{
    // compute cholesky decompositon a = l \cdot l^T  
    // only upper triangle of the symmetric matrix a need to be provided
    // on output l is returned in lower triangle of a and 
    // its diagonale elements in p[1..ND]       
    double sum;
    for (int i=0; i<ND; ++i) {
        for (int j=i; j<ND; ++j) {
            
            int k;
            for (sum=a[i][j], k=i-1; k>-1; --k) sum -= a[i][k]*a[j][k]; 
            
            if (i==j) {
                if (sum<0.0) {
                    std::cerr << "choleskyDecomposition non positive definite matrix." 
                              << std::endl;
                    assert(false);
                }
                p[i] = std::sqrt(sum);
                
            } else a[j][i] = sum/p[i];
        }       
    }    
    return;
}
        
    
template <int ND>
void choleskySolve(double a[ND][ND], double p[ND], double rhs[ND], double x[ND]) 
{
    // a and p are the outputs of choleskyDecomposition()
    // b is the right hand side and x the unknown 
    double sum;
    int k;
    
    // solve l \cdot y = b storing y in x 
    for (int i=0; i<ND; ++i) {
        for (sum=rhs[i], k=i-1; k>-1; --k) sum -= a[i][k]*x[k];
        x[i] = sum/p[i];
    }
        
    // solve l^T \cdot x = y
    for (int i=(ND-1); i>-1; --i) {
        for (sum =x[i], k=i+1; k<ND; ++k) sum -= a[k][i]*x[k];
        x[i] = sum/p[i];
    }   
    return;
}

} // namespace utilities

#endif

---