Workaround fixing aliasing issue in x = SparseLU::solve(x)

(transplanted from 46dd1bb1be
)
This commit is contained in:
Gael Guennebaud 2013-11-15 11:19:19 +01:00
parent 4c5da3b03a
commit 516304cd90
2 changed files with 5 additions and 2 deletions

View File

@ -553,7 +553,8 @@ struct permut_matrix_product_retval
template<typename Dest> inline void evalTo(Dest& dst) const
{
const Index n = Side==OnTheLeft ? rows() : cols();
// FIXME we need an is_same for expression that is not sensitive to constness. For instance
// is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
{
// apply the permutation inplace

View File

@ -229,8 +229,10 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ
// Permute the right hand side to form X = Pr*B
// on return, X is overwritten by the computed solution
X.resize(B.rows(),B.cols());
// this ugly const_cast_derived() helps to detect aliasing when applying the permutations
for(Index j = 0; j < B.cols(); ++j)
X.col(j) = rowsPermutation() * B.col(j);
X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
//Forward substitution with L
this->matrixL().solveInPlace(X);