cholesky.h

高效的c++科学算法库

/*                                                           -*- c++ -*-
 ***********************************************************************
 * Scientific Library (GNU Public Licence)
 *
 * module sparse
 *
 * Author: Pierre Aubert [e-mail]
 *
 * $Id: cholesky.h,v 1.6 1998/07/16 14:45:46 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_SSS_H
#define SL_ALGEBRA_CHOLESKY_SSS_H

#ifdef HAVE_NAMESPACE
namespace sl 
{
#endif

    template <typename T_value, typename T_storage>
    void
    cholesky( Matrix<T_value,SymmetricSkyline,T_storage> & m )
    {
        assert( m.isFactored() == false );
        m.setIsFactored();

        Index i,j,k;
        UInt di,dj,dk,ind,nmax;
        T_value my_datadk;

        Matrix<UInt,TVector> const& my_diag = m.diag();
        Matrix<UInt,TVector> const& my_haut = m.haut();
    
        for( i=1 ; i<=my_diag.rows() ; ++i )
            {
                di = my_diag(i);
                for( j=i-my_haut(i)+1 ; j<=i-1  ; ++j )
                    {
                        dj  = my_diag(j);
                        ind = di+j-i;
                        nmax= max(j-my_haut(j),i-my_haut(i));
                        for( k=nmax+j ; k<=j-1  ; ++k )
                            {
                                m(i,j) -= m(i,k)*m(k,j);
                            } // end for k 
                          // STOP
                        assert( abs(m(j,j)) > FLT_EPSILON );
                        m(i,j) /= m(j,j);
                    } // end for j 
                for( k=i+1-my_haut(i) ; k<=i-i  ; ++k )
                    {
                        dk = di+k-i;
                        my_datadk = m(i,k);
                        m(i,i) -= my_datadk*my_datadk;
                    } // end for k 
                  // STOP
                assert( m(i,i) > FLT_EPSILON );
                m(i,i)=sqrt(m(i,i));
            } // end for i 

    }

    template <typename T_value, typename T_storage>
    void choleskyLowerSolve( Matrix<T_value,SymmetricSkyline,T_storage> const& m,
                     Matrix<T_value,TVector> const& B,
                     Matrix<T_value,TVector>      & X )
    {
        assert( m.isFactored() == true );

        Index i,j;
        Matrix<UInt,TVector> const& my_diag = m.diag();
        Matrix<UInt,TVector> const& my_haut = m.haut();
    
    
        for( i=1 ; i<=my_diag.rows() ; ++i )
            {
                X(i) = B(i);
                for( j=i-my_haut(i)+1 ; j<=i-1  ; ++j )
                    {
                        X(i) -= m(i,j)*X(j);
                    } // end for j 
                X(i) /= m(i,i);
            } // end for i 

    }

    template <typename T_value, typename T_storage>
    void choleskyUpperSolve( Matrix<T_value,SymmetricSkyline,T_storage> const& m ,
                     Matrix<T_value,TVector> & X )
    {
        assert( m.isFactored() == true );

        Index i,j,n;
        Matrix<UInt,TVector> const& my_diag = m.diag();
        Matrix<UInt,TVector> const& my_haut = m.haut();
    

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

    }

    template <typename T_value, typename T_storage>
    void choleskySolve( Matrix<T_value,SymmetricSkyline,T_storage> & m ,
                        Matrix<T_value,TVector> const& B,
                        Matrix<T_value,TVector>      & X )
    {
        try {
            if( ! m.isFactored() )
                cholesky(m);

        }
        catch( ErrorCode eMatrixIsNotPositiveDefinite )
            {
                throw;
            }
        
        choleskyLowerSolve(m,B,X);
        choleskyUpperSolve(m,X);
    }

#ifdef HAVE_NAMESPACE
}
#endif

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