Hi everyone,

warning: I am new to Eigen (and C++) and am probably doing something stupid! ;-)

I implemented the most trivial algorithm (eGCD based) to compute the Hermite normal form of a matrix. It uses column transformations of a ColMajor SparseMatrix with mpz_class (GMP) entries. Code is here: https://github.com/1k5/hnf/blob/master/SparseHNF.cc

Now in the last commit, I changed from doing column transformation by multiplying with a matrix from the right to doing the column transformations by hand. See diff https://github.com/1k5/hnf/commit/2b979 ... 96c1b62189 . Obviously this is much faster. However, now this function eats memory like crazy. In the test directory lies a matrix of size about 750x1500 for which previously the memory usage during HNF computation was around 50MB, now it grows to 30GB and beyond.

Do any of you Eigen experts see from looking at the diff where my memory goes?

hm, which Eigen version are you using? Please give a try to the devel branch or 3.2 branch.

Also, your algorithm is modifying in-place the sparsity pattern of your sparse matrices which is really inefficient. First, you can avoid the swaps using a PermutationMatrix initialized to the identity and swaping the indices instead of the columns. Then; every time you want to access to column j, you first fetch its actual position through the permutation:

std::swap(perm.indices()(j),perm.indices()(k));
jp = perm.indices()(j)
H.col(jp)

At the complete end you can apply the permutation once: H = H*perm;

Then operations like H.col(j) = v+w are also rather costly in terms of memory reallocation/copies. You can limit it by reserving enough room per column, e.g.:

VectorXi c(H.cols());
for(...) c[j] = H.col(j).nonZeros();
H.reserve(3*c); // replace 3 with a good heuristic about the fill-in ratio.

ggael wrote:hm, which Eigen version are you using? Please give a try to the devel branch or 3.2 branch.

I'm using devel branch as of 2 days ago as mirrored on github here: https://github.com/RLovelett/eigen . I'll try again with direct download from the mercurial repo.

Also, your algorithm is modifying in-place the sparsity pattern of your sparse matrices which is really inefficient. First, you can avoid the swaps using a PermutationMatrix initialized to the identity and swaping the indices instead of the columns. Then; every time you want to access to column j, you first fetch its actual position through the permutation:

std::swap(perm.indices()(j),perm.indices()(k));
jp = perm.indices()(j)
H.col(jp)

At the complete end you can apply the permutation once: H = H*perm;

Hmm, okay I see. Then I'd need to apply the (col) permutation for H to the rows of the matrix in step (3). That should be easy to do.
So far the algorithm is the most naive implementation of the algorithm described on wikipedia on Smith normal form, slightly changed for Hermite NF. At least for dense matrices much better algorithms exist, so that'd be more interesting to look at first.

Then operations like H.col(j) = v+w are also rather costly in terms of memory reallocation/copies. You can limit it by reserving enough room per column, e.g.:

VectorXi c(H.cols());
for(...) c[j] = H.col(j).nonZeros();
H.reserve(3*c); // replace 3 with a good heuristic about the fill-in ratio.

Yes. The reason I do H.prune(0,0) now is to get rid of explicit zeros so that it is easy to find a pivot point. Probably replacing the pivot code with a loop with innerIterator() and ignoring explicit zeros would be faster since it allows keeping some extra space around.

But, none of the above explains why the changes I made (see diff from my first post) cause such a huge memory leak. I'll try with the latest Eigen version from mercurial and post an update.

And thanks a lot for helping (again)!

yes, my suggestions were only about making your code more efficient, not related to the memory leaks you are observing.

First think to check is how you defined Eigen::NumTraits< mpz_class >, make sure that RequireInitialization=1. Then, running your program within a memory debugger like valgring could help to spot the issue.