diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 7f2884769..855a5a648 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -368,7 +368,7 @@ class SparseMatrix m_innerNonZeros[j] = innerNNZ; } if(m_outerSize>0) - m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize-1] + m_innerNonZeros[m_outerSize-1] + reserveSizes[m_outerSize-1]; + m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize-1] + m_innerNonZeros[m_outerSize-1] + internal::convert_index(reserveSizes[m_outerSize-1]); m_data.resize(m_outerIndex[m_outerSize]); } @@ -472,7 +472,65 @@ class SparseMatrix } } - //--- + // remove outer vectors j, j+1 ... j+num-1 and resize the matrix + void removeOuterVectors(Index j, Index num = 1) { + eigen_assert(num >= 0 && j >= 0 && j + num <= m_outerSize && "Invalid parameters"); + + const Index newRows = IsRowMajor ? m_outerSize - num : rows(); + const Index newCols = IsRowMajor ? cols() : m_outerSize - num; + + const Index begin = j + num; + const Index end = m_outerSize; + const Index target = j; + + // if the removed vectors are not empty, uncompress the matrix + if (m_outerIndex[j + num] > m_outerIndex[j]) uncompress(); + + // shift m_outerIndex and m_innerNonZeros [num] to the left + internal::smart_memmove(m_outerIndex + begin, m_outerIndex + end + 1, m_outerIndex + target); + if (!isCompressed()) + internal::smart_memmove(m_innerNonZeros + begin, m_innerNonZeros + end, m_innerNonZeros + target); + + // if m_outerIndex[0] > 0, shift the data within the first vector while it is easy to do so + if (m_outerIndex[0] > StorageIndex(0)) { + uncompress(); + const Index from = internal::convert_index(m_outerIndex[0]); + const Index to = Index(0); + const Index chunkSize = internal::convert_index(m_innerNonZeros[0]); + m_data.moveChunk(from, to, chunkSize); + m_outerIndex[0] = StorageIndex(0); + } + + // truncate the matrix to the smaller size + conservativeResize(newRows, newCols); + } + + // insert empty outer vectors at indices j, j+1 ... j+num-1 and resize the matrix + void insertEmptyOuterVectors(Index j, Index num = 1) { + EIGEN_USING_STD(fill_n); + eigen_assert(num >= 0 && j >= 0 && j < m_outerSize && "Invalid parameters"); + + const Index newRows = IsRowMajor ? m_outerSize + num : rows(); + const Index newCols = IsRowMajor ? cols() : m_outerSize + num; + + const Index begin = j; + const Index end = m_outerSize; + const Index target = j + num; + + // expand the matrix to the larger size + conservativeResize(newRows, newCols); + + // shift m_outerIndex and m_innerNonZeros [num] to the right + internal::smart_memmove(m_outerIndex + begin, m_outerIndex + end + 1, m_outerIndex + target); + // m_outerIndex[begin] == m_outerIndex[target], set all indices in this range to same value + fill_n(m_outerIndex + begin, num, m_outerIndex[begin]); + + if (!isCompressed()) { + internal::smart_memmove(m_innerNonZeros + begin, m_innerNonZeros + end, m_innerNonZeros + target); + // set the nonzeros of the newly inserted vectors to 0 + fill_n(m_innerNonZeros + begin, num, StorageIndex(0)); + } + } template void setFromTriplets(const InputIterators& begin, const InputIterators& end); diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index ddde617db..9e102b742 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -38,6 +38,7 @@ template void sparse_basic(const SparseMatrixType& re double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix DenseMatrix; typedef Matrix DenseVector; + typedef Matrix CompatibleDenseMatrix; Scalar eps = 1e-6; Scalar s1 = internal::random(); @@ -162,6 +163,74 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2,m1); } + // test removeOuterVectors / insertEmptyOuterVectors + { + for (int mode = 0; mode < 4; mode++) { + CompatibleDenseMatrix m1(rows, cols); + m1.setZero(); + SparseMatrixType m2(rows, cols); + Vector reserveSizes(outer); + for (Index j = 0; j < outer; j++) reserveSizes(j) = internal::random(1, inner - 1); + m2.reserve(reserveSizes); + for (Index j = 0; j < outer; j++) { + Index i = internal::random(0, inner - 1); + Scalar val = internal::random(); + m1.coeffRefByOuterInner(j, i) = val; + m2.insertByOuterInner(j, i) = val; + } + if (mode % 2 == 0) m2.makeCompressed(); + + if (mode < 2) { + Index num = internal::random(0, outer - 1); + Index start = internal::random(0, outer - num); + + Index newRows = SparseMatrixType::IsRowMajor ? rows - num : rows; + Index newCols = SparseMatrixType::IsRowMajor ? cols : cols - num; + + CompatibleDenseMatrix m3(newRows, newCols); + m3.setConstant(Scalar(NumTraits::quiet_NaN())); + + if (SparseMatrixType::IsRowMajor) { + m3.topRows(start) = m1.topRows(start); + m3.bottomRows(newRows - start) = m1.bottomRows(newRows - start); + } else { + m3.leftCols(start) = m1.leftCols(start); + m3.rightCols(newCols - start) = m1.rightCols(newCols - start); + } + + SparseMatrixType m4 = m2; + m4.removeOuterVectors(start, num); + + VERIFY_IS_CWISE_EQUAL(m3, m4.toDense()); + } else { + Index num = internal::random(0, outer - 1); + Index start = internal::random(0, outer - 1); + + Index newRows = SparseMatrixType::IsRowMajor ? rows + num : rows; + Index newCols = SparseMatrixType::IsRowMajor ? cols : cols + num; + + CompatibleDenseMatrix m3(newRows, newCols); + m3.setConstant(Scalar(NumTraits::quiet_NaN())); + + if (SparseMatrixType::IsRowMajor) { + m3.topRows(start) = m1.topRows(start); + m3.middleRows(start, num).setZero(); + m3.bottomRows(rows - start) = m1.bottomRows(rows - start); + } else { + m3.leftCols(start) = m1.leftCols(start); + m3.middleCols(start, num).setZero(); + m3.rightCols(cols - start) = m1.rightCols(cols - start); + } + + SparseMatrixType m4 = m2; + m4.insertEmptyOuterVectors(start, num); + + VERIFY_IS_CWISE_EQUAL(m3, m4.toDense()); + } + } + } + + // test sort if (inner > 1) { bool StorageOrdersMatch = int(DenseMatrix::IsRowMajor) == int(SparseMatrixType::IsRowMajor);