5

Suppose that I have a large sparse matrix with the following pattern:

  • the number of nonzeros per column and their locations are fixed
  • only matrix block A and B will change and the rest of the matrix stays static; (blocks A and B themselves are also sparse with fixed nonzero locations)

As instructed in the document, i've initialized the above matrix by

  • reserving the exact number of nonzeros per column for the column major sparse matrix
  • inserting column by column
  • inserting from the smallest row index per column

In later part of the program, it's natural to reuse the matrix and only updates the A, B blocks inplace. Possible ways are:

  1. accessing existing entries by coeffRef, would introduce binary search so not preferred here.
  2. iterating over the outer and inner dimensions as documented here

However, it seems a bit unnecessary to iterate over all nonzero entries since most part of the sparse matrix stays the same.

Is it possible to update A, B inplace without iterating over all nonzeros in the matrix?

wuyang3
  • 55
  • 7

1 Answers1

2

From what I can tell, the InnerIterator can be used used for this and runs in constant time.

Eigen::Index col = 1;
Eigen::Index offset_in_col = 1;

using SparseMatrixD = Eigen::SparseMatrix<double>;
SparseMatrixD mat = ...;

SparseMatrixD::InnerIterator i =
      SparseMatrixD::InnerIterator(mat, col) + offset_in_col;
assert(i.row() == 1);
assert(i.col() == 1);
assert(i.value() == C);

This should access the value C. All you need to know is how many nonzero elements are per column (or inner dimension in general). You don't need to know how many nonzero columns (outer dimensions) are stored because that array (SparseMatrix.outerIndexPtr()) has one entry per column.

Homer512
  • 9,144
  • 2
  • 8
  • 25
  • it actually worked, thanks! In the end i find that this is not the bottleneck for my application, this might be useful if the unchanged blocks are really large and with lots of nonzeros. btw, adding `Eigen::Index` to `InnerIterator` won't compile, i'm directly incrementing the iterator. – wuyang3 Jun 15 '22 at 05:09