I am interested in implementing a resolvent (input/output) analysis framework using the FreeFEM/PETSc interface.

To do this, I need to construct a matrix that is defined by a series of products and of several sparse matrices and matrix inverses. The matrix of interest is Hermitian, with a dimension of n\times n. It is defined as L=B^HR^HM_qRB.

Here, B is a m\times n matrix of 1s and 0s with m\geq n, M_q is a positive semi-definite m\times m matrix, and R is an m\times m matrix defined by an inverse as R=(i\omega M_q+J)^{-1}.

With the exception of R, I can construct the Mat objects for each of the basic components of L, and I understand how to perform the necessary MatMatMult() operations. However, I do not know how to use MUMPS to find R. It seems like this would require the MatMatSolve(), but I couldnâ€™t find any documentation on this in FreeFEM. Has anyone encountered a similar problem and found a solution?

For context, I need to construct this matrix in order to solve the eigenvalue problem (using the SLEPc interface):

L\tilde{\mathbf{f}}=\lambda M_f\tilde{\mathbf{f}}, where M_f is a positive definite n\times n matrix.

A good overview with more detail on this subject is given here (see Section 4).

You cannot access R from MUMPS, nobody can, except MUMPS developers. So youâ€™ll need to solve (i\omega M_q+J)C = B using KSPSolve() in your .edp, by first converting B to a dense Mat since KSPMatSolve() (in PETSc library) only handles dense right-hand sides. This will be extremely costly, these inversions make your L dense, and depending on which \lambda you are looking for, SLEPc may have to â€śinvertâ€ť L, so maybe it would be best to think of an approximation of L^{-1} which would not need an explicit representation of your L.

You are correct that this would be very expensive. Thanks to your comment, I think I see the solution. Iâ€™ll state it briefly here in case it helps anyone in the future.

Define a func which computes the action of the operator L on \tilde{\mathbf{f}} i.e. \mathbf{x}=L\tilde{\mathbf{f}}.

Use this func to create the Arnoldi basis in the EVP in a matrix-free approach (Ă la FreeFem-tutorial - Section 8 - example11.edp). This way the LU decomposition is only factored once and applied only as many times as is needed to build the Arnoldi subspace.

That would indeed be orders of magnitude cheaper to compute.
This is like for PCD preconditioning, cf. FreeFem-sources/oseen-2d-PETSc.edp at develop Â· FreeFem/FreeFem-sources Â· GitHub, the PCD operator is never assembled, but is computed instead by a sequence of MatMult() or KSPSolve(). In your case, youâ€™ll need to implement:

MatMult() in a func as a sequence of MatMult(B, ...) + KSPSolve(R, ...) + MatMult(Mq, ...) + KSPSolveHermitianTranspose(R, ...) + MatMultHermitianTranspose(B, ...)

PCApply() â€“ for L^{-1} if you are doing shift-and-invert in SLEPc for smallest \lambda â€“ in a func (maybe an initial approximation could be simply M_q^{-1})

LU decomposition of what? L? You wonâ€™t have that. R or M_q? Yes, whatever you choose as a preconditioner for these two operators, PETSc is smart enough to reuse them throughout the Arnoldi iterations.

Thank you for your suggestions. I believe I am close to having this working, but I am a bit confused on the required syntax.

I have defined a func complex[int] Lop(complex[int]& inPETSC) {...} to handle the sequence of MatMult()s and KSPsolve()s that define the operator L. I have tested this function independently (i.e. outside of the EPSSolve()) and can confirm that it works as intended.

My struggle right now has to do with how to get EPSSolve() to do what I want. In particular, I donâ€™t understand how to use the precon argument used in your examples. My current code is similar to what Iâ€™ve given below. It runs without error, but the results Iâ€™m getting donâ€™t make sense yet. This could be due to BCâ€™s or something else, but I wondered if you see an obvious error based on my lack of experience with PETSc.

Mat<complex> L(Mf, Lop); // Mf is a Mat<complex>, Lop is a func
int nev = getARGV("-nev",1);
string EPSparams = " -eps_nev " + nev + " " +
" -eps_type krylovschur " +
" -eps_largest_real " +
" -st_pc_type none " +
" -eps_gen_hermitian ";
complex[int] val(nev);
Xh<complex>[int] deff(vec)(nev); // Xh is a fespace
int k = EPSSolve(L, Mf, vectors = vec, values = val, sparams = EPSparams);

I believe you donâ€™t need -st_pc_type none with -eps_largest_real. What do you mean by â€śthe results [â€¦] donâ€™t make senseâ€ť? Are you getting a meaningful -eps_view? What about -eps_monitor and -eps_converged_reason?