I am using the Eigen library in C++ for solving sparse linear equations: Ax=b where, A is a square sparse matrix and b is a dense vector with ILU-preconditioned BiCGSTAB. I am initializing the matrix A using the setFromTriplets function. The linear system is generated from discretization of partial differential equations in space and time.

My application changes the matrix slightly at every time-step. I want to modify a small number of rows (around 1% rows) in the matrix in the beginning of each time-step. I am storing the matrix in the row-major format so that I can access the row directly. I don't want to reassemble the entire matrix from triplets since the number of rows to be modified are around 1%. Moreover, the modification is such that the number of nonzeros in the row are exactly identical. I just want to change the column indices and values. Hence, I do not need to allocate extra memory for the row. After going through the Eigen documentation, I found the functions coeffRef and insert. Both of them will allocate extra memory if the element does not exist. I would like to avoid this since the number of nonzeros are not changing.

Any help is appreciated.

Let's say your sparse matrix declaration looks like this:

Code: Select all

Eigen::SparseMatrix<double, Eigen::RowMajor> A


Internally, eigen stores this matrix in Compressed Sparse Row (CSR) format.
An m by n CSR matrix with nnz nonzero entries has vectors, VALUE and COLUMN, both length nnz.
The rows are represented by a vector, ROWPTR, of length m + 1.
If you want all elements in row i, then you look at all elements in VALUE and COLUMN starting inclusively at ROWPTR[i] and ending exclusively at RWPTR[i + 1].

With this CSR datatype, you can get

Code: Select all

const double * ROWPTR = A.outerIndexPtr();
const double * COLUMN = A.innerIndexPtr();
const double * VALUE = A.valuePtr();


As an example of a column-value modification (assuming there aren't any elements explicitly defined on the last column),

Code: Select all

    const int Inclusive_Row_5_Begin = ROWPTR[5];
    const int Exclusive_Row_5_End = ROWPTR[5 + 1];

    for(int i = Inclusive_Row_5_Begin; i < Exclusive_Row_5_End; ++i)
    {
        VALUE[i] = 2;
        COLUMN[i] += 1;
    }