bug #1423: fix LSCG\'s Jacobi preconditioner for row-major matrices.

This commit is contained in:
Gael Guennebaud 2017-06-08 15:06:27 +02:00
parent 4bbc320468
commit 682b2ef17e
2 changed files with 29 additions and 6 deletions

View File

@ -152,13 +152,28 @@ class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar>
{ {
// Compute the inverse squared-norm of each column of mat // Compute the inverse squared-norm of each column of mat
m_invdiag.resize(mat.cols()); m_invdiag.resize(mat.cols());
for(Index j=0; j<mat.outerSize(); ++j) if(MatType::IsRowMajor)
{ {
RealScalar sum = mat.innerVector(j).squaredNorm(); m_invdiag.setZero();
if(sum>0) for(Index j=0; j<mat.outerSize(); ++j)
m_invdiag(j) = RealScalar(1)/sum; {
else for(typename MatType::InnerIterator it(mat,j); it; ++it)
m_invdiag(j) = RealScalar(1); m_invdiag(it.index()) += it.value();
}
for(Index j=0; j<mat.cols(); ++j)
if(m_invdiag(j)>0)
m_invdiag(j) = RealScalar(1)/m_invdiag(j);
}
else
{
for(Index j=0; j<mat.outerSize(); ++j)
{
RealScalar sum = mat.innerVector(j).squaredNorm();
if(sum>0)
m_invdiag(j) = RealScalar(1)/sum;
else
m_invdiag(j) = RealScalar(1);
}
} }
Base::m_isInitialized = true; Base::m_isInitialized = true;
return *this; return *this;

View File

@ -14,12 +14,20 @@ template<typename T> void test_lscg_T()
{ {
LeastSquaresConjugateGradient<SparseMatrix<T> > lscg_colmajor_diag; LeastSquaresConjugateGradient<SparseMatrix<T> > lscg_colmajor_diag;
LeastSquaresConjugateGradient<SparseMatrix<T>, IdentityPreconditioner> lscg_colmajor_I; LeastSquaresConjugateGradient<SparseMatrix<T>, IdentityPreconditioner> lscg_colmajor_I;
LeastSquaresConjugateGradient<SparseMatrix<T,RowMajor> > lscg_rowmajor_diag;
LeastSquaresConjugateGradient<SparseMatrix<T,RowMajor>, IdentityPreconditioner> lscg_rowmajor_I;
CALL_SUBTEST( check_sparse_square_solving(lscg_colmajor_diag) ); CALL_SUBTEST( check_sparse_square_solving(lscg_colmajor_diag) );
CALL_SUBTEST( check_sparse_square_solving(lscg_colmajor_I) ); CALL_SUBTEST( check_sparse_square_solving(lscg_colmajor_I) );
CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_diag) ); CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_diag) );
CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_I) ); CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_I) );
CALL_SUBTEST( check_sparse_square_solving(lscg_rowmajor_diag) );
CALL_SUBTEST( check_sparse_square_solving(lscg_rowmajor_I) );
CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_rowmajor_diag) );
CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_rowmajor_I) );
} }
void test_lscg() void test_lscg()