commit 682b2ef17ee360ac11cebbe7286dc4edd9accfa3
author Gael Guennebaud <g.gael@free.fr> Thu Jun 08 15:06:27 2017 +0200
committer Gael Guennebaud <g.gael@free.fr> Thu Jun 08 15:06:27 2017 +0200
tree 32711f0ce5bf53ecb87d0578fcd3ea61fa4b6c36
parent 4bbc32046810f65bb0f77f6dbe538abad51de281

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

Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h

diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
index 358444a..279c917 100644
--- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
+++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
@@ -152,13 +152,28 @@
     {
       // Compute the inverse squared-norm of each column of mat
       m_invdiag.resize(mat.cols());
-      for(Index j=0; j<mat.outerSize(); ++j)
+      if(MatType::IsRowMajor)
       {
-        RealScalar sum = mat.innerVector(j).squaredNorm();
-        if(sum>0)
-          m_invdiag(j) = RealScalar(1)/sum;
-        else
-          m_invdiag(j) = RealScalar(1);
+        m_invdiag.setZero();
+        for(Index j=0; j<mat.outerSize(); ++j)
+        {
+          for(typename MatType::InnerIterator it(mat,j); it; ++it)
+            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;
       return *this;

test/lscg.cpp

diff --git a/test/lscg.cpp b/test/lscg.cpp
index daa62a9..d49ee00 100644
--- a/test/lscg.cpp
+++ b/test/lscg.cpp
@@ -14,12 +14,20 @@
 {
   LeastSquaresConjugateGradient<SparseMatrix<T> > lscg_colmajor_diag;
   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_I)     );
   
   CALL_SUBTEST( check_sparse_leastsquare_solving(lscg_colmajor_diag)  );
   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()