Speeding up eigenvalue computation

Linear Algebra
1

Hello!
I am trying to solve a 2D eigenvalue problem in Fenics using SLEPc (the Schrodinger equation, to be more precise). I’d like to use a uniform mesh of size say 500 x 500, with polynomials of degree 5, for example. I’m looking for the smallest, say 100 eigenpairs. I find that this computation takes extremely long (at least 14 hours), or worse, does not happen at all! (It works for smaller mesh sizes). Is there a way to speed up this computation by using MPI perhaps or some other technique?

Below is a snippet of my code that outlines the eigensolver settings.

#Set up variational formulation
u = TrialFunction(Vh)
v = TestFunction(Vh)
a = inner(grad(u), grad(v))*dx + V*u*v*dx
b = u*v*dx

#Assemble stiffness and mass matrices.
A = PETScMatrix()
assemble(a, tensor=A)
B = PETScMatrix()
assemble(b, tensor=B)

#Set up SLEPc eigenvalue solver
eps = SLEPc.EPS()
eps.create()
eps.setOperators(A.mat(), B.mat())
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setTarget(0.)
eps.setDimensions(nreq)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
ksp = st.getKSP()
ksp.setType(PETSc.KSP.Type.PREONLY)
pc = ksp.getPC()
pc.setType(PETSc.PC.Type.CHOLESKY)

#Solve the eigenvalue problem
start = time.time()
eps.solve()
nconv = eps.getConverged()
print("Number of eigenpair converged = {}".format(nconv))
end = time.time()
print("Time elapsed "+str(end-start)+" s")

Thanks in advance!

2

Look into scalable methods rather than direct factorisation. E.g. Jacobi-Davidson iteration with an appropriate scalable preconditioner.

3

I’ve been quite unsuccessful trying to improve the speed of eigenvalue computations. See e.g. the following minimal example. Is there any standard configuration that works similarly to the default solvers in commercial packages?

shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('precond')  # spectral transform
shift.setShift(0)  # spectral shift
ksp = shift.getKSP()
ksp.setType(PETSc.KSP.Type.GMRES) 
pc = ksp.getPC()
#pc.setType(PETSc.PC.Type.CHOLESKY)
pc.setType(PETSc.PC.Type.ILU)

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=5)

eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eigensolver.setST(shift)
#eigensolver.setType('jd')  # Solver type
eigensolver.setType('lobpcg')  # Solver type
eigensolver.setWhichEigenpairs(eigensolver.Which.SMALLEST_REAL)  # Target Real
eigensolver.setOperators(stiffnessMatrix,massMatrix)
eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = stiffnessMatrix.getVecs()
print( "Number of converged eigenpairs %d" % evs )
threshold = 1
if evs > 0:
    for i in range (evs):
            l = eigensolver.getEigenpair(i ,vr, vi)