> El 31 oct 2016, a las 22:40, Peetz, Darin T <[email protected]> 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
ex13.c
Description: Binary data
