[petsc-users] MUMPS icntl for petsc4py

[Jose E. Roman]

Here is the explanation. With shift-and-invert, two main things must be done at STSetUp: 1) build the matrix A-sigma*B, (2) factorize it. Normally this is done at the beginning of EPSSolve. Before that you can set PC options, but the problem is that MUMPS options belong to Mat, not PC, so step 1) must be done beforehand. But you cannot call PCSetUp because you have not yet configured MUMPS options. Around version 3.12 we split the implementation of STSetUp so that 1) and 2) can be done separately. STGetOperator is what triggers 1).

Jose


> El 23 ene 2023, a las 13:13, Quentin Chevalier <quentin.chevalier at polytechnique.edu> escribió:
> 
> Many thanks Jose, it works beautifully !
> 
> I'm at a loss as to why, but thanks for the quick fix !
> 
> Quentin
> 
> 
> 
> Quentin CHEVALIER – IA parcours recherche
> 
> LadHyX - Ecole polytechnique
> 
> __________
> 
> 
> 
> On Sun, 22 Jan 2023 at 10:58, Jose E. Roman <jroman at dsic.upv.es> wrote:
>> 
>> You have to call ST.getOperator() as is done in this C example:
>> https://slepc.upv.es/documentation/current/src/eps/tutorials/ex43.c.html
>> 
>> Jose
>> 
>> 
>>> El 22 ene 2023, a las 10:40, Quentin Chevalier <quentin.chevalier at polytechnique.edu> escribió:
>>> 
>>> Hello PETSc users,
>>> 
>>> I'm getting an INFOG(1)=-9 and INFO(2)=27 error on an eigenvalue code based on dolfinx run in a docker container. Based on https://mumps-solver.org/doc/userguide_5.5.1.pdf, I figured the fix would be to increase ICNTL(14).
>>> 
>>> I'm coding in python through the petsc4py/slepc4py wrapper. I found a Mat.setMumpsIcntl method but I can't seem to place it properly and always obtain another error : "Operation done in wrong order and the like".
>>> 
>>> Here's the code snippet that is failing :
>>> # Solver
>>> EPS = SLEPc.EPS().create(COMM_WORLD)
>>> EPS.setOperators(-A,M) # Solve Ax=sigma*Mx
>>> EPS.setProblemType(SLEPc.EPS.ProblemType.PGNHEP) # Specify that A is not hermitian, but M is semi-definite
>>> EPS.setWhichEigenpairs(EPS.Which.TARGET_MAGNITUDE) # Find eigenvalues close to sigma
>>> EPS.setTarget(sigma)
>>> EPS.setDimensions(2,10) # Find k eigenvalues only with max number of Lanczos vectors
>>> EPS.setTolerances(1e-9,100) # Set absolute tolerance and number of iterations
>>> # Spectral transform
>>> ST = EPS.getST(); ST.setType('sinvert')
>>> # Krylov subspace
>>> KSP = ST.getKSP()
>>> KSP.setTolerances(rtol=1e-6, atol=1e-9, max_it=100)
>>> # Krylov subspace
>>> KSP.setType('preonly')
>>> # Preconditioner
>>> PC = KSP.getPC(); PC.setType('lu')
>>> PC.setFactorSolverType('mumps')
>>> KSP.setFromOptions()
>>> EPS.setFromOptions()
>>> PC.getFactorMatrix().setMumpsIcntl(14,50)
>>> print(f"Solver launch for sig={sigma:.1f}...",flush=True)
>>> EPS.solve()
>>> n=EPS.getConverged()
>>> 
>>> For context, matrix A is complex, size 500k x 500k but AIJ sparse, and I'm running this code on 36 nodes.
>>> 
>>> I'd appreciate any insight on how to fix this issue, it's not clear to me what the order of operations should be. Funnily enough, it's very shift-dependent.
>>> 
>>> Cheers,
>>> 
>>> Quentin
>>