GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add internal method _solve_impl_transposed() to LU decomposition classes that solves A^T x = b or A^* x = b.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2015-11-30 13:39:24 -08:00
parent 274b2272b7
commit 1663d15da7
3 changed files with 148 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
Eigen/src/LU/FullPivLU.h
View File

@ -389,6 +389,10 @@ template<typename _MatrixType> class FullPivLU
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;
#endif
protected:
@ -753,6 +757,70 @@ void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
for(Index i = nonzero_pivots; i < m_lu.cols(); ++i)
dst.row(permutationQ().indices().coeff(i)).setZero();
}
template<typename _MatrixType>
template<bool Conjugate, typename RhsType, typename DstType>
void FullPivLU<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
{
/* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1},
* and since permutations are real and unitary, we can write this
* as A^T = Q U^T L^T P,
* So we proceed as follows:
* Step 1: compute c = Q^T rhs.
* Step 2: replace c by the solution x to U^T x = c. May or may not exist.
* Step 3: replace c by the solution x to L^T x = c.
* Step 4: result = P^T c.
* If Conjugate is true, replace "^T" by "^*" above.
*/
const Index rows = this->rows(), cols = this->cols(),
nonzero_pivots = this->rank();
eigen_assert(rhs.rows() == cols);
const Index smalldim = (std::min)(rows, cols);
if(nonzero_pivots == 0)
{
dst.setZero();
return;
}
typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
// Step 1
c = permutationQ().inverse() * rhs;
if (Conjugate) {
// Step 2
m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()
.adjoint()
.solveInPlace(c.topRows(nonzero_pivots));
// Step 3
m_lu.topLeftCorner(smalldim, smalldim)
.template triangularView<UnitLower>()
.adjoint()
.solveInPlace(c.topRows(smalldim));
} else {
// Step 2
m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()
.transpose()
.solveInPlace(c.topRows(nonzero_pivots));
// Step 3
m_lu.topLeftCorner(smalldim, smalldim)
.template triangularView<UnitLower>()
.transpose()
.solveInPlace(c.topRows(smalldim));
}
// Step 4
PermutationPType invp = permutationP().inverse().eval();
for(Index i = 0; i < smalldim; ++i)
dst.row(invp.indices().coeff(i)) = c.row(i);
for(Index i = smalldim; i < rows; ++i)
dst.row(invp.indices().coeff(i)).setZero();
}
#endif
namespace internal {

27
Eigen/src/LU/PartialPivLU.h
View File

@ -208,6 +208,33 @@ template<typename _MatrixType> class PartialPivLU
// Step 3
m_lu.template triangularView<Upper>().solveInPlace(dst);
}
template<bool Conjugate, typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC
void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
/* 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;
}
#endif
protected:

36
test/lu.cpp
View File

@ -92,6 +92,20 @@ template<typename MatrixType> void lu_non_invertible()
// test that the code, which does resize(), may be applied to an xpr
m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
// test solve with transposed
m3 = MatrixType::Random(rows,cols2);
m2 = m1.transpose()*m3;
m3 = MatrixType::Random(rows,cols2);
lu.template _solve_impl_transposed<false>(m2, m3);
VERIFY_IS_APPROX(m2, m1.transpose()*m3);
// test solve with conjugate transposed
m3 = MatrixType::Random(rows,cols2);
m2 = m1.adjoint()*m3;
m3 = MatrixType::Random(rows,cols2);
lu.template _solve_impl_transposed<true>(m2, m3);
VERIFY_IS_APPROX(m2, m1.adjoint()*m3);
}
template<typename MatrixType> void lu_invertible()
@ -124,6 +138,12 @@ template<typename MatrixType> void lu_invertible()
m2 = lu.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
VERIFY_IS_APPROX(m2, lu.inverse()*m3);
// test solve with transposed
lu.template _solve_impl_transposed<false>(m3, m2);
VERIFY_IS_APPROX(m3, m1.transpose()*m2);
// test solve with conjugate transposed
lu.template _solve_impl_transposed<true>(m3, m2);
VERIFY_IS_APPROX(m3, m1.adjoint()*m2);
// Regression test for Bug 302
MatrixType m4 = MatrixType::Random(size,size);
@ -136,14 +156,24 @@ template<typename MatrixType> void lu_partial_piv()
PartialPivLU.h
*/
typedef typename MatrixType::Index Index;
Index rows = internal::random<Index>(1,4);
Index cols = rows;
Index size = internal::random<Index>(1,4);
MatrixType m1(cols, rows);
MatrixType m1(size, size), m2(size, size), m3(size, size);
m1.setRandom();
PartialPivLU<MatrixType> plu(m1);
VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
m3 = MatrixType::Random(size,size);
m2 = plu.solve(m3);
VERIFY_IS_APPROX(m3, m1*m2);
VERIFY_IS_APPROX(m2, plu.inverse()*m3);
// test solve with transposed
plu.template _solve_impl_transposed<false>(m3, m2);
VERIFY_IS_APPROX(m3, m1.transpose()*m2);
// test solve with conjugate transposed
plu.template _solve_impl_transposed<true>(m3, m2);
VERIFY_IS_APPROX(m3, m1.adjoint()*m2);
}
template<typename MatrixType> void lu_verify_assert()