Hi Tea,
I will let you know when the repo goes public. For now, the relevant portion of the code is below.
Cheers!
// 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");