# HG changeset patch
# User Joshua N. Pritikin <jpritikin@pobox.com>
# Date 1449686307 18000
#      Wed Dec 09 13:38:27 2015 -0500
# Node ID 0d4145475a392efdb0b3664b266c789bd6f9831c
# Parent  19ad4aee60c82f931fbd9b9f4f4e40b56acac708
Add sparse, in-place filterOuter

diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h
--- a/Eigen/src/SparseCore/SparseMatrix.h
+++ b/Eigen/src/SparseCore/SparseMatrix.h
@@ -493,16 +493,46 @@ class SparseMatrix
         return; 
       m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
       for (Index i = 0; i < m_outerSize; i++)
       {
         m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; 
       }
     }
     
+    /** Filter outer dimension of matrix */
+    template <typename T>
+    void filterOuter(std::vector<T> &toInclude) {
+      eigen_assert(Index(toInclude.size()) == outerSize() &&
+		   "filter mask must be the same size as the outer dimension");
+      uncompress();
+      Index dx = 0;
+      for (Index cx=0; cx < outerSize(); ++cx) {
+	if (!toInclude[cx]) continue;
+	m_outerIndex[dx] = m_outerIndex[cx];
+	m_innerNonZeros[dx] = m_innerNonZeros[cx];
+	++dx;
+      }
+      if (dx == outerSize()) return; // didn't filter anything
+      m_outerIndex[dx] = m_outerIndex[outerSize()];
+
+      StorageIndex *newOuterIndex = static_cast<StorageIndex*>(std::realloc(m_outerIndex, (dx + 1) * sizeof(StorageIndex)));
+      if (!newOuterIndex) internal::throw_std_bad_alloc();
+      m_outerIndex = newOuterIndex;
+      if (dx == 0) {
+	std::free(m_innerNonZeros);
+	m_innerNonZeros = 0;
+      } else {
+	StorageIndex *newInnerNonZeros = static_cast<StorageIndex*>(std::realloc(m_innerNonZeros, dx * sizeof(StorageIndex)));
+	if (!newInnerNonZeros) internal::throw_std_bad_alloc();
+	m_innerNonZeros = newInnerNonZeros;
+      }
+      m_outerSize = dx;
+    }
+
     /** Suppresses all nonzeros which are \b much \b smaller \b than \a reference under the tolerence \a epsilon */
     void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits<RealScalar>::dummy_precision())
     {
       prune(default_prunning_func(reference,epsilon));
     }
     
     /** Turns the matrix into compressed format, and suppresses all nonzeros which do not satisfy the predicate \a keep.
       * The functor type \a KeepFunc must implement the following function:
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -429,16 +429,59 @@ template<typename SparseMatrixType> void
           m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1;
           
         VERIFY_IS_APPROX(m1, refMat1);
           
           
       }
   }
 
+  // test filterOuter
+  {
+    DenseMatrix ref = DenseMatrix::Zero(rows, cols);
+    SparseMatrixType m1(rows, cols);
+    for(int k=0; k<std::max(rows*cols/4, Index(1000)); ++k)
+    {
+      Index i = internal::random<Index>(0,rows-1);
+      Index j = internal::random<Index>(0,cols-1);
+      Scalar v = internal::random<Scalar>();
+      m1.coeffRef(i,j) = v;
+      ref.coeffRef(i,j) = v;
+    }
+    if(internal::random<Index>(0,10)<5) {
+      m1.makeCompressed();
+    }
+    std::vector<bool> filter(m1.outerSize());
+    Index retained = 0;
+    for (int fx=0; fx < m1.outerSize(); ++fx) {
+      bool inc = internal::random<bool>();
+      filter[fx] = inc;
+      retained += inc;
+    }
+    m1.filterOuter(filter);
+    VERIFY(m1.outerSize() == retained);
+    if (retained == (Index) filter.size()) {
+      VERIFY_IS_APPROX(m1, ref);
+    } else {
+      if (retained == m1.cols()) {
+	for (int fx=0, dx=0; fx < ref.cols(); ++fx) {
+	  if (!filter[fx]) continue;
+	  VERIFY_IS_APPROX(m1.col(dx), ref.col(fx));
+	  dx++;
+	}
+      } else {
+	for (int fx=0, dx=0; fx < ref.rows(); ++fx) {
+	  if (!filter[fx]) continue;
+	  VERIFY_IS_APPROX(m1.row(dx), ref.row(fx));
+	  dx++;
+	}
+      }
+    }
+  }
+
   // test Identity matrix
   {
     DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows);
     SparseMatrixType m1(rows, rows);
     m1.setIdentity();
     VERIFY_IS_APPROX(m1, refMat1);
     for(int k=0; k<rows*rows/4; ++k)
     {