Not an answer but too long for a comment.
Well, this additional memory consumption is actually not as bad as you might think, and I am not so sure whether this should be considered a memory leak. The additive assembler has
to sort the entries of Transpose[e[[1;;2]]]
;
to reorder the entries of e[[3]]
;
to compute the arrays for the row pointers;
to compute the arrays for the column indices; and
to accumulate the values from (the reordered) array e[[3]]
into the list of nonzeroes.
I guess one would fuse steps 3. to 5. together. Then one has to hold at the same time: arrays of the Dimensions
{r,2}
(for e[[1;;2]]
), {r}
(for e[[3]]
), {m+1}
(for the row pointers), {nnz,1}
(for the column indices) and {nnz}
(for the list of nonzero entries). Here nnz
is the number of nonzeroes in the final matrix. Since this is not known in the beginning, probable r
is used as upper bound so that memory of {r,1}
(for the column indices) and {r}
(for the list of nonzero entries) is allocated for the computations (and the tails are discarded when computations are done).
Also the sorting needs some scratch space - and here one can basically trade runtime for memory consumption (e.g. bucket sort has complexity $O(r)$ but requires quite a lot of memory while quicksort is in-place but has complexity $O(r \, \log(r))$. I don't know which particular sorting algorithm is used here. The memory consumption for the sorting should actually not add up to the one for the rest of the assembly. But maybe the scratch space for the sort was not deallocated directly after the sort but only later? (I find it unlikely, though.)
All this does probably not account for all the memory insufficiency, but for quite a lot of it. Further memory consumption may be due to the behavior of LibraryLink
: Typically it has to form a copy for returning the results. But that pretty much depends on the implementation details into which I have to insight.
Moreover, some parts of the assembly might be parallelized, e.g. with OpenMP. Then each thread probably needs some further scratch space to avoid write locks into shared memory. And that may also blow up the memory consumption dramatically.
Anyways, all this suggetes that one should try to figure out how to compute the column indices and row pointers directly and to generate the SparseArray
with the undocumented syntax
SparseArray @@ {Automatic, {m, n}, 0, {1, {rp, ci}, vals}}
SparseArray @@ {Automatic, {m, n}, 0., {1, {rp, ci}, vals}}
SparseArray @@ {Automatic, {m, n}, 0.+0.I, {1, {rp, ci}, vals}}
(depending on which scalar type you want), where rp
corresponds to b["RowPointers"]
, ci
corresponds to b["ColumnIndices"]
, and vals
corresponds to b["NonzeroValues"]
, provided that matrix b
is in proper CRS format (i.e., if SparseArray`SparseArraySortedQ[b]
returns True
).
Finally, I'd like to say that one should be careful when interpreting the result of MaxMemoryUsed
. Actually, the documentation is super vague about what it does and the sentence
On most computer systems, MaxMemoryUsed[] will give results close to
those obtained from external process status requests.
does not really raise any confidence.