Creating a PETSc Matrix from given Matrix, solve eigenvalue problem (petsc4py?)

Linear Algebra
1

Hi everyone,

I’m using macOS Mojave (10.14.6), miniconda, and FEniCS 2018 with Python.
I have a Matrix (to be more precise the stiffness and the mass matrix) as a Numpy array and now I want to plug it into a PETSc matrix to use the eigenvalue solver. What I did so far is to use petsc4py and (I think that I) was able to create a PETSc matrix from my given Numpy arrays. Now there is a problem with the SLEPcEigenSolver:

from fenics import *
import numpy as np
import petsc4py.PETSc as pet

# preliminaries
N = 25
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, 'Lagrange', 1)
Nh = V.dim()
bc = DirichletBC(V, Constant(0.0), 'on_boundary')

bdry_dofs = bc.get_boundary_values().keys()
first, last = V.dofmap().ownership_range()
all_dofs = range(last-first)
interior_dofs = list(set(all_dofs) - set(bdry_dofs))
boundary_dofs = list(bdry_dofs)

v = TrialFunction(V)
z = TestFunction(V)
a = dot(grad(v), grad(z))*dx
m = v*z*dx

# construct matrix
A, M = PETScMatrix(), PETScMatrix()
A = assemble(a, tensor = A);
M = assemble(m, tensor = M);

# delete zero boundary conditions
A_mat = np.delete(np.delete(A.array(), boundary_dofs, axis = 0), boundary_dofs, axis = 1)
M_mat = np.delete(np.delete(M.array(), boundary_dofs, axis = 0), boundary_dofs, axis = 1)

A = pet.Mat().create()
A.setSizes(A_mat.shape)
A.setType('aij')
A.setUp()
A.setValues(range(0, A_mat.shape[0]), range(0, A_mat.shape[0]), A_mat)
A.assemble()
#A.setType('dense')

M = pet.Mat().create()
M.setSizes(M_mat.shape)
M.setType('aij')
M.setUp()
M.setValues(range(0, M_mat.shape[0]), range(0, M_mat.shape[0]), M_mat)
M.assemble()
#M.setType('dense')

# eigensolver specs; solve eigenvalue problem
eigensolver = SLEPcEigenSolver(A, M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.solve(K)
conv = eigensolver.get_number_converged()
print('number of converged eigenvalues: {:3d}'.format(conv))

The error I get is:

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfin.cpp.la.SLEPcEigenSolver(arg0: dolfin.cpp.la.PETScMatrix)
    2. dolfin.cpp.la.SLEPcEigenSolver(arg0: dolfin.cpp.la.PETScMatrix, arg1: dolfin.cpp.la.PETScMatrix)
    3. dolfin.cpp.la.SLEPcEigenSolver(arg0: MPICommWrapper)

Invoked with: <petsc4py.PETSc.Mat object at 0x7f9b3037ba98>, <petsc4py.PETSc.Mat object at 0x7f9b11ce7eb8>

So I somehow have to convert the petsc4py into a dolphin PETScMatrix? Does anybody know how to do that. Or, even better: Can I bypass the construction of a petsc4py Matrix and use the built-in dolphin PETSc Matrix and plug my Matrix in?

Thanks for any hints/help.

2

Inserting

A = PETScMatrix(A)
M = PETScMatrix(M)

after you created your matrices seems to work.

3

Thanks, that solved the problem. Do you know if I can do similar things without the petsc4py module and use the built in PETSc matrix? The setValues takes very long.

4

No, sorry I can’t help you with that.

5

A_petsc = as_backend_type(A).mat()