linearSolvers.hpp

/* $Id: linearSolver.hpp  $ */

#ifndef LINEAR_SOLVERS_HPP
#define LINEAR_SOLVERS_HPP

#include "jafarConfig.h"

#ifdef HAVE_BOOST_SANDBOX
#ifdef HAVE_LAPACK

#include "jmath/jblas.hpp"

namespace jafar {
  namespace jmath {

    namespace LinearSolvers {
      
      int solve_QR(jblas::mat_column_major& A, const jblas::mat_column_major& B, jblas::mat_column_major& X);

      int solve_QR_noQ(jblas::mat_column_major A, jblas::mat B, jblas::mat& X);

      /*
       * This function returns the solution of Ax=b
       *
       * The function assumes that A is symmetric & positive definite and employs
       * the Cholesky decomposition:
       * If A=U^T U with U upper triangular, the system to be solved becomes
       * (U^T U) x = b
       * This amount to solving U^T y = b for y and then U x = y for x
       *
       * The function returns 0 in case of error, 1 if successfull
       */
      int solve_Cholesky(jblas::mat_column_major A, jblas::mat_column_major B, jblas::mat& X);
      int solve_Cholesky(jblas::mat_column_major A, jblas::vec B, jblas::vec& x);

      int solve_LU(jblas::mat_column_major A, jblas::vec b, jblas::vec& x);
      //int solve_LU(const jblas::mat& A, const jblas::vec& b, jblas::vec& x);

      /*
       * This function returns the solution of Ax = b
       *
       * The function is based on SVD decomposition:
       * If A=U D V^T with U, V orthogonal and D diagonal, the linear system becomes
       * (U D V^T) x = b or x=V D^{-1} U^T b
       * Note that V D^{-1} U^T is the pseudoinverse A^+
       *
       * The function returns 0 in case of error, 1 if successfull
       */
      int solve_SVD(const jblas::mat& A, const jblas::vec& b, jblas::vec& x);

      /*
       * This function returns the solution of Ax = b for a real symmetric matrix A
       *
       * The function is based on Bunch-Kaufman factorization:
       * A is factored as U*D*U^T where U is upper triangular and
       * D symmetric and block diagonal
       *
       * The function returns 0 in case of error, 1 if successfull
       */
      int solve_BK(const jblas::mat& A, const jblas::vec& b, jblas::vec& x);
    } // namespace LinearSolvers
  } // namespace jmath
} // namespace jafar


#endif // HAVE_LAPACK
#endif // HAVE_BOOST_SANDBOX

#endif // JMATH_LINEAR_LEAST_SQUARES_HPP