@vv From Pernet & Stain "...There are implementations of Hermite normal form algorithms in NTL [Sho], PARI [PAR], GAP [GAP], Maple and Mathematica. The algorithm in this paper is asymptotically better than these standard implementations..."

@vv Still I would like to have some references. I am not going to implement anything, by the way. In the case of Sage Math it is well documented and you can find that the default algorithm for Hermite normal form in Sage is the one proposed by Clément Pernet and William Stein (a variant of the algorithm proposed by Daniele Micciancio and Bogdan Warinschi, in fact).

@vv Sure, but the documentation, as I see, does not say which algorithms are implemented.

@roman_pearce Long time I did not write anything. Anyway, the nice solution is to use the fgb_gbasis_elim procedure from FGb library.

@tomleslie This solution works! Thanks!

@tomleslie Sorry but this does not work as it should since I do not want to flatten an array.

OK I do not know how but I pusted wrong system. Rigth now is updated. Sorry guys!

Sorry guys I posted this from a hotel room and I think, I made a mistake somewhere. I will check it later and repost if necessery.

@acer Thanks! It was so simple.

OK I think I found the source of the problem. First of all I have in my code these lines to save data to a file

fileID := fopen(sprintf("sam_%d.csv", id), APPEND, TEXT):
writedata(fileID, Threads:-Map(proc (x) return [midpoint, op(ExtractSamplePoints(x))] end proc, samplePoints), string, proc (f, x) fprintf(f, %a, x) end proc):
fclose(fileID):

After I commented them out memory usage was fine. This made me think if somehow printer procedure used in Grid:-Launch is not related the problem. Indeed, after I defined it as a procedure which returns NULL my memory problem seems to be solved. Maybe someone could provide a nicer solution or more details why it even occurs?

Right now

%MEM     TIME+ COMMAND

1.9  29:26.28 mserver 

and not 49.1 like before.

@roman_pearce Thank you for your nice reply. I've worked on this issue a bit more and after I replaced UnivariatePolynomial by another variable elimination strategy which uses resultants (see below) it starts to work. I thought at some point that the problem is caused by remember tables but I added calls to forget and this did not work out. So, I guess C libraries are responsible for the problem. On the other hand I sent a mail to the technical support to notice possible issue with memory. Well, we will see how this will end up. I will post here any answer I will get.

EliminationResultant := proc(S, vars)
local r1, r2, r3, rr1, rr2, rr3;
r1 := resultant(S[1], S[2], vars[3]):
r2 := resultant(S[1], S[3], vars[3]):
r3 := resultant(S[2], S[3], vars[3]):
rr1 := resultant(r1, r2, vars[2]):
rr2 := resultant(r1, r3, vars[2]):
rr3 := resultant(r2, r2, vars[2]):
return foldl(gcd, rr1, rr2, rr3):
end proc:

PS With the input which I gave you the original code needed around 11h on the machine with 20 dual core processors. With EliminationResultant  it needs around 5h if I remember correctly.

@roman_pearce On the other hand, could you tell me why UnivariatePolynomial procedure can possible causes memory problems. Should not the garbage collector take care of that even we have millions of calls to it?

@roman_pearce I just realized that I forgot to past the output from CodeTools:-Usage for some data set similar to what I've posted here.

memory used=2.13GiB, alloc change=4.43GiB, cpu time=31.81m, real time=10.62h, gc time=31.09m

This is interesting because does not present what I've observed. I suppose that CodeTools does not play nice with Grid framework. By the way I use maple 2015.1 64bit for Linux.

@roman_pearce Thanks for your feedback. Here some data dataQ.txt

@Thomas Richard Thanks for recommendation of the package PolyhedralSets. It looks very promising.