Possible bug in SLEPcEigenSolver

[Christian Waluga]

Today we stumbled upon a possible bug in the SLEPcEigenSolver. A minimal example to reproduce this is attached.

When we call eigensolver.solve() in the code listed below, PETSc seems to try to allocate 16 Exabytes, which I unfortunately don't have :) .

Python(13477) malloc: *** mmap(size=18446744059458596864) failed (error code=12)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Out of memory. This could be due to allocating
[0]PETSC ERROR: too large an object or bleeding by not properly
[0]PETSC ERROR: destroying unneeded objects.
[0]PETSC ERROR: Memory allocated 0 Memory used by process 0
[0]PETSC ERROR: Try running with -malloc_dump or -malloc_log for info.
[0]PETSC ERROR: Memory requested 18446744059458592768!
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 3, Wed Aug 29 11:26:24 CDT 2012 
[0]PETSC ERROR: See docs/changes/index.html for recent updates.
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[0]PETSC ERROR: See docs/index.html for manual pages.
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Unknown Name on a darwin11. named feynman.local by waluga Wed Jul 17 13:13:39 2013
[0]PETSC ERROR: Libraries linked from /Users/johannr/fenics-snapshot/local/lib
[0]PETSC ERROR: Configure run at Tue May 14 10:31:34 2013
[0]PETSC ERROR: Configure options --prefix=/Users/johannr/fenics-snapshot/local COPTFLAGS=-O2 --with-debugging=0 --with-clanguage=cxx --with-c-support=1 --with-blas-lapack-dir=/usr --with-umfpack=1 --with-umfpack-include=/Users/johannr/fenics-snapshot/local/include/suitesparse --with-umfpack-lib="[/Users/johannr/fenics-snapshot/local/lib/libumfpack.a,/Users/johannr/fenics-snapshot/local/lib/libamd.a]" --with-spooles=1 --with-spooles-include=/Users/johannr/fenics-snapshot/local/include --with-spooles-lib=/Users/johannr/fenics-snapshot/local/lib/libspooles.a --with-ptscotch=1 --with-ptscotch-dir=/Users/johannr/fenics-snapshot/local --with-ml=1 --with-ml-include=/Users/johannr/fenics-snapshot/local/include/trilinos --with-ml-lib=/Users/johannr/fenics-snapshot/local/lib/libml.dylib --with-hdf5=1 --with-hdf5-dir=/Users/johannr/fenics-snapshot/local --with-x=0 -with-x11=0 --with-fortran=0 --with-shared-libraries=1 PETSC_DIR=/Users/johannr/fenics-snapshot/fenics-superbuild/build-fenics/CMakeExternals/src/PETSc PETSC_ARCH=darwin11.4.0-cxx-opt
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: PetscMallocAlign() line 49 in /Users/johannr/fenics-snapshot/fenics-superbuild/build-fenics/CMakeExternals/src/PETSc/src/sys/memory/mal.c
[0]PETSC ERROR: VecDuplicateVecs_Contiguous() line 61 in src/vec/contiguous.c
[0]PETSC ERROR: VecDuplicateVecs() line 623 in /Users/johannr/fenics-snapshot/fenics-superbuild/build-fenics/CMakeExternals/src/PETSc/src/vec/vec/interface/vector.c
[0]PETSC ERROR: EPSAllocateSolution() line 51 in src/eps/interface/mem.c
[0]PETSC ERROR: EPSSetUp_KrylovSchur() line 137 in src/eps/impls/krylov/krylovschur/krylovschur.c
[0]PETSC ERROR: EPSSetUp() line 137 in src/eps/interface/setup.c
[0]PETSC ERROR: EPSSolve() line 109 in src/eps/interface/solve.c

The strange thing is that the unit square mesh of comparable size does not cause such problems. Also, when I change eigensolver.solve() to eigensolver.solve(0), everything seems to work fine. I suspect that PETSc takes care about the allocation of the missing memory automatically then.

My guess would be that in the call to solve(), the size is not properly obtained from the PETScMatrix before passing it to solve(size_t n). Then again, why does this occur only with some meshes? The coarser refinement levels of this particular mesh did not cause any problems...

BTW: The error could be produced with FEniCS 1.2.0+ in Mac OS X 10.8, as well as with 1.1 and 1.2 on two different Ubuntu systems.

from dolfin import *

lmbda = 10
mu = 8

#mesh = UnitSquareMesh(400,400)
mesh = Mesh('mesh.xml.gz')
print mesh.size(2)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0))

def sigma(v):
  return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)

a = inner(sigma(u), grad(v))*dx
L = inner(u,v)*dx

A = PETScMatrix()
M = PETScMatrix()
assemble(a, tensor=A)
assemble(L, tensor=M)

eigensolver = SLEPcEigenSolver(A,M)

print 'solving'
eigensolver.solve()
#eigensolver.solve(0)

[Prof Garth Wells]

Your problem is quite big, so you need to be careful with the solver parameters. Using

eigensolver = SLEPcEigenSolver(A, M)
eigensolver.parameters["problem_type"] = "gen_hermitian"
eigensolver.parameters["tolerance"] = 1.0e-6

print 'solving'
num = 4
eigensolver.solve(num)
for i in range(num):
    print("Eigenvalue:", i, eigensolver.get_eigenvalue(i))

I get four eigenvalues in less than 20s.