[petsc-users] Matrix-free generalised eigenvalue problem

Tue Jul 18 02:57:45 CDT 2023

Matrix to matrix products are taking much longer than expected... My
snippet is below. m and n are quite large, >1M each. I'm running this on 35
procs. As you can see U, S and V are quite sparse SVD matrices (only their
first 4 columns are dense, plus a chop). I expected therefore approximate R
to have only rank 4 and computations to run smoothly once the cross
products between U and V are computed... Right now my bottle neck is not
preconditioner but getting that approximation of M2. What do you think ?

def approximate(self,k:int):
m, n = self.R.getSize()
m_local,n_local = self.R.getLocalSize()
I,J=self.tmp3.vector.getOwnershipRange()
S=pet.Mat().create(comm=comm); S.setSizes([[n_local,n],[n_local,n]])
U=pet.Mat().create(comm=comm); U.setSizes([[m_local,m],[n_local,n]])
V=pet.Mat().create(comm=comm); V.setSizes([[n_local,n],[n_local,n]])
for A in (U,S,V): A.setUp()
g=np.loadtxt("gains/txt/dummy.txt")
for i in range(k):
S.setValue(i,i,g[i])
loadStuff("left-vecs/", params,self.tmp3)
U.setValues(self.reindex,[i],self.tmp3.x.array)
loadStuff("right-vecs/",params,self.tmp3)
V.setValues(range(I,J), [i],self.tmp3.vector)
for A in (U,S,V): A.assemble()
U.chop(1e-6); V.chop(1e-6)

V.hermitianTranspose(V)
self.P.matMult(U.matMult(S.matMult(V,S),U),U) # Multiplications in place
(everyone is square except U)
M2=self.N.duplicate()
M2.matMult(U,U)
U.hermitianTranspose().matMult(U,M2)
return self.M+M2

Quentin

On Tue, 18 Jul 2023 at 07:48, Quentin Chevalier <
quentin.chevalier at polytechnique.edu> wrote:
>
> Thank you for that pointer ! I have on hand a partial SVD of R, so I used
that to build the approximate matrix instead.
>
> It's great that so many nice features of PETSc like
STSetPreconditionerMat are accessible through petsc4py !
>
> Good day,
>
> Quentin
>
>
>
> On Mon, 17 Jul 2023 at 17:29, Jose E. Roman <jroman at dsic.upv.es> wrote:
> >
> > It is possible to pass a different matrix to build the preconditioner.
That is, the shell matrix for B (EPSSetOperators) and an explicit matrix
(that approximates B) for the preconditioner. For instance, you can try
passing M for building the preconditioner. Since M is an explicit matrix,
you can try the default preconditioner (block Jacobi with ILU as local
solver) or even a full LU decomposition. The effectiveness of the
preconditioner will depend on how the update M+R^H P M P R moves the
eigenvalues around.
> >
> > You can do this with STSetSplitPreconditioner() or
STSetPreconditionerMat(). In your case any of them will do.
> >
> > Jose
> >
> >
> > > El 17 jul 2023, a las 15:50, Quentin Chevalier <
quentin.chevalier at polytechnique.edu> escribió:
> > >
> > > Thank you for this suggestion, I tried to implement that but it's
> > > proven pretty hard to implement MATOP_GET_DIAGONAL without completely
> > > tanking performance. After all, B is a shell matrix for a reason : it
> > > looks like M+R^H P M P R with R itself a shell matrix.
> > >
> > > Allow me to point out that I have no shift. My eigenvalue problem is
> > > purely about the largest ones out there. Section 8.2 and 3.4.3 led me
> > > to think that there was a way to avoid computing (or writing a shell
> > > matrix about it) B^-1... But you seem to stress that there's no way
> > > around it.
> > >
> > > Quentin
> > >
> > >
> > >
> > > On Mon, 17 Jul 2023 at 11:56, Jose E. Roman <jroman at dsic.upv.es>
wrote:
> > >>
> > >> The B-inner product is independent of the ST operator. See Table
3.2. In generalized eigenproblems you always have an inverse.
> > >>
> > >> If your matrix is diagonally dominant, try implementing the
MATOP_GET_DIAGONAL operation and using PCJACOBI. Apart from this, you have
to build your own preconditioner.
> > >>
> > >> Jose
> > >>
> > >>
> > >>> El 17 jul 2023, a las 11:48, Quentin Chevalier <
quentin.chevalier at polytechnique.edu> escribió:
> > >>>
> > >>> Hello Jose,
> > >>>
> > >>> I guess I expected B to not be inverted but instead used as a mass
for a problem-specific inner product since I specified GHEP as a problem
type. p50 of the same user manual seems to imply that that would indeed be
the case. I don't see what problem there would be with using a shell B
matrix as a weighting matrix, as long as a mat utility is provided of
course.
> > >>>
> > >>> I tried the first approach - I set up my KSP as CG since B is
hermitian positive-definite (I made a mistake in my first email), but I'm
getting a KSPSolve has not converged, reason DIVERGED_ITS error. I'm
letting it run for 1000 iterations already so it seems suspiciously slow
for a CG solver.
> > >>>
> > >>> I'm grappling with a shell preconditioner now to try and speed it
up, but I'm unsure which one allows for shell matrices.
> > >>>
> > >>> Thank you for your time,
> > >>>
> > >>> Quentin
> > >>>
> > >>>
> > >>> On Wed, 12 Jul 2023 at 19:24, Jose E. Roman <jroman at dsic.upv.es>
wrote:
> > >>>>
> > >>>> By default, it is solving the problem as B^{-1}*A*x=lambda*x (see
chapter on Spectral Transformation). That is why A can be a shell matrix
without problem. But B needs to be an explicit matrix in order to compute
an LU factorization. If B is also a shell matrix then you should set an
iterative solver for the associated KSP (see examples in the chapter).
> > >>>>
> > >>>> An alternative is to create a shell matrix M that computes the
action of B^{-1}*A, then pass M to the EPS solver as a standard
eigenproblem.
> > >>>>
> > >>>> Jose
> > >>>>
> > >>>>
> > >>>>> El 12 jul 2023, a las 19:04, Quentin Chevalier <
quentin.chevalier at polytechnique.edu> escribió:
> > >>>>>
> > >>>>> Hello PETSc Users,
> > >>>>>
> > >>>>> I have a generalised eigenvalue problem : Ax= lambda Bx
> > >>>>> I used to have only A as a matrix-free method, I used mumps and
an LU preconditioner, everything worked fine.
> > >>>>>
> > >>>>> Now B is matrix-free as well, and my solver is returning an error
: "MatSolverType mumps does not support matrix type python", which is
ironic given it seem to handle A quite fine.
> > >>>>>
> > >>>>> I have read in the user manual here that there some methods may
require additional methods to be supplied for B like MATOP_GET_DIAGONAL but
it's unclear to me exactly what I should be implementing and what is the
best solver for my case.
> > >>>>>
> > >>>>> A is hermitian, B is hermitian positive but not positive-definite
or real. Therefore I have specified a GHEP problem type to the EPS object.
> > >>>>>
> > >>>>> I use PETSc in complex mode through the petsc4py bridge.
> > >>>>>
> > >>>>> Any help on how to get EPS to work for a generalised matrix-free
case would be welcome. Performance is not a key issue here - I have a
tractable high value case on hand.
> > >>>>>
> > >>>>> Thank you for your time,
> > >>>>>
> > >>>>> Quentin
> > >>>>
> > >>
> >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230718/43c091d0/attachment.html>