Hi,

is it possible to share the matrix structure (not the data) between multiple matrices?

Application: In a finite element code I have multiple sparse matrices, which however have exactly the same sparsity pattern. To save memory I would like to generate one 'master' matrix (sparsity structure + data) and have all other matrices share the structure with it.

Any help is highly appreciated.

Matthias

This can be accomplished with some manual effort via Map<SparseMatrix<double >:

https://eigen.tuxfamily.org/dox-devel/c ... _01_4.html

So the idea is to create one "reference" SparseMatrix for the sparsity pattern, and as many non-zeros vectors as you need (e.g., using VectorXd), and "assemble" them Map<SparseMatrix<double > and the SparseMatrix accessors like outerIndexPtr(), innerIndexPtr().

Make sure the SparseMatrix is is compressed mode by calling makeCompress() once.

Hi,

Thank you very much for the information on sharing the structure between multiple matrices and assembling them.

As a follow-up question please consider the following expression

res = A*f(u,v,w) + B*g(u,v,w) + C*h(u,v,w)

which corresponds for instance to the discretized divergence operator applied to the flux functions F=([f,g,h](u,v,w). All three matrices A, B, and C have the same sparsity structure derived from the finite element connectivity pattern.

Is it possible, ideally without implementing the sparse matrix-vector multiplication from scratch, to exploit the fact of a shared matrix structure and to prevent that temporary objects are created and the vectors u, v, and w are fetched multiple times from memory? To be more precise, is it possible to instruct Eigen that the three sparse matrix-vector multiplications will lead to exactly the same for-loops and can therefore be fused into a single pass thriugh the matrix? As an alternative it would also be possible to store both the matrix cofficients and the vector data as multiblock matrix/vector with a slightly adjusted matrix-vector multiplication. Do multiblock matrices/vectors exist in Eigen?

Any further help is appreciated.

Kind regards
Matthias

One way to go is to pack A, B, and C in a SparseMatrix<Array3d,RowMajor>, pack u, v, w in a Matrix<Array3d,Dynamic,1> and glue f, g, and h in a single function from Array3d to Array3d.

I proposed a RowMajor sparse matrix to get advantage of OpenMP if you compile with -fopenmp (or the equivalent compiler switch).

Thanks very much for the quick response.

Just a final question. Will the sparse matrix-vector product of the block version you propose result in a block vector or will it take care of the internal summation over the three components, that is, will it yield a scale sparse matrix?

Good point, so you probably need to wrap the Array3d into a small custom type with operator* overloaded to return a double, and then specialize Eigen::ScalarBinaryOpTraits<YourWrapper,YourWrapper>:

https://eigen.tuxfamily.org/dox-devel/s ... raits.html

to let Eigen knows.