[petsc-users] Provide Matrix Factorization to EPS for Generalized Eigenvalue Problem

[Peetz, Darin T]

Thank you for providing the example.  I have managed to get it working properly; however, after looking at the example, I have one more question.

First, let me provide a little background on how the code is set up.  I have one function that takes care of solving the linear system of equations, and a separate one that takes care of the eigenvalue problem.  The KSP object for the linear system is passed between both functions, and kept for the duration of the program to save some of the cost of factorizing the matrix each time.  I have been keeping the EPS object local to the function that takes care of the eigenvalue problem since it is only relevant within that function, and I haven't seen a need to reuse the same EPS object in every iteration (I'm assuming the cost of creating the EPS object itself is relatively small).  The reason for the extra symbolic factorizations seems to have been that when I destroyed the EPS object before it goes out of scope, it also destroys the underlying KSP object that I provided.

Basically, I'm wondering if my assumption about the cost of creating the EPS object is correct.  Is it okay to recreate the EPS object every time the eigenvalue function is called, or does this create a significant overhead?  It seems that if I try to reuse the same EPS object I have to make sure I'm reusing the same A-matrix as well, or else the ST object redoes both the symbolic and numeric factorizations (so the same K matrix is numerically factorized twice in each iteration) each time.  Thus I would be able to save a significant amount of memory if I could let the EPS object and A-matrix exist in a local scope.

As a side note, I worked around the KSP destruction issue by providing an empty (uninitialized) KSP object to the EPS object after calling EPSSolve and extracting the eigenvalues/eigenvectors.  This is functional but seems like a bad/crude approach.  Is there a more elegant way to "free" the KSP object from the EPS object once it is no longer needed, or should I stick with this approach?

Thanks a lot for your help,
Darin

________________________________________
From: Jose E. Roman [jroman at dsic.upv.es]
Sent: Tuesday, November 01, 2016 6:14 AM
To: Peetz, Darin T
Cc: petsc-users at mcs.anl.gov
Subject: Re: [petsc-users] Provide Matrix Factorization to EPS for Generalized Eigenvalue Problem

> El 31 oct 2016, a las 22:40, Peetz, Darin T <peetz2 at illinois.edu> escribió:
>
> Hello,
>
> I'm wondering how I could go about providing a matrix factorization calculated in Petsc to the eigenvalue routines in Slepc.  I'm trying to solve the eigenvalue problem for stability, where the solution to KU=F is needed to construct the A-matrix (K_sigma) in the eigenvalue problem.  Since the eigenvalue problem is generalized, it seems like the best way to solve it is to factorize the B-Matrix (K, same as in KU=F) with a package like MUMPS and use a method in Slepc such as Krylov-Schur.  Since I need to solve both KU=F and the eigenvalue problem, I'd like to compute the factorization of K first to solve KU=F, and then reuse it in the EPS routines.
>
> I've tried using EPSGetST() and STSetKSP() to provide the KSP object that I used to solve KU=F, but then for some reason when I change nonzero values in K (but not nonzero locations) Petsc redoes the symbolic factorization when I go to solve KU=F again (it's part of an optimization routine, so I'm solving both problems, updating parts of K, and repeating).  This does provide the correct solution, and allows me to use the same factorization for KU=F and the eigenvalue problem, but the extra symbolic factorizations, while comparatively cheap, are unnecessary and ideally should be eliminated.  If I skip the calls to EPSGetST() and STSetKSP(), the symbolic factorization for the KSP object associated with KU=F is only performed once, as it should be.  Is there some option I'm overlooking, or maybe a better way to go about this?
>
> Thanks,
> Darin

I don't know what your code is doing exactly, but I tried to reproduce the behaviour in a SLEPc example and it works as expected. The attached example calls KSPSolve()+EPSSolve()+KSPSolve()+EPSSolve(), modifying B in the middle, and it computes just one symbolic factorization and two numeric factorizations.

An alternative scheme would be the following: pass the matrices to EPS, call EPSSetUp(), then use EPSGetST()+STGetKSP() to extract the KSP object from EPS. This way you could solve linear systems without having to create the KSP object yourself.

A second alternative would be to solve a standard eigenproblem with a shell matrix instead of a generalized eigenproblem. That is, instead of EPSSetOperators(A,B) do EPSSetOperators(S,NULL), where S is a shell matrix whose MATOP_MULT operation does a MatMult(A) followed by a KSPSolve(B).

Jose