Hi, here I did something similar to what you are asking. Blocks are extracted from cbc.block_mat
and PETSc.MatNest
object is created that is then used by the eigenvalue solver from slepc4py
. Alternatively, for many eigenvalue solvers only the action of the operator is required. This is illustrated here (A is a matrix, W is a list of PETSc.Vec)
class EnergyNorm(object):
def __init__(self, A, W):
self.A = A
self.W = W
def mult(self, mat, x, y):
alphas = [wi.dot(x) for wi in self.W]
self.A.mult(x, y)
for alpha, wi in zip(alphas, self.W):
y.axpy(alpha, wi)
energy_norm = EnergyNorm(A, W)
# need dummy vector for size
vec = assemble(inner(Constant((1, 1, 1)), v)*dx)
vec = as_backend_type(vec).vec()
from petsc4py import PETSc
E = PETSc.Mat().createPython([vec.getSizes(), vec.getSize()],
comm=mpi_comm_world())
E.setPythonContext(energy_norm)
E.setUp()
from slepc4py import SLEPc
# Setup the eigensolver
solver = SLEPc.EPS().create()
solver.setOperators(E, B.mat())
solver.setType(solver.Type.KRYLOVSCHUR)
solver.setDimensions(3, PETSc.DECIDE)