I am aware that in order to efficiently fill sparse matrices it is necessary to reserve room for the non-zeros per column (column major matrix), as one can appreciate in http://eigen.tuxfamily.org/dox/group__TutorialSparse.html.

Although I was able to come up with a good number for the problem sizes I was running up until now (see the previous post https://forum.kde.org/viewtopic.php?f=74&t=128919), as I increase the size of the problem - but not the sparsity pattern - the matrix filling has become extremely slow.

It is sort of a dummy question, however I wouldn't expect the behaviour I am experiencing now. Although I have tried different strategies to come up with a suitable reserve size, I am sure I am missing something...

Because the sparsity of the problem I am running is quite standard - an hepta-diagonal matrix arising from the discretization of a 3D PDE - I wonder if someone could help me with a good calculation that gives me a good reserve size, perhaps based on the the size of the problem (i.e. Nx x Ny x Nz).

Thanks in advance!

What is your current strategy?

The following should do the job:
mat.reserve(VectorXi::Constant(n,7)));
If not, then this means your matrix is not hepta-diagonal.

Thank you very much for the reply, ggael!

My current strategy is to over estimate the reserve size, basically because I am solving 1D, 2D, and 3D problems that result in, respectively, tri-, penta-, and hepta-diagonal matrices. Of course it is a quite naive approach as I should be able to come up with a more precise reserve size by looking at the dimensionality of the problem. However this strategy was working just fine for the size of problems I was solving up until now (aprox. 15000 unknows). But now that I am going higher number of unknows things became really slow. Below I share a code snippet that illustrates the type of slowness I am experiencing:

Code: Select all

            size_t ni = 1000;
            size_t nj = 1000;
            size_t nk = 1;
            size_t n = ni*nj*nk;
            size_t elementsPerRow = 1;

            Eigen::SparseMatrix<double> matrix(n, n);
            matrix.reserve(Eigen::VectorXi::Constant(n, elementsPerRow));

            for (size_t i = 0; i < n; ++i)
            {
                matrix.coeffRef(i, i) = 1.0;
            }


In the code above I am simply filling the main diagonal of the matrix. So, according to ggael suggestion, I believe the reserve size should be just 1. However, it takes so much time to run that I am not even able to report how much time it would take to fill this matrix up.

I am sure I am still missing something very basic...

Here, this code takes 0.01s using either Eigen 3.2 or 3.3. Don't forget the -O3 flag or equivalent for MSVC.

Hmmm... It now seems that I oversimplified my test case by keeping out too much from the original code.

Although I found that I was not compiling using the proper compilation flags - something I thought would also have a major influence in the performance of the code in the smaller tests as well - I fixed that and the code runs much faster. Thanks for the suggestion.

However, I was able to find another "hot-spot" by profiling my code and it has to do with my actual algorithm to fill the matrix up. Apparently too many cache misses.... Not sure at this point... I will try out some other things and post my findings here.

I turns out that the major impact on the code performance was the logic to fill the matrix... which became obvious after a clean up on my compilation flags.