linear_doolittle.hpp

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: doolittle.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_LINEAR_DOOLITTLE_HPP
#define realroot_SOLVE_SBDSLV_LINEAR_DOOLITTLE_HPP
//--------------------------------------------------------------------
#include <realroot/system_epsilon.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace linear
{

  template< typename real_t >
  void Lsolve(real_t * x , real_t const  * const  L, real_t const *const  b,  unsigned n )
  {
    for ( unsigned i = 0; i  < n ; i ++ )
      {
    real_t sum = 0.0;
    for ( unsigned j = 0; j < i; j ++ ) sum += L[i*n+j]*x[j];
    x[i] = (b[i]-sum);
      };
  };
  
  template< typename real_t >
  bool Usolve( real_t * x, real_t const * const U, real_t const *const  b,  unsigned n )
  {
    for ( int  i = n-1; i  >=0  ; i -- )
      {
        real_t sum = 0.0;
        for ( int j = n-1; j > i; j -- ) sum += U[i*n+j]*x[j];
    if ( std::abs(U[i*n+i]) < numerics::epsilon<real_t>::result*n*n) return false;
    x[i] = (b[i]-sum)/U[i*n+i];
      };
    return true;
  };
  
  template< typename real_t >
  bool doolittle( real_t const * const  A, real_t * L, real_t * U, unsigned n )
  {
    unsigned k,i,p;
    for ( k = 0; k < n; k++ )
      {
    //  valueof(k);
    for ( i = k; i < n; i ++ )
      {
        U[ k*n + i ] = A[ k*n + i ];
        for ( p = 0; p < k; p ++ ) U[ k*n + i ] -= L[ k*n + p ]*U[ p*n + i ];
      };
    for ( i = k + 1; i < n; i ++  )
      {
        if ( std::abs( U[ k*n + k ] ) < numerics::epsilon<real_t>::result*n*n ) return false;
        real_t sum = 0.0;
        for ( p = 0; p < k; p ++ ) sum += L[i*n+p]*U[p*n+k];
        L[i*n+k] = (A[i*n + k] - sum)/U[ k*n + k ];
      };
      };
    return true;
  };
  
  template< typename real_t >
  bool LUsolve( real_t * x, real_t const * const M, real_t const * const  b, unsigned n )
  {
    real_t tmp[ n*n ];
    real_t xtmp[n] ;
    if ( !doolittle( M, tmp, tmp, n ) ) return false;
    Lsolve( xtmp, tmp, b, n );
    return Usolve( x, tmp, xtmp, n );
  };

  template< typename real_t >
  bool LUinverse( real_t * iM, real_t const * const M, unsigned n )
  {
    real_t lu[n*n];
    real_t xtmp[n] ;
    real_t tmp[n];
    real_t _b[n+1];
    real_t * b = _b+1;
    std::fill( b, b + n , (real_t)0.0 );
    if ( ! doolittle( M, lu, lu, n ) ) return false;
    for ( unsigned c = 0; c < n; c ++ )
      {
    if (c > 0) b[c-1] = 0;
    b[c] = 1;
    Lsolve( xtmp, lu, b, n );
    if ( !Usolve( tmp, lu, xtmp, n ) ) return false;
    for ( unsigned i = 0; i < n; i ++ ) iM[ i*n + c ] = tmp[i]; 
      };
    return true;
  }; 

};
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //