How to get sparse complex matrices from my code to PETSc efficiently

Question

What is the most efficient way to get a complex sparse matrix from my Fortran code to PETSc? I understand that this is problem dependent, so I tried to give as many relevant details as possible below.

I've been playing with the FEAST eigenvalue solver [1] for problems of the type $A x = \lambda B x$, the dimension of the matrices $A$ and $B$ is $N$, and pretty much all the time is spent solving $N \times N$ complex linear system with M0 right hand sides. N is large (number of FE basis functions, in 3D), M0 is small (in my case I am interested in M0 ~ 20). The matrices $A$ and $B$ are real, symmetric and sparse, and the complex problem that needs solving is $zA-B$, where $z$ is a complex number. The author of FEAST seems to suggest, that the accuracy of the solution to this linear system doesn't have to be very high in order to get high accurate eigenvalues and eigenvectors, so some fast iterative solvers might be a great solution to this.

So far I've been just using Lapack for the complex system, and that works great for $N < 1500$ on my computer. For larger $N$, I don't know yet what the best solver is, so I wanted to just use PETSc and play with the iterative solvers there.

I wrote a simple C driver, and call it from Fortran, see [2] for all the code, but the problem is just with this part (update: I've put here all the lines to create the matrix, as I just realized that this is relevant):

ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
for (i=0; i<n; i++) col[i] = i;
ierr = MatSetValues(A,n,col,n,col,A_,INSERT_VALUES);CHKERRQ(ierr);

Which is extremely slow (i.e. for N~1500, this takes maybe 2s, while actually solving it is immediate), in fact the MatSetValues line takes pretty much 100% of all the time for the whole eigenvalue calculation... The matrix A_ is a 2D matrix that comes from Fortran. I tried to disable the MAT_IGNORE_ZERO_ENTRIES but it didn't make any difference. So I think that the problem simply is that even for moderate N like 1500, I need to use some sparse matrix format, is that correct?

So I implemented the CSR format in my Fortran code for the $A$ and $B$ matrices (or for the $zA-B$) and now I am trying to figure out how to efficiently give it to PETSc. For now, I just want to get something working sequentially, that beats Lapack. Should I use the MatCreateSeqAIJWithArrays function for this?

Is this the most efficient way to do that? Since the matrices $A$ and $B$ don't change, only the complex number $z$ changes in the FEAST algorithm, and in an FE calculation, I think that both $A$ and $B$ have the same sparse structure, one can probably improve things further by preallocating the sparse structure and then just quickly evaluate $zA_x - B_x$ in each FEAST iteration (the $A_x$, $B_x$ are the values arrays of the CSR format), I can do this easily in Fortran, but maybe it's slow to always call MatCreateSeqAIJWithArrays, so it might be faster to do all this in PETSc once the matrices $A$ and $B$ are transferred over.

I would like to know whether this CSR approach to this problem is correct, or whether I am doing it wrong (clearly my original approach with the MatSetValues is not optimal). Thanks for any tips.

[1] http://www.ecs.umass.edu/~polizzi/feast/

[2] https://github.com/certik/hfsolver/pull/14

[3] http://www.mcs.anl.gov/petsc/petsc-3.1/docs/manualpages/Mat/MatCreateSeqAIJWithArrays.html

Answer 1

It is important to preallocate correctly. This is almost certainly the reason why your assembly was slow. If you are starting with a dense matrix representation, just scan through it once counting the number of nonzeros per row, then call MatSeqAIJSetPreallocation(). See this FAQ. The option MAT_IGNORE_ZERO_ENTRIES is really intended to be used when there is some mild sparsity in those otherwise dense blocks rather than building the entire matrix in one call from a dense matrix. For this reason, it doesn't automatically do preallocation based on the sparsity of that one block.

Creating a dense intermediate matrix is not memory scalable, so you will eventually want to avoid it. MatSetValues() is really meant to be used for logically dense blocks in a sparse matrix. You typically either call once per row (or group of rows, typical of FD methods) or once per element (typical of FEM methods). If you are translating an existing assembled sparse matrix, just call MatSetValues() once per row. If you would prefer to skip the intermediate matrix (better performance and lower memory), just call MatSetValues() once per element.

Note that you can call PETSc directly from Fortran, though there are users of every Fortran dialect between basic Fortran 77 and the most recent versions. The interfaces will look rather crude to you as an eager adopter of the latest features. Suggestions for maintainable ways to improve support for the latest dialects would be appreciated.