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