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

[Jose E. Roman]

> 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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex13.c
Type: application/octet-stream
Size: 7671 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161101/382eb09d/attachment.obj>