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 = 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


1 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.

  • $\begingroup$ Thanks for the excellent answer. I can see now the idea behind MatSetValues(). For translating an existing matrix, isn't it faster to just call MatCreateSeqAIJWithArrays once rather than calling MatSeqAIJSetPreallocation and then MatSetValues for each row? $\endgroup$ Oct 24, 2012 at 3:06
  • $\begingroup$ I will provide feedback about the Fortran wrappers. For now it is the fastest for me to stick to writing the driver that I need in C, and only wrap this driver myself. $\endgroup$ Oct 24, 2012 at 3:10
  • 1
    $\begingroup$ Sure, but that requires you to start with a CSR matrix assembled using the same convention. If you use some other format, just pack one row at a time. If you preallocate, the time spent in MatSetValues() will be very small compared to computing the entries and solving the system. Also, assembling using MatSetValues[Blocked][Local]() prevents your code from depending on a particular matrix format, allowing you to choose storage formats at run-time. $\endgroup$
    – Jed Brown
    Oct 24, 2012 at 3:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.