/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/math/big/big_backward_gauss_seidel.h

Go to the documentation of this file.

#ifndef OPENTISSUE_CORE_MATH_BIG_BACKWARD_GAUSS_SEIDEL_H
#define OPENTISSUE_CORE_MATH_BIG_BACKWARD_GAUSS_SEIDEL_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/core/math/big/big_types.h>
#include <OpenTissue/core/math/math_is_number.h>
#include <OpenTissue/core/math/math_value_traits.h>

#include <stdexcept>

namespace OpenTissue
{
  namespace math
  {
    namespace big
    {

      template<typename T>
      inline void backward_gauss_seidel(
        boost::numeric::ublas::compressed_matrix<T> const & A
        , boost::numeric::ublas::vector<T>       & x
        , boost::numeric::ublas::vector<T> const & b
        )
      {
        typedef boost::numeric::ublas::compressed_matrix<T> matrix_type;
        typedef boost::numeric::ublas::vector<T>            vector_type;
        typedef typename vector_type::size_type             size_type;
        typedef typename vector_type::value_type            value_type;
        typedef OpenTissue::math::ValueTraits<T>            value_traits;

        using std::fabs;

        if(A.size1() <= 0 || A.size2() <= 0)
          throw std::invalid_argument("backward_gauss_seidel(): A was empty");

        if(b.size() != A.size1())
          throw std::invalid_argument("backward_gauss_seidel(): The size of b must be the same as the number of rows in A");

        if(x.size() != A.size2())
          throw std::invalid_argument("backward_gauss_seidel(): The size of x must be the same as the number of columns in A");

        size_type const N = A.filled1() - 1;
        for (size_type i = 0; i < N; ++i)
        {
          size_type row = N - i - 1;

          assert( ( row >= 0 && row< N )        || !"backward_gauss_seidel::do_apply(): row index were out of range");

          size_type begin  = A.index1_data()[row];
          size_type end    = A.index1_data()[row + 1];

          value_type sum   = b(row);
          value_type diag  = value_traits::zero();

          for (size_type j = begin; j < end; ++ j)
          {
            size_type  const & col  = A.index2_data()[j];
            value_type const & A_ij = A.value_data()[j];

            assert( ( col >= 0 && col< A.size2() )       || !"backward_gauss_seidel::do_apply(): column index were out of range");
            assert( is_number( A_ij )                    || !"backward_gauss_seidel::do_apply(): A_ij value was not a number?");

            // 2007-06-04 kenny: Performance issue, the if-statement should
            //                   be on the outside of the for-loop!
            if(row==col)
              diag = A_ij;
            else
              sum -= A_ij * x(  col   );
          }
          // 2007-06-10 kenny: Yikes how should we handle a diagonal zero-value?
          assert( fabs(diag) > value_traits::zero() || !"backward_gauss_seidel::do_apply(): Diagonal were zero");
          assert( is_number( sum )                  || !"backward_gauss_seidel::do_apply(): sum value was not a number?");
          assert( is_number( diag )                 || !"backward_gauss_seidel::do_apply(): diag value was not a number?");
          x (row) = sum / diag;
          assert( is_number( x(row ) )              || !"backward_gauss_seidel::do_apply(): updated value was not a number?");
        }

      }

      template<typename T>
      inline void backward_gauss_seidel(
        boost::numeric::ublas::compressed_matrix<T> const & A
        , boost::numeric::ublas::vector<T>       & x
        , boost::numeric::ublas::vector<T> const & b
        , size_t                           const & max_iterations
        , size_t                                 & iterations
        )
      {
        if(max_iterations < 1)
          throw std::invalid_argument("backward_gauss_seidel(): max_iterations must be a positive number");

        iterations = 0;
        while(iterations<max_iterations)
        {
          ++iterations;
          backward_gauss_seidel(A,x,b);
        }
      }

      class BackwardGaussSeidelFunctor
      {
      public:

        template<typename T>
        void operator()(
          boost::numeric::ublas::compressed_matrix<T> const & A
          , boost::numeric::ublas::vector<T>       & x
          , boost::numeric::ublas::vector<T> const & b
          , size_t                           const & max_iterations
          , size_t                                 & iterations
          )
        {
          backward_gauss_seidel(A,x,b,max_iterations,iterations);
        }

        template<typename T>
        void operator()(
          boost::numeric::ublas::compressed_matrix<T> const & A
          , boost::numeric::ublas::vector<T>       & x
          , boost::numeric::ublas::vector<T> const & b
          )
        {
          backward_gauss_seidel(A,x,b);
        }

      };

    } // end of namespace big
  } // end of namespace math
} // end of namespace OpenTissue

// OPENTISSUE_CORE_MATH_BIG_BACKWARD_GAUSS_SEIDEL_H
#endif