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
> <[email protected]> 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 <[email protected]> wrote:
>>
>> On Fri, May 6, 2022 at 9:28 AM Quentin Chevalier
>> <[email protected]> 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 <[email protected]> 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
>>>>> <[email protected]> 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 <[email protected]> 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
>>>>>> <[email protected]> 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
>>>>>> <[email protected]> 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/
