[petsc-users] Mysterious error code 77
Matthew Knepley
knepley at gmail.com
Fri May 6 08:34:59 CDT 2022
On Fri, May 6, 2022 at 9:28 AM Quentin Chevalier <
quentin.chevalier at polytechnique.edu> wrote:
> Sorry for forgetting the list. Making two matrices was more of a
> precaution then a carefully thought strategy.
>
> It would seem the MWE as I provided it above (with a setDimensions to
> reduce calculation time) does work in sequential and gives the eigenvalue
> 118897.88711586884.
>
> I had a working MUMPS solver config for a similar problem so I simply put
> it there. This gives the following code.
>
> However, running it even in sequential gives a "MatSolverType mumps does
> not support matrix type python" error.
> petsc4py.PETSc.Error: error code 92
> [0] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:136
> [0] EPSSetUp() at /usr/local/slepc/src/eps/interface/epssetup.c:350
> [0] STSetUp() at
> /usr/local/slepc/src/sys/classes/st/interface/stsolve.c:582
> [0] STSetUp_Sinvert() at
> /usr/local/slepc/src/sys/classes/st/impls/sinvert/sinvert.c:123
> [0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:408
> [0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:1016
> [0] PCSetUp_LU() at /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c:82
> [0] MatGetFactor() at /usr/local/petsc/src/mat/interface/matrix.c:4779
> [0] See https://petsc.org/release/overview/linear_solve_table/ for
> possible LU and Cholesky solvers
> [0] MatSolverType mumps does not support matrix type python
>
You would need type AIJ for MUMPS.
Thanks,
Matt
> Did I do something wrong ? I'm a bit confused by your FAQ entry, I never
> had a problem with this configuration in parallel.
>
> from petsc4py import PETSc as pet
> from slepc4py import SLEPc as slp
> from mpi4py.MPI import COMM_WORLD
>
> dir="./sanity_check_mats/"
>
> # Global sizes
> m=980143; m_local=m
> n=904002; n_local=n
> if COMM_WORLD.size>1:
> m_local=COMM_WORLD.rank*(490363-489780)+489780
> n_local=COMM_WORLD.rank*(452259-451743)+451743
> dir+="par/"
> else:
> m_local=m
> n_local=n
> dir+="seq/"
>
> QE=pet.Mat().createAIJ([[m_local,m],[n_local,n]])
> L=pet.Mat().createAIJ([[m_local,m],[m_local,m]])
> Mf=pet.Mat().createAIJ([[n_local,n],[n_local,n]])
>
> viewerQE = pet.Viewer().createMPIIO(dir+"QE.dat", 'r', COMM_WORLD)
> QE.load(viewerQE)
> viewerL = pet.Viewer().createMPIIO(dir+"L.dat", 'r', COMM_WORLD)
> L.load(viewerL)
> viewerMf = pet.Viewer().createMPIIO(dir+"Mf.dat", 'r', COMM_WORLD)
> Mf.load(viewerMf)
>
> QE.assemble()
> L.assemble()
> Mf.assemble()
>
> KSPs = []
> # Useful solvers (here to put options for computing a smart R)
> for Mat in [L,L.hermitianTranspose()]:
> KSP = pet.KSP().create()
> KSP.setOperators(Mat)
> KSP.setFromOptions()
> KSPs.append(KSP)
> class LHS_class:
> def mult(self,A,x,y):
> w=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)
> z=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)
> QE.mult(x,w)
> KSPs[0].solve(w,z)
> KSPs[1].solve(z,w)
> QE.multTranspose(w,y)
>
> # Matrix free operator
> LHS=pet.Mat()
> LHS.create(comm=COMM_WORLD)
> LHS.setSizes([[n_local,n],[n_local,n]])
> LHS.setType(pet.Mat.Type.PYTHON)
> LHS.setPythonContext(LHS_class())
> LHS.setUp()
>
> x, y = LHS.createVecs()
> x.set(1)
> LHS.mult(x,y)
> print(y.norm())
>
> # Eigensolver
> EPS = slp.EPS(); EPS.create()
> EPS.setOperators(LHS,Mf) # Solve QE^T*L^-1H*L^-1*QE*f=sigma^2*Mf*f
> (cheaper than a proper SVD)
> EPS.setProblemType(slp.EPS.ProblemType.GHEP) # Specify that A is
> hermitian (by construction), but M is semi-positive
> EPS.setWhichEigenpairs(EPS.Which.LARGEST_MAGNITUDE) # Find largest
> eigenvalues
> # Spectral transform
> ST = EPS.getST(); ST.setType('sinvert')
> # Krylov subspace
> KSP = ST.getKSP(); KSP.setType('preonly')
> # Preconditioner
> PC = KSP.getPC(); PC.setType('lu')
> PC.setFactorSolverType('mumps')
> EPS.setDimensions(1,5)
> EPS.setFromOptions()
> EPS.solve()
> print(EPS.getEigenvalue(0))
>
> Quentin CHEVALIER – IA parcours recherche
> LadHyX - Ecole polytechnique
> __________
>
>
>
> On Fri, 6 May 2022 at 14:47, Jose E. Roman <jroman at dsic.upv.es> wrote:
> >
> > Please respond to the list.
> >
> > Why would you need different matrices for sequential and parallel? PETSc
> takes care of loading matrices from binary files in parallel codes. If I
> use the 'seq' matrices for both sequential and parallel runs it works. I
> get the eigenvalue 128662.745858
> >
> > Take into account that in parallel you need to use a parallel linear
> solver, e.g., configure PETSc with MUMPS. See the FAQ #10
> https://slepc.upv.es/documentation/faq.htm#faq10
> >
> > Jose
> >
> > > El 5 may 2022, a las 14:57, Quentin Chevalier <
> quentin.chevalier at polytechnique.edu> escribió:
> > >
> > > Thank you for your answer, Jose.
> > >
> > > It would appear my problem is slightly more complicated than it
> appears. With slight modification of the original MWE to account for
> functioning in serial or parallel, and computing the matrices in either
> sequential or parallel (apologies, it's a large file), and including a
> slightly modified version of your test case, I obtain two different results
> in serial and parallel (commands python3 MWE.py and mpirun -n 2 python3
> MWE.py).
> > >
> > > Serial gives my a finite norm (around 28) and parallel gives me two
> norms of 0 and a code 77.
> > >
> > > This is still a problem as I would really like my code to be parallel.
> I saved my matrices using :
> > > viewerQE = pet.Viewer().createMPIIO("QE.dat", 'w', COMM_WORLD)
> > > QE.view(viewerQE)
> > > viewerL = pet.Viewer().createMPIIO("L.dat", 'w', COMM_WORLD)
> > > L.view(viewerL)
> > > viewerMq = pet.Viewer().createMPIIO("Mq.dat", 'w', COMM_WORLD)
> > > Mq.view(viewerMq)
> > > viewerMf = pet.Viewer().createMPIIO("Mf.dat", 'w', COMM_WORLD)
> > > Mf.view(viewerMf)
> > >
> > > New MWE (also attached for convenience) :
> > > from petsc4py import PETSc as pet
> > > from slepc4py import SLEPc as slp
> > > from mpi4py.MPI import COMM_WORLD
> > >
> > > dir="./mats/"
> > >
> > > # Global sizes
> > > m=980143; m_local=m
> > > n=904002; n_local=n
> > > if COMM_WORLD.size>1:
> > > m_local=COMM_WORLD.rank*(490363-489780)+489780
> > > n_local=COMM_WORLD.rank*(452259-451743)+451743
> > > dir+="par/"
> > > else:
> > > m_local=m
> > > n_local=n
> > > dir+="seq/"
> > >
> > > QE=pet.Mat().createAIJ([[m_local,m],[n_local,n]])
> > > L=pet.Mat().createAIJ([[m_local,m],[m_local,m]])
> > > Mf=pet.Mat().createAIJ([[n_local,n],[n_local,n]])
> > >
> > > viewerQE = pet.Viewer().createMPIIO(dir+"QE.dat", 'r', COMM_WORLD)
> > > QE.load(viewerQE)
> > > viewerL = pet.Viewer().createMPIIO(dir+"L.dat", 'r', COMM_WORLD)
> > > L.load(viewerL)
> > > viewerMf = pet.Viewer().createMPIIO(dir+"Mf.dat", 'r', COMM_WORLD)
> > > Mf.load(viewerMf)
> > >
> > > QE.assemble()
> > > L.assemble()
> > > Mf.assemble()
> > >
> > > KSPs = []
> > > # Useful solvers (here to put options for computing a smart R)
> > > for Mat in [L,L.hermitianTranspose()]:
> > > KSP = pet.KSP().create()
> > > KSP.setOperators(Mat)
> > > KSP.setFromOptions()
> > > KSPs.append(KSP)
> > >
> > > class LHS_class:
> > > def mult(self,A,x,y):
> > > w=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)
> > > z=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)
> > > QE.mult(x,w)
> > > KSPs[0].solve(w,z)
> > > KSPs[1].solve(z,w)
> > > QE.multTranspose(w,y)
> > >
> > > # Matrix free operator
> > > LHS=pet.Mat()
> > > LHS.create(comm=COMM_WORLD)
> > > LHS.setSizes([[n_local,n],[n_local,n]])
> > > LHS.setType(pet.Mat.Type.PYTHON)
> > > LHS.setPythonContext(LHS_class())
> > > LHS.setUp()
> > >
> > > x, y = LHS.createVecs()
> > > x.set(1)
> > > LHS.mult(x,y)
> > > print(y.norm())
> > >
> > > # Eigensolver
> > > EPS = slp.EPS(); EPS.create()
> > > EPS.setOperators(LHS,Mf) # Solve QE^T*L^-1H*L^-1*QE*f=sigma^2*Mf*f
> (cheaper than a proper SVD)
> > > EPS.setProblemType(slp.EPS.ProblemType.GHEP) # Specify that A is
> hermitian (by construction), but M is semi-positive
> > > EPS.setWhichEigenpairs(EPS.Which.LARGEST_MAGNITUDE) # Find largest
> eigenvalues
> > > EPS.setFromOptions()
> > > EPS.solve()
> > >
> > >
> > > Quentin CHEVALIER – IA parcours recherche
> > > LadHyX - Ecole polytechnique
> > > __________
> > >
> > >
> > >
> > >
> > >
> > > Quentin CHEVALIER – IA parcours recherche
> > > LadHyX - Ecole polytechnique
> > >
> > > __________
> > >
> > >
> > >
> > > On Thu, 5 May 2022 at 12:05, Jose E. Roman <jroman at dsic.upv.es> wrote:
> > > Your operator is not well formed. If you do this:
> > >
> > > x, y = LHS.createVecs()
> > > x.set(1)
> > > LHS.mult(x,y)
> > > y.view()
> > >
> > > you will see that the output is all zeros. That is why SLEPc complains
> that "Initial vector is zero or belongs to the deflation space".
> > >
> > > Jose
> > >
> > >
> > > > El 5 may 2022, a las 10:46, Quentin Chevalier <
> quentin.chevalier at polytechnique.edu> escribió:
> > > >
> > > > Just a quick amend on the previous statement ; the problem arises in
> > > > sequential and parallel. The MWE as is is provided for the parallel
> > > > case, but imposing m_local=m makes it go sequential.
> > > >
> > > > Cheers,
> > > >
> > > > Quentin CHEVALIER – IA parcours recherche
> > > > LadHyX - Ecole polytechnique
> > > > __________
> > > >
> > > >
> > > >
> > > > On Thu, 5 May 2022 at 10:34, Quentin Chevalier
> > > > <quentin.chevalier at polytechnique.edu> wrote:
> > > >>
> > > >> Hello all and thanks for your great work in bringing this very
> helpful package to the community !
> > > >>
> > > >> That said, I wouldn't need this mailing list if everything was
> running smoothly. I have a rather involved eigenvalue problem that I've
> been working on that's been throwing a mysterious error :
> > > >>>
> > > >>> petsc4py.PETSc.Error: error code 77
> > > >>>
> > > >>> [1] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:149
> > > >>> [1] EPSSolve_KrylovSchur_Default() at
> /usr/local/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c:289
> > > >>> [1] EPSGetStartVector() at
> /usr/local/slepc/src/eps/interface/epssolve.c:824
> > > >>> [1] Petsc has generated inconsistent data
> > > >>> [1] Initial vector is zero or belongs to the deflation space
> > > >>
> > > >>
> > > >> This problem occurs in parallel with two processors, using the
> petsc4py library using the dolfinx/dolfinx docker container. I have PETSc
> version 3.16.0, in complex mode, python 3, and I'm running all of that on a
> OpenSUSE Leap 15.2 machine (but I think the docker container has a Ubuntu
> OS).
> > > >>
> > > >> I wrote a minimal working example below, but I'm afraid the process
> for building the matrices is involved, so I decided to directly share the
> matrices instead :
> https://seminaris.polytechnique.fr/share/s/ryJ6L2nR4ketDwP
> > > >>
> > > >> They are in binary format, but inside the container I hope someone
> could reproduce my issue. A word on the structure and intent behind these
> matrices :
> > > >>
> > > >> QE is a diagonal rectangular real matrix. Think of it as some sort
> of preconditioner
> > > >> L is the least dense of them all, the only one that is complex, and
> in order to avoid inverting it I'm using two KSPs to compute solve problems
> on the fly
> > > >> Mf is a diagonal square real matrix, its on the right-hand side of
> the Generalised Hermitian Eigenvalue problem (I'm solving
> QE^H*L^-1H*L^-1*QE*x=lambda*Mf*x
> > > >>
> > > >> Full MWE is below :
> > > >>
> > > >> from petsc4py import PETSc as pet
> > > >> from slepc4py import SLEPc as slp
> > > >> from mpi4py.MPI import COMM_WORLD
> > > >>
> > > >> # Global sizes
> > > >> m_local=COMM_WORLD.rank*(490363-489780)+489780
> > > >> n_local=COMM_WORLD.rank*(452259-451743)+451743
> > > >> m=980143
> > > >> n=904002
> > > >>
> > > >> QE=pet.Mat().createAIJ([[m_local,m],[n_local,n]])
> > > >> L=pet.Mat().createAIJ([[m_local,m],[m_local,m]])
> > > >> Mf=pet.Mat().createAIJ([[n_local,n],[n_local,n]])
> > > >>
> > > >> viewerQE = pet.Viewer().createMPIIO("QE.dat", 'r', COMM_WORLD)
> > > >> QE.load(viewerQE)
> > > >> viewerL = pet.Viewer().createMPIIO("L.dat", 'r', COMM_WORLD)
> > > >> L.load(viewerL)
> > > >> viewerMf = pet.Viewer().createMPIIO("Mf.dat", 'r', COMM_WORLD)
> > > >> Mf.load(viewerMf)
> > > >>
> > > >> QE.assemble()
> > > >> L.assemble()
> > > >>
> > > >> KSPs = []
> > > >> # Useful solvers (here to put options for computing a smart R)
> > > >> for Mat in [L,L.hermitianTranspose()]:
> > > >> KSP = pet.KSP().create()
> > > >> KSP.setOperators(Mat)
> > > >> KSP.setFromOptions()
> > > >> KSPs.append(KSP)
> > > >> class LHS_class:
> > > >> def mult(self,A,x,y):
> > > >> w=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)
> > > >> z=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)
> > > >> QE.mult(x,w)
> > > >> KSPs[0].solve(w,z)
> > > >> KSPs[1].solve(z,w)
> > > >> QE.multTranspose(w,y)
> > > >>
> > > >> # Matrix free operator
> > > >> LHS=pet.Mat()
> > > >> LHS.create(comm=COMM_WORLD)
> > > >> LHS.setSizes([[n_local,n],[n_local,n]])
> > > >> LHS.setType(pet.Mat.Type.PYTHON)
> > > >> LHS.setPythonContext(LHS_class())
> > > >> LHS.setUp()
> > > >>
> > > >> # Eigensolver
> > > >> EPS = slp.EPS(); EPS.create()
> > > >> EPS.setOperators(LHS,Mf) # Solve QE^T*L^-1H*L^-1*QE*x=lambda*Mf*x
> > > >> EPS.setProblemType(slp.EPS.ProblemType.GHEP) # Specify that A is
> hermitian (by construction), and B is semi-positive
> > > >> EPS.setWhichEigenpairs(EPS.Which.LARGEST_MAGNITUDE) # Find largest
> eigenvalues
> > > >> EPS.setFromOptions()
> > > >> EPS.solve()
> > > >>
> > > >> Quentin CHEVALIER – IA parcours recherche
> > > >>
> > > >> LadHyX - Ecole polytechnique
> > > >>
> > > >> __________
> > >
> > > <image003.jpg><MWE.py>
> >
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220506/05e83b72/attachment-0001.html>
More information about the petsc-users
mailing list