cholesky.h

高效的c++科学算法库

/*                                                           -*- c++ -*-
 ***********************************************************************
 * Scientific Library (GNU Public Licence)
 *
 * module sparse
 *
 * Author: Pierre Aubert [e-mail]
 *
 * $Id: cholesky.h,v 1.2 1998/07/17 15:09:47 paubert Exp $
 *
 * Suggestions: [e-mail]
 * Bugs:   [e-mail] 
 *
 * For more information, please see the sl++ Home Page:
 * http://wwwinfo.cern.ch/~ldeniau/sl.html 
 *
 ***********************************************************************
 */

#ifndef SL_ALGEBRA_CHOLESKY_SYMMETRIC_H
#define SL_ALGEBRA_CHOLESKY_SYMMETRIC_H

#ifdef HAVE_NAMESPACE
namespace sl 
{
#endif

    template <typename T_value, typename T_storage>
    void cholesky( Matrix<T_value,Symmetric,T_storage> & m )
    {
        Index i,j,k,n,p;

        n = m.n();
        p = m.p();
        
        assert( n == p );

        for( i=1 ; i<=n ; ++i )
            {
                for( j=1 ; j<=i-1  ; ++j )
                    {
                        for( k=1 ; k<=j-1  ; ++k )
                            m(i,j) -= m(i,k)*m(k,j);
                        m(i,j) /= m(j,j);
                    }
                for( k=1 ; k<=i-i  ; ++k )
                    m(i,i) -= m(i,k)*m(i,k);
                m(i,i)=sqrt(m(i,i));
            } 

    } // end cholkesky()
        

    template <typename T_value, typename T_storage>
    void choleskyLowerSolve( Matrix<T_value,Symmetric,T_storage> const& m,
                     Matrix<T_value,TVector> const& B,
                     Matrix<T_value,TVector>      & X )
    {
        Index i,j;
    
        for( i=1 ; i<=m.n() ; ++i )
            {
                X(i) = B(i);
                for( j=1 ; j<=i-1  ; ++j )
                    X(i) -= m(i,j)*X(j);
                X(i) /= m(i,i);
            }
    }

    template <typename T_value, typename T_storage>
    void choleskyUpperSolve( Matrix<T_value,Symmetric,T_storage> const& m ,
                     Matrix<T_value,TVector> & X )
    {
        Index i,j,n=m.n();

        for( i=n ; i>=1 ; --i )
            {
                for( j=i+1 ; j<=n  ; ++j )
                    X(i) -= m(i,j)*X(j);
                X(i) /= m(i,i);
            } // end for i 

    }


    template <typename T_value, typename T_storage>
    void choleskySolve( Matrix<T_value,Symmetric,T_storage> & m ,
                        Matrix<T_value,TVector> const& B,
                        Matrix<T_value,TVector>      & X )
    {
        cholesky(m);

        choleskyLowerSolve(m,B,X);
        choleskyUpperSolve(m,X);
    }

#ifdef HAVE_NAMESPACE
}
#endif

#endif
// end of SL_ALGEBRA_CHOLESKY_SYMMETRIC_H definition(s)
// do not write anything after this line!