Tools/Math/LA_Lu_decomp.h

Go to the documentation of this file.

//------------------------------------------------------------------------------
/**
 * @file LA_Lu_decomp.h
 * Contains class LinAlg::LU
 *
 * @author <a href="mailto:stefanuhrig@gmx.net">Stefan Uhrig</a>
 */
//------------------------------------------------------------------------------
#ifndef LINALG_LU_H_INCLUDED
#define LINALG_LU_H_INCLUDED
//------------------------------------------------------------------------------
#include <vector>
#include "LA_Matrix.h"
//------------------------------------------------------------------------------
namespace LA
{
//------------------------------------------------------------------------------
/**
 * Performs a LU-decomposition of a matrix.
 * Adapted from Template Numerical Toolkit (TNT) (http://math.nist.gov/tnt/)
 *
 * <P>
 * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
 * unit lower triangular matrix L, an n-by-n upper triangular matrix U,
 * and a permutation vector piv of length m so that A(piv,:) = L*U.
 * If m < n, then L is m-by-m and U is m-by-n.
 * <P>
 * The LU decompostion with pivoting always exists, even if the matrix is
 * singular, so the constructor will never fail.  The primary use of the
 * LU decomposition is in the solution of square systems of simultaneous
 * linear equations.
 * 
 * @author <a href="mailto:stefanuhrig@gmx.net">Stefan Uhrig</a>
 */
class LU
{
public:
  /**
   * Constructs the LU-decomposition of a the passed matrix.
   * @param M Matrix to decompose.
   */
  LU(const Matrix& M);

private:
  LU(const LU& lu) {}
  LU& operator=(const LU& lu) { return *this; }

public:
  /**
   * True if the matrix is singular.
   * @return True if the matrix is singular.
   */
  bool singular() const;
  
  /**
   * Returns the determinant of the matrix.
   * @return The determinant of the matrix.
   */
  double det() const;
  
  /**
   * Solves the equation Ax=b.
   * @param b A vector.
   * @return The solution vector x.
   */
  Vector solve(const Vector& b) const;
  
  /**
   * Calculates the inverse of the matrix.
   * @return The inverse of the matrix.
   */
  Matrix inverse() const;

private:
  Vector permute_copy(const Vector& A, const std::vector<int>& piv) const;

private:
  Matrix LU_;
  int m, n;
  int pivsign;
  std::vector<int> piv;
};
//------------------------------------------------------------------------------
}
//------------------------------------------------------------------------------
#endif
//------------------------------------------------------------------------------

---