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 <[email protected]> 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 > > <[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/ >
