[petsc-users] Mysterious error code 77
Jose E. Roman
jroman at dsic.upv.es
Fri May 6 08:34:33 CDT 2022
MUMPS is for factorizing matrix Mf, i.e., EPS.Which.LARGEST_MAGNITUDE.
But if you do shift-and-invert you have to do it in a different way. LHS cannot be factorized because it is a shell matrix.
> El 6 may 2022, a las 15:28, Quentin Chevalier <quentin.chevalier at polytechnique.edu> escribió:
>
> 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
>
> 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>
> >
More information about the petsc-users
mailing list