@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/
