I could reproduce it and I think I found it, @joachim I sent you a merge request for the fix
The problem was corrupt memory, valgrind gave some like these:
==14882== Invalid write of size 8
==14882== at 0x1B45829B: operator= (vector.hpp:740)
==14882== by 0x1B45829B: ngmg::LinearProlongation::ProlongateInline(int, ngla::BaseVector&) const (prolongation.cpp:55)
==14882== by 0x1B45456C: ngmg::MultigridPreconditioner::MGM(int, ngla::BaseVector&, ngla::BaseVector const&, int) const (mgpre.cpp:265)
==14882== by 0x1B453D70: ngmg::MultigridPreconditioner::Mult(ngla::BaseVector const&, ngla::BaseVector&) const (mgpre.cpp:173)
==14882== by 0x1B457042: void ngla::VMatVecExpr::AssignTo<double>(double, ngla::BaseVector&) const (basematrix.hpp:209)
==14882== by 0x1B456C38: void ngla::VVecExpr<ngla::VMatVecExpr>::AssignTo<double>(double, ngla::BaseVector&) const (basevector.hpp:64)
==14882== by 0x1B45686D: ngla::BaseVector& ngla::AutoVector::operator=<ngla::VMatVecExpr>(ngla::VVecExpr<ngla::VMatVecExpr> const&) (basevector.hpp:370)
==14882== by 0x1B454C1E: ngmg::TwoLevelMatrix::Mult(ngla::BaseVector const&, ngla::BaseVector&) const (mgpre.cpp:328)
==14882== by 0x1B457042: void ngla::VMatVecExpr::AssignTo<double>(double, ngla::BaseVector&) const (basematrix.hpp:209)
==14882== by 0x1B456C38: void ngla::VVecExpr<ngla::VMatVecExpr>::AssignTo<double>(double, ngla::BaseVector&) const (basevector.hpp:64)
==14882== by 0x1B45686D: ngla::BaseVector& ngla::AutoVector::operator=<ngla::VMatVecExpr>(ngla::VVecExpr<ngla::VMatVecExpr> const&) (basevector.hpp:370)
==14882== by 0x261170C1: ngla::CGSolver<std::complex<double> >::Mult(ngla::BaseVector const&, ngla::BaseVector&) const (cg.cpp:588)
==14882== by 0x1CAB010A: ngsolve::NumProcBVP::Do(ngstd::LocalHeap&)::{lambda()#1}::operator()() const (bvp.cpp:389)
==14882== Address 0x1a1131b0 is 0 bytes after a block of size 2,064 alloc'd
==14882== at 0x4C3089F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==14882== by 0x25FE2484: ngla::S_BaseVectorPtr<std::complex<double> >::S_BaseVectorPtr(unsigned long, int) (vvector.hpp:33)
==14882== by 0x25FE317E: ngla::VVector<std::complex<double> >::VVector(unsigned long) (vvector.hpp:150)
==14882== by 0x26351461: void __gnu_cxx::new_allocator<ngla::VVector<std::complex<double> > >::construct<ngla::VVector<std::complex<double> >, int const&>(ngla::VVector<std::complex<double> >*, int const&) (new_allocator.h:136)
==14882== by 0x263511BA: void std::allocator_traits<std::allocator<ngla::VVector<std::complex<double> > > >::construct<ngla::VVector<std::complex<double> >, int const&>(std::allocator<ngla::VVector<std::complex<double> > >&, ngla::VVector<std::complex<double> >*, int const&) (alloc_traits.h:475)
==14882== by 0x26350347: std::_Sp_counted_ptr_inplace<ngla::VVector<std::complex<double> >, std::allocator<ngla::VVector<std::complex<double> > >, (__gnu_cxx::_Lock_policy)2>::_Sp_counted_ptr_inplace<int const&>(std::allocator<ngla::VVector<std::complex<double> > >, int const&) (shared_ptr_base.h:526)
==14882== by 0x2634B326: std::__shared_count<(__gnu_cxx::_Lock_policy)2>::__shared_count<ngla::VVector<std::complex<double> >, std::allocator<ngla::VVector<std::complex<double> > >, int const&>(std::_Sp_make_shared_tag, ngla::VVector<std::complex<double> >*, std::allocator<ngla::VVector<std::complex<double> > > const&, int const&) (shared_ptr_base.h:637)
==14882== by 0x263414B1: std::__shared_ptr<ngla::VVector<std::complex<double> >, (__gnu_cxx::_Lock_policy)2>::__shared_ptr<std::allocator<ngla::VVector<std::complex<double> > >, int const&>(std::_Sp_make_shared_tag, std::allocator<ngla::VVector<std::complex<double> > > const&, int const&) (shared_ptr_base.h:1295)
==14882== by 0x263381F0: std::shared_ptr<ngla::VVector<std::complex<double> > >::shared_ptr<std::allocator<ngla::VVector<std::complex<double> > >, int const&>(std::_Sp_make_shared_tag, std::allocator<ngla::VVector<std::complex<double> > > const&, int const&) (shared_ptr.h:344)
==14882== by 0x2633092B: std::shared_ptr<ngla::VVector<std::complex<double> > > std::allocate_shared<ngla::VVector<std::complex<double> >, std::allocator<ngla::VVector<std::complex<double> > >, int const&>(std::allocator<ngla::VVector<std::complex<double> > > const&, int const&) (shared_ptr.h:691)
==14882== by 0x26316303: std::shared_ptr<ngla::VVector<std::complex<double> > > std::make_shared<ngla::VVector<std::complex<double> >, int const&>(int const&) (shared_ptr.h:707)
==14882== by 0x26389962: ngla::SparseMatrix<std::complex<double>, std::complex<double>, std::complex<double> >::CreateVector() const (sparsematrix.cpp:1357)
==14882==
The problem was in the wrong size of the Range of a FlatSysVector, vector.hpp:766 ff:
INLINE FlatSysVector<T> Range(size_t first, size_t last)
{ return FlatSysVector<T> (last-first+1, blocksize, data+(first*blocksize)); }
INLINE /* const */ FlatSysVector<T> Range(size_t first, size_t last) const
{ return FlatSysVector<T> (last-first+1, blocksize, data+(first*blocksize)); }
Best
Christopher