MOAB
4.9.3pre
|
LU decomposition of a matrix with partial pivoting, and related features. More...
#include <PartialPivLU.h>
Public Types | |
enum | { MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime } |
typedef _MatrixType | MatrixType |
typedef SolverBase< PartialPivLU > | Base |
typedef PermutationMatrix < RowsAtCompileTime, MaxRowsAtCompileTime > | PermutationType |
typedef Transpositions < RowsAtCompileTime, MaxRowsAtCompileTime > | TranspositionType |
typedef MatrixType::PlainObject | PlainObject |
Public Member Functions | |
PartialPivLU () | |
Default Constructor. | |
PartialPivLU (Index size) | |
Default Constructor with memory preallocation. | |
template<typename InputType > | |
PartialPivLU (const EigenBase< InputType > &matrix) | |
template<typename InputType > | |
PartialPivLU & | compute (const EigenBase< InputType > &matrix) |
const MatrixType & | matrixLU () const |
const PermutationType & | permutationP () const |
template<typename Rhs > | |
const Solve< PartialPivLU, Rhs > | solve (const MatrixBase< Rhs > &b) const |
const Inverse< PartialPivLU > | inverse () const |
internal::traits< MatrixType > ::Scalar | determinant () const |
MatrixType | reconstructedMatrix () const |
Index | rows () const |
Index | cols () const |
template<typename RhsType , typename DstType > | |
EIGEN_DEVICE_FUNC void | _solve_impl (const RhsType &rhs, DstType &dst) const |
template<bool Conjugate, typename RhsType , typename DstType > | |
EIGEN_DEVICE_FUNC void | _solve_impl_transposed (const RhsType &rhs, DstType &dst) const |
Static Protected Member Functions | |
static void | check_template_parameters () |
Protected Attributes | |
MatrixType | m_lu |
PermutationType | m_p |
TranspositionType | m_rowsTranspositions |
Index | m_det_p |
bool | m_isInitialized |
LU decomposition of a matrix with partial pivoting, and related features.
_MatrixType | the type of the matrix of which we are computing the LU decomposition |
This class represents a LU decomposition of a square invertible matrix, with partial pivoting: the matrix A is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P is a permutation matrix.
Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.
The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided by class FullPivLU.
This is not a rank-revealing LU decomposition. Many features are intentionally absent from this class, such as rank computation. If you need these features, use class FullPivLU.
This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses in the general case. On the other hand, it is not suitable to determine whether a given matrix is invertible.
The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP().
Definition at line 62 of file PartialPivLU.h.
typedef SolverBase<PartialPivLU> Eigen::PartialPivLU< _MatrixType >::Base |
Reimplemented from Eigen::SolverBase< PartialPivLU< _MatrixType > >.
Definition at line 68 of file PartialPivLU.h.
typedef _MatrixType Eigen::PartialPivLU< _MatrixType >::MatrixType |
Definition at line 67 of file PartialPivLU.h.
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> Eigen::PartialPivLU< _MatrixType >::PermutationType |
Definition at line 75 of file PartialPivLU.h.
typedef MatrixType::PlainObject Eigen::PartialPivLU< _MatrixType >::PlainObject |
Definition at line 77 of file PartialPivLU.h.
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> Eigen::PartialPivLU< _MatrixType >::TranspositionType |
Definition at line 76 of file PartialPivLU.h.
anonymous enum |
Definition at line 71 of file PartialPivLU.h.
{ MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
Eigen::PartialPivLU< MatrixType >::PartialPivLU | ( | ) |
Default Constructor.
The default constructor is useful in cases in which the user intends to perform decompositions via PartialPivLU::compute(const MatrixType&).
Definition at line 254 of file PartialPivLU.h.
: m_lu(), m_p(), m_rowsTranspositions(), m_det_p(0), m_isInitialized(false) { }
Eigen::PartialPivLU< MatrixType >::PartialPivLU | ( | Index | size | ) | [explicit] |
Default Constructor with memory preallocation.
Like the default constructor but with preallocation of the internal data according to the specified problem size.
Definition at line 264 of file PartialPivLU.h.
: m_lu(size, size), m_p(size), m_rowsTranspositions(size), m_det_p(0), m_isInitialized(false) { }
Eigen::PartialPivLU< MatrixType >::PartialPivLU | ( | const EigenBase< InputType > & | matrix | ) | [explicit] |
Constructor.
matrix | the matrix of which to compute the LU decomposition. |
Definition at line 275 of file PartialPivLU.h.
: m_lu(matrix.rows(), matrix.rows()), m_p(matrix.rows()), m_rowsTranspositions(matrix.rows()), m_det_p(0), m_isInitialized(false) { compute(matrix.derived()); }
EIGEN_DEVICE_FUNC void Eigen::PartialPivLU< _MatrixType >::_solve_impl | ( | const RhsType & | rhs, |
DstType & | dst | ||
) | const [inline] |
Definition at line 191 of file PartialPivLU.h.
{ /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. * So we proceed as follows: * Step 1: compute c = Pb. * Step 2: replace c by the solution x to Lx = c. * Step 3: replace c by the solution x to Ux = c. */ eigen_assert(rhs.rows() == m_lu.rows()); // Step 1 dst = permutationP() * rhs; // Step 2 m_lu.template triangularView<UnitLower>().solveInPlace(dst); // Step 3 m_lu.template triangularView<Upper>().solveInPlace(dst); }
EIGEN_DEVICE_FUNC void Eigen::PartialPivLU< _MatrixType >::_solve_impl_transposed | ( | const RhsType & | rhs, |
DstType & | dst | ||
) | const [inline] |
Definition at line 213 of file PartialPivLU.h.
{ /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. * So we proceed as follows: * Step 1: compute c = Pb. * Step 2: replace c by the solution x to Lx = c. * Step 3: replace c by the solution x to Ux = c. */ eigen_assert(rhs.rows() == m_lu.cols()); if (Conjugate) { // Step 1 dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs); // Step 2 m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst); } else { // Step 1 dst = m_lu.template triangularView<Upper>().transpose().solve(rhs); // Step 2 m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst); } // Step 3 dst = permutationP().transpose() * dst; }
static void Eigen::PartialPivLU< _MatrixType >::check_template_parameters | ( | ) | [inline, static, protected] |
Definition at line 241 of file PartialPivLU.h.
Index Eigen::PartialPivLU< _MatrixType >::cols | ( | void | ) | const [inline] |
Reimplemented from Eigen::EigenBase< PartialPivLU< _MatrixType > >.
Definition at line 186 of file PartialPivLU.h.
{ return m_lu.cols(); }
PartialPivLU< MatrixType > & Eigen::PartialPivLU< MatrixType >::compute | ( | const EigenBase< InputType > & | matrix | ) |
Definition at line 462 of file PartialPivLU.h.
{ check_template_parameters(); // the row permutation is stored as int indices, so just to be sure: eigen_assert(matrix.rows()<NumTraits<int>::highest()); m_lu = matrix.derived(); eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); const Index size = matrix.rows(); m_rowsTranspositions.resize(size); typename TranspositionType::StorageIndex nb_transpositions; internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); m_det_p = (nb_transpositions%2) ? -1 : 1; m_p = m_rowsTranspositions; m_isInitialized = true; return *this; }
internal::traits< MatrixType >::Scalar Eigen::PartialPivLU< MatrixType >::determinant | ( | ) | const |
Definition at line 487 of file PartialPivLU.h.
{ eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return Scalar(m_det_p) * m_lu.diagonal().prod(); }
const Inverse<PartialPivLU> Eigen::PartialPivLU< _MatrixType >::inverse | ( | ) | const [inline] |
Definition at line 162 of file PartialPivLU.h.
{ eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return Inverse<PartialPivLU>(*this); }
const MatrixType& Eigen::PartialPivLU< _MatrixType >::matrixLU | ( | ) | const [inline] |
Definition at line 115 of file PartialPivLU.h.
{ eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return m_lu; }
const PermutationType& Eigen::PartialPivLU< _MatrixType >::permutationP | ( | ) | const [inline] |
Definition at line 123 of file PartialPivLU.h.
{ eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return m_p; }
MatrixType Eigen::PartialPivLU< MatrixType >::reconstructedMatrix | ( | ) | const |
Definition at line 497 of file PartialPivLU.h.
{ eigen_assert(m_isInitialized && "LU is not initialized."); // LU MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix() * m_lu.template triangularView<Upper>(); // P^{-1}(LU) res = m_p.inverse() * res; return res; }
Index Eigen::PartialPivLU< _MatrixType >::rows | ( | void | ) | const [inline] |
Reimplemented from Eigen::EigenBase< PartialPivLU< _MatrixType > >.
Definition at line 185 of file PartialPivLU.h.
{ return m_lu.rows(); }
const Solve<PartialPivLU, Rhs> Eigen::PartialPivLU< _MatrixType >::solve | ( | const MatrixBase< Rhs > & | b | ) | const [inline] |
This method returns the solution x to the equation Ax=b, where A is the matrix of which *this is the LU decomposition.
b | the right-hand-side of the equation to solve. Can be a vector or a matrix, the only requirement in order for the equation to make sense is that b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition. |
Example:
Output:
Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution theoretically exists and is unique regardless of b.
Reimplemented from Eigen::SolverBase< PartialPivLU< _MatrixType > >.
Definition at line 149 of file PartialPivLU.h.
{ eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return Solve<PartialPivLU, Rhs>(*this, b.derived()); }
Index Eigen::PartialPivLU< _MatrixType >::m_det_p [protected] |
Definition at line 249 of file PartialPivLU.h.
bool Eigen::PartialPivLU< _MatrixType >::m_isInitialized [protected] |
Definition at line 250 of file PartialPivLU.h.
MatrixType Eigen::PartialPivLU< _MatrixType >::m_lu [protected] |
Definition at line 246 of file PartialPivLU.h.
PermutationType Eigen::PartialPivLU< _MatrixType >::m_p [protected] |
Definition at line 247 of file PartialPivLU.h.
TranspositionType Eigen::PartialPivLU< _MatrixType >::m_rowsTranspositions [protected] |
Definition at line 248 of file PartialPivLU.h.