From 682b2ef17ee360ac11cebbe7286dc4edd9accfa3 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Thu, 8 Jun 2017 15:06:27 +0200
Subject: [PATCH] bug #1423: fix LSCG\'s Jacobi preconditioner for row-major
 matrices.

---
 .../BasicPreconditioners.h                    | 27 ++++++++++++++-----
 test/lscg.cpp                                 |  8 ++++++
 2 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
index 358444aff..279c9173c 100644
--- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
+++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
@@ -152,13 +152,28 @@ class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar>
     {
       // 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;
diff --git a/test/lscg.cpp b/test/lscg.cpp
index daa62a954..d49ee00c3 100644
--- a/test/lscg.cpp
+++ b/test/lscg.cpp
@@ -14,12 +14,20 @@ template<typename T> void test_lscg_T()
 {
   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()