lu.h

高效的c++科学算法库

/*                                                           -*- c++ -*-
 ***********************************************************************
 * Scientific Library (GNU Public Licence)
 *
 * module sparse
 *
 * Author: Pierre Aubert [e-mail]
 *
 * $Id: lu.h,v 1.2 1998/07/28 11:40:17 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_LU_DIAGONAL_H
#define SL_ALGEBRA_LU_DIAGONAL_H

#ifdef HAVE_NAMESPACE
namespace sl 
{
#endif

    template <typename T_value, typename T_storage>
    void lu( Matrix<T_value,Diagonal,T_storage> & m, Matrix<Index,TVector> & indx)
    {
        typedef T_value value_t;
        
          // declarations
        Index j;
        Index n=m.n();
        Index p=m.p();
        
          // check for size
        if( n != p )
            {
                cerr << "sl/algebra/diagonal/lu.h: lu() matrix must be square" << endl;
                throw eMatrixMustBeSquare;
            }
        
        for( j=1; j<= n; ++j )
            {
                if( m(j,j) == 0 )
                    {
                        cerr << "sl/algebra/diagonal/lu.h: lu() matrix must be positive definite" << endl;
                        throw eMatrixIsNotPositiveDefinite;
                    }
                else
                    m(j,j) = 1.0/m(j,j);
            }
        
    } // end lu()
        

    template <typename T_value, typename T_storage>
    void luLowerSolve( Matrix<T_value,Diagonal,T_storage> const& A,
                       Matrix<Index,TVector>              const& indx,
                       Matrix<T_value,TVector>            const& B,
                       Matrix<T_value,TVector>                 & X )
    {
        X=A*B;
    }

    template <typename T_value, typename T_storage>
    void luUpperSolve( Matrix<T_value,Diagonal,T_storage> const& A ,
                       Matrix<Index,TVector>              const& indx,
                       Matrix<T_value,TVector>                 & X )
    {
    }


    template <typename T_value, typename T_storage>
    void luSolve( Matrix<T_value,Diagonal,T_storage> & A ,
                  Matrix<Index,TVector>              & indx,
                  Matrix<T_value,TVector>       const& B,
                  Matrix<T_value,TVector>            & X )
    {
        try {
            lu(A,indx);
        }
        catch( ErrorCode error ) 
            {
                if( error == eMatrixIsNotPositiveDefinite )
                    cerr << "sl/algebra/diagonal/lu.h: luSolve() matrix is not positive definite" << endl;
                else if( error == eMatrixMustBeSquare )
                    cerr << "sl/algebra/diagonal/lu.h: luSolve() matrix must be square" << endl;
                else
                    cerr << "sl/algebra/diagonal/lu.h: luSolve() unknown exception" << endl;
                throw;
            }
        
        luLowerSolve(A,indx,B,X);
        luUpperSolve(A,indx,X);
    }

#ifdef HAVE_NAMESPACE
}
#endif

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