From ff46ec0f240ef84e2293b33b265c703e9b765c2e Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 22 Sep 2014 23:33:28 +0200
Subject: [PATCH] bug #881: make SparseMatrixBase::isApprox(SparseMatrixBase)
 exploits sparse computations instead of converting the operands to dense
 matrices.

---
 Eigen/src/SparseCore/SparseFuzzy.h      | 30 ++++++++++++++-----------
 Eigen/src/SparseCore/SparseMatrixBase.h |  3 +--
 test/sparse_basic.cpp                   |  4 ++++
 3 files changed, 22 insertions(+), 15 deletions(-)

diff --git a/Eigen/src/SparseCore/SparseFuzzy.h b/Eigen/src/SparseCore/SparseFuzzy.h
index 45f36e9eb..3e67cbf5f 100644
--- a/Eigen/src/SparseCore/SparseFuzzy.h
+++ b/Eigen/src/SparseCore/SparseFuzzy.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -10,17 +10,21 @@
 #ifndef EIGEN_SPARSE_FUZZY_H
 #define EIGEN_SPARSE_FUZZY_H
 
-// template<typename Derived>
-// template<typename OtherDerived>
-// bool SparseMatrixBase<Derived>::isApprox(
-//   const OtherDerived& other,
-//   typename NumTraits<Scalar>::Real prec
-// ) const
-// {
-//   const typename internal::nested<Derived,2>::type nested(derived());
-//   const typename internal::nested<OtherDerived,2>::type otherNested(other.derived());
-//   return    (nested - otherNested).cwise().abs2().sum()
-//          <= prec * prec * (std::min)(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum());
-// }
+namespace Eigen {
+  
+template<typename Derived>
+template<typename OtherDerived>
+bool SparseMatrixBase<Derived>::isApprox(const SparseMatrixBase<OtherDerived>& other, const RealScalar &prec) const
+{
+  using std::min;
+  const typename internal::nested_eval<Derived,2,PlainObject>::type actualA(derived());
+  typename internal::conditional<IsRowMajor==OtherDerived::IsRowMajor,
+    const typename internal::nested_eval<OtherDerived,2,PlainObject>::type,
+    const PlainObject>::type actualB(other.derived());
+
+  return (actualA - actualB).squaredNorm() <= prec * prec * (min)(actualA.squaredNorm(), actualB.squaredNorm());
+}
+
+} // end namespace Eigen
 
 #endif // EIGEN_SPARSE_FUZZY_H
diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h
index b5c50d93a..a438c6487 100644
--- a/Eigen/src/SparseCore/SparseMatrixBase.h
+++ b/Eigen/src/SparseCore/SparseMatrixBase.h
@@ -324,8 +324,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
 
     template<typename OtherDerived>
     bool isApprox(const SparseMatrixBase<OtherDerived>& other,
-                  const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const
-    { return toDense().isApprox(other.toDense(),prec); }
+                  const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
 
     template<typename OtherDerived>
     bool isApprox(const MatrixBase<OtherDerived>& other,
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index c86534bad..e5b6d5a0a 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -304,6 +304,10 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
 
     VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint());
+    
+    // check isApprox handles opposite storage order
+    typename Transpose<SparseMatrixType>::PlainObject m3(m2);
+    VERIFY(m2.isApprox(m3));
   }