[petsc-users] Mysterious error code 77

Fri May 6 09:43:31 CDT 2022

Ok thank you for the clarification. It would appear my MWE was poorly
chosen indeed, I got it to work too. I still get code 77 inside my
main code, I'll keep you posted.

Indeed the matrices are stored for testing purposes only.

Of course you are completely right about shift invert, I was confusing
with another case where I had set a target eigenvalue.

Thanks for your time,

Quentin CHEVALIER – IA parcours recherche
LadHyX - Ecole polytechnique
__________

On Fri, 6 May 2022 at 15:55, Jose E. Roman <jroman at dsic.upv.es> wrote:
>
> I used your script. The only change I did was
>
>     dir+="seq/"
>
> instead of
>
>     dir+="par/"
>
> That is, read the sequential matrices. I would suggest writing matrix files from a sequential program. I assume that matrix files are for testing only, for production runs you should generate the matrices and solve in the same program, not save to file.
>
> The original script is doing EPS.setWhichEigenpairs(EPS.Which.LARGEST_MAGNITUDE). This computes largest eigenvalues.
> You cannot do ST.setType('sinvert') in this example. This would compute eigenvalues close to zero. Read chapter 3 of the manual to understand how shift-and-invert works.
>
> Jose
>
>
> > El 6 may 2022, a las 15:46, Quentin Chevalier <quentin.chevalier at polytechnique.edu> escribió:
> >
> > @Matthew, I'm sure you've noticed the L.solve inside the LHS class.
> > That's preventing me from going AIJ ; I would have to invert L, which
> > sounds like a bad idea.
> >
> > @Jose, I've had troubles reading the matrix files in parallel when
> > they were created in sequential or parallel with a different number of
> > processors - that's why I resorted to 2, and hard-coded the local
> > sizes. If I naively try to create AIJ with global sizes only, I get
> > Incompatible sizes every time. Would you care to enlighten me as to
> > how to read the same matrix files with an arbitrary number of
> > processors ?
> >
> > More importantly, you managed to run the MWE sent on May 5th without
> > any change uin the same config as me and couldn't reproduce a code 77
> > ?
> >
> > You're also telling me that Mf needs to be factorised with MUMPS, but
> > that LHS can't. So did you specify an additional solver/preconditioner
> > inside the MWE or did it work as is ?
> >
> > Cheers,
> >
> > Quentin CHEVALIER – IA parcours recherche
> > LadHyX - Ecole polytechnique
> > __________
> >
> >
> > On Fri, 6 May 2022 at 15:35, Matthew Knepley <knepley at gmail.com> wrote:
> >>
> >> 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/
>