# Resolvent operator with FF/PETSc (MatMatSolve?)

**URL:** <https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603>\
**Category:** General Discussion\
**Created:** [March 16, 2022, 4:05pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603 "2022-03-16T16:05:45Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [March 16, 2022, 4:05pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/1 "2022-03-16T16:05:45Z")

</div>

Hello FF developers,

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](https://hal.archives-ouvertes.fr/hal-00756811/document) (see Section 4).

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [March 17, 2022, 7:22am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/2 "2022-03-17T07:22:47Z")

</div>

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.

See [FreeFem-tutorial - Section 8 - example14.edp](http://joliv.et/FreeFem-tutorial/section_8/example14.edp.html) for an example of solves with multiple right-hand sides.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [March 17, 2022, 10:21am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/3 "2022-03-17T10:21:24Z")

</div>

@prj thank you for the quick reply.

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.

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

2. Use this `func` to create the Arnoldi basis in the EVP in a matrix-free approach (à la [FreeFem-tutorial - Section 8 - example11.edp](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/section_8/example11.edp.html)). This way the LU decomposition is only factored once and applied only as many times as is needed to build the Arnoldi subspace.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [March 17, 2022, 11:14am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/4 "2022-03-17T11:14:45Z")

</div>

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](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/oseen-2d-PETSc.edp#L110), 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})

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [March 17, 2022, 11:22am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/5 "2022-03-17T11:22:18Z")

</div>

> [@cmd](#):
>
> This way the LU decomposition

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.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [March 24, 2022, 6:03pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/6 "2022-03-24T18:03:22Z")

</div>

Thanks, yes I was referring the the LU decomposition of R. I can see that PETSc is automatically reusing that factorization.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [March 24, 2022, 6:56pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/7 "2022-03-24T18:56:18Z")

</div>

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.

```auto
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);

```

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [March 24, 2022, 8:04pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/8 "2022-03-24T20:04:07Z")

</div>

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`?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [March 24, 2022, 9:41pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/9 "2022-03-24T21:41:52Z")

</div>

Yes, these are all meaningful. I think the issue is with my varfs/BCs. Thank you for your helpful insights as always @prj 🙂

---

<div class="post-metadata">

**Author:** ![tea](https://avatars.discourse-cdn.com/v4/letter/t/7ea924/32.png) [@tea](https://community.freefem.org/u/tea)\
**Post date:** [December 15, 2023, 9:17am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/10 "2023-12-15T09:17:53Z")

</div>

Hello Chris,

could you maybe share the final code?

Regards,  
Tea

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [December 15, 2023, 12:03pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/11 "2023-12-15T12:03:32Z")

</div>

Hi there, yes. It is available via Github. What is your username?

---

<div class="post-metadata">

**Author:** ![tea](https://avatars.discourse-cdn.com/v4/letter/t/7ea924/32.png) [@tea](https://community.freefem.org/u/tea)\
**Post date:** [December 15, 2023, 1:20pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/12 "2023-12-15T13:20:09Z")

</div>

Hi Chris,

Thanks for the response!  
You mean my Github username? It’s TeaVoj. You can also share the link.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [December 16, 2023, 12:07am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/13 "2023-12-16T00:07:25Z")

</div>

Hi Tea,

I will let you know when the repo goes public. For now, the relevant portion of the code is below.

Cheers!

```auto
// construct matrices
M = vMq(XMh, XMh); // Response Norm
Mf = vMf(Xh, Xh); // Forcing Norm
matrix<complex> LocPQ = vP(Xh, XMh); // Forcing/Response Correspondence
Mat<complex> PQ(M, Mf, LocPQ);

func complex[int] LHSop(complex[int]& inPETSc) {
  complex[int] temp(XMh.ndof), outPETSc(inPETSc.n);
  MatMult(PQ, inPETSc, outPETSc);
  KSPSolve(J, outPETSc, temp);
  MatMult(M, temp, outPETSc);
  KSPSolveHermitianTranspose(J, outPETSc, temp);
  MatMultHermitianTranspose(PQ, temp, outPETSc);
  return outPETSc;
}

Mat<complex> LHS(Mf, LHSop);

J = vJ(XMh, XMh, tgv = -1); //Linear operator
set(J, sparams = KSPparams);
int k = EPSSolve(LHS, Mf, vectors = fvec, values = val,
                 sparams = "-eps_type krylovschur -eps_largest_real -eps_monitor_conv -options_left no -eps_gen_hermitian");

```

---

<div class="post-metadata">

**Author:** ![tea](https://avatars.discourse-cdn.com/v4/letter/t/7ea924/32.png) [@tea](https://community.freefem.org/u/tea)\
**Post date:** [December 19, 2023, 3:12pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/14 "2023-12-19T15:12:44Z")

</div>

Hi Chris,

thanks a lot. This is very helpful.

I am still having some issues with my code, it is probably due to my inexperience with FreeFem.  
Could you maybe share which KSPparams you use?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [December 19, 2023, 3:24pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/15 "2023-12-19T15:24:35Z")

</div>

Glad that helps. A good rule of thumb is to always start with: `"-ksp_type preonly -pc_type lu"`

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [December 21, 2023, 11:06am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/16 "2023-12-21T11:06:00Z")

</div>

Dear Chris,

I am considering using the resolvent framework, and your code is very helpful. Could you share the variational formulation of `vMq`, `vMf` and `LocPQ`, and the corresponding FE spaces?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [December 21, 2023, 11:30am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/17 "2023-12-21T11:30:36Z")

</div>

Hi @aszaboa and @tea, I’ve opened the [ff-bifbox repo](https://github.com/cmd8/ff-bifbox) to the public now. Note that the code is still under development and some things may not work perfectly yet.

You can see/test my resolvent analysis implementation by following the commands in `examples/garnaud_2012`.

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [December 22, 2023, 9:10am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/18 "2023-12-22T09:10:03Z")

</div>

Thank you so much for sharing the code!

---

<div class="post-metadata">

**Author:** ![clthu](https://avatars.discourse-cdn.com/v4/letter/c/f04885/32.png) [@clthu](https://community.freefem.org/u/clthu)\
**Post date:** [August 9, 2024, 2:02am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/19 "2024-08-09T02:02:43Z")

</div>

Hi, Chris.

I have some small puzzle about the rslvcompute.edp in your ff-bifbox repo. First, the tgv = -2 for linear operator and tgv = -20 for response norm. Is this better than tgv = -1 for linear operator and no tgv for response norm ( I used to apply ARPACK to solve the problem with this way to construct matrix in FreeFem++ and don’t encounter problems)? Second, to compute the response modes, I think the complex array gm or qm is in PETSc indexing, and its length should be sized as qm(J.n) instead of qm(XMh.n).

Thanks very much for your codes and I learn much from it!

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 9, 2024, 2:26am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/20 "2024-08-09T02:26:53Z")

</div>

Hi @clthu,

Thanks for the message. I’m glad you found the codes helpful! To answer your post:

1. Yes, `tgv=-2` is needed here instead of `tgv=-1`, because we are solving both with `KSPSolve()` (i.e. needs row elimination to enforce BCs) and with `KSPSolveHermitianTranspose()` (i.e. needs column elimination to enforce BCs). Another option would be to use the normal `tgv=1e30` approach, but I avoided that here to avoid potential issues with iterative solvers.
2. Yes, you are correct – [good catch](https://github.com/cmd8/ff-bifbox/commit/0280ccd68d7d9f4856aed8d840120610f453eba5). However, this does not cause problems in this case since `ChangeNumbering()` automatically resizes the arrays.

[Next page](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603.md?page=2)