I'm trying to define a matrix-free solver. Here is a very simple example with the identity matrix:
from dolfin import *
class IdentityMatrix(LinearOperator):
def __init__(self, V):
u = Function(V)
LinearOperator.__init__(self, u.vector(), u.vector())
def mult(self, x, y):
y.zero()
y.axpy(1.0, x)
if __name__ == "__main__":
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, 'Lagrange', 1)
Id = IdentityMatrix(V)
solver = PETScKrylovSolver("cg", "amg")
solver.set_operator(Id)
b = interpolate(Constant(1.0), V)
x = Function(V)
solver.solve(x.vector(), b.vector())
This crashes with the following error message:
$ python test_matrixfree.py
Solving linear system of size 121 x 121 (PETSc Krylov solver).
Traceback (most recent call last):
File "test_matrixfree.py", line 23, in <module>
solver.solve(x.vector(), b.vector())
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 56.
*** Where: This error was encountered inside /build/dolfin-k_QrtL/dolfin-1.6.0/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 1.6.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
Same error occurs with other preconditioners. However, this works no preconditioner is used,
solver = PETScKrylovSolver("cg", "none")
Isn't there any way to define a preconditioner in that case?
I use the Ubuntu binaries for fenics 1.6.0, but I observed the same error with a compiled version of fenics 1.6.0.
Thanks for any help!
