I'm using H2lib for preconditioning my FEM matrices, and I noticed that for a time the memory usage is very high, even though (after compressing) the matrix is quite small. The problem appears because in order to compress, I first have to allocate the uncompressed matrix which has too many dense entries.
In order to work around this issue, I have now added two small functions to my local copy of H2lib (mostly copied from the build_from_block and copy_sparsematrix functions already there), which may be of interest for others. The idea is to first initialize the matrix without allocating memory for the dense subblocks, and then do the copying of the sparse entires and the compression at the same step. That way only one of the uncompressed subblocks needs to fit into memory after the other. This reduces peak memory usage by large factor in my usecase.
phmatrix
build_from_block_hmatrix_noinit(pcblock b, uint k)
{
phmatrix h, h1;
pcblock b1;
int rsons, csons;
int i, j;
h = NULL;
if (b->son) {
rsons = b->rsons;
csons = b->csons;
h = new_super_hmatrix(b->rc, b->cc, rsons, csons);
for (j = 0; j < csons; j++) {
for (i = 0; i < rsons; i++) {
b1 = b->son[i + j * rsons];
h1 = build_from_block_hmatrix_noinit(b1, k);
ref_hmatrix(h->son + i + j * rsons, h1);
}
}
}
else if (b->a > 0)
h = new_rk_hmatrix(b->rc, b->cc, k);
else{
//printf("initializing empty\n");
h = new_hmatrix(b->rc, b->cc);
h->f = new_amatrix(0,0);
h->desc = 1;
}
update_hmatrix(h);
return h;
}
void
copy_sparsematrix_coarse_hmatrix(psparsematrix sp, phmatrix hm, ptruncmode tm,double eps)
{
uint i, j, k, l, m, rsize, csize;
uint *ridx, *cidx, *row, *col;
pfield coeff;
pamatrix f;
if (hm->son) {
for (j = 0; j < hm->csons; j++) {
for (i = 0; i < hm->rsons; i++) {
copy_sparsematrix_coarse_hmatrix(sp, hm->son[i + j * (hm->rsons)],tm,eps);
}
}
}
else if (hm->r) { /*hm->r is defined, but hm->r->A = hm->r->B = NulL => all entries are 0 */
}
else if (hm->f) {
rsize = hm->rc->size;
csize = hm->cc->size;
ridx = hm->rc->idx;
cidx = hm->cc->idx;
row = sp->row;
col = sp->col;
coeff = sp->coeff;
if(hm->f->rows==0 && hm->f->cols==0) {
del_amatrix(hm->f);
hm->f = new_amatrix(rsize, csize);
}
f = hm->f;
for (j = 0; j < f->cols; j++)
for (i = 0; i < f->rows; i++)
f->a[i + j * f->rows] = 0.0;
for (k = 0; k < rsize; k++) {
i = ridx[k];
for (m = row[i]; m < row[i + 1]; m++) {
for (l = 0; l < csize; l++) {
if (col[m] == cidx[l]) {
setentry_amatrix(hm->f, k, l, coeff[m]);
}
}
}
}
}
//printf("eps:%.6e",eps);
coarsen_hmatrix(hm, tm, eps, false);
}