diff -N -r -c HLib-1.3old/Examples/example_bem3d.c HLib-1.3/Examples/example_bem3d.c *** HLib-1.3old/Examples/example_bem3d.c 2004-12-12 17:42:30.000000000 +0100 --- HLib-1.3/Examples/example_bem3d.c 2005-05-09 11:16:23.281225144 +0200 *************** *** 295,300 **** --- 295,309 ---- M = onthefly_hca_coarsen_supermatrix(root, root, bfactory, (operator=='d' ? aca_eps/30.0 : aca_eps/4.0), acak, aca_eps, 1, HLIB_HCAII,lower,eta,0); + /* + blk = build_blockcluster(root, root, HLIB_MAXADMISSIBILITY, + HLIB_BLOCK_INHOMOGENEOUS, eta, nmin); + M = coarsen_hca_from_blockcluster(blk, bfactory, + (operator=='d' ? aca_eps/30.0 : aca_eps/4.0), + aca_eps, + 1, acak); + */ + break; case 5: blk = build_blockcluster(root, root, HLIB_MAXADMISSIBILITY, diff -N -r -c HLib-1.3old/Library/blockcluster.c HLib-1.3/Library/blockcluster.c *** HLib-1.3old/Library/blockcluster.c 2004-12-12 17:42:31.000000000 +0100 --- HLib-1.3/Library/blockcluster.c 2005-05-09 11:16:34.165570472 +0200 *************** *** 85,90 **** --- 85,106 ---- return nr; } + int + getnrleaves_blockcluster(pcblockcluster bc) + { + int i; + int nr; + + nr = 0; + if(bc->son != NULL) + for(i=0; iblock_rows * bc->block_cols; i++) + nr += getnrleaves_blockcluster(bc->son[i]); + else + nr = 1; + + return nr; + } + pblockcluster build_blockcluster(pccluster row, pccluster col, BlockAdmissiblityCriterion adm, diff -N -r -c HLib-1.3old/Library/blockcluster.h HLib-1.3/Library/blockcluster.h *** HLib-1.3old/Library/blockcluster.h 2004-12-12 17:42:31.000000000 +0100 --- HLib-1.3/Library/blockcluster.h 2005-05-09 11:16:34.165570472 +0200 *************** *** 62,67 **** --- 62,71 ---- int getnrnodes_blockcluster(pcblockcluster bc); + /* Returns the number of leaves in the block tree structure */ + int + getnrleaves_blockcluster(pcblockcluster bc); + /* Create a block cluster tree from the roots row/col of two cluster trees. row: Root row cluster col: Root column cluster diff -N -r -c HLib-1.3old/Library/hcoarsening.c HLib-1.3/Library/hcoarsening.c *** HLib-1.3old/Library/hcoarsening.c 2004-12-12 17:42:32.000000000 +0100 --- HLib-1.3/Library/hcoarsening.c 2005-05-09 11:16:45.189894520 +0200 *************** *** 785,788 **** --- 785,1006 ---- return sc; } + psupermatrix + coarsen_supermatrix_from_blockcluster(pcblockcluster bc, + double eps, + int diag, + prkmatrix (*buildrk)(pccluster row, + pccluster col, + void *data), + pfullmatrix (*buildfull)(pccluster row, + pccluster col, + void *data), + void *data) + { + psupermatrix s; + prkmatrix r; + pfullmatrix f; + size_t sz_current, sz_rk, sz_full; + int rows, cols, block_rows, block_cols; + int try_rk, try_full; + int ad; + int i, j; + + rows = bc->row->size; + cols = bc->col->size; + block_rows = bc->block_rows; + block_cols = bc->block_cols; + s = NULL; + try_rk = !diag; + try_full = 1; + + if(bc->son) { + assert(block_rows > 0); + assert(block_cols > 0); + + s = new_supermatrix(block_rows, block_cols, rows, cols, + NULL, NULL, NULL); + for(j=0; js[i+j*block_rows] = + coarsen_supermatrix_from_blockcluster(bc->son[i+j*block_rows], eps, + (i == j ? diag : 0), + buildrk, buildfull, data); + if(s->s[i+j*block_rows]->s && s->rows > 128 && s->cols > 128) + try_rk = 0; + } + } + else { + if(bc->type & HLIB_BLOCK_WEAKADM) { + r = buildrk(bc->row, bc->col, data); + s = new_supermatrix(1, 1, rows, cols, NULL, r, NULL); + } + else { + f = buildfull(bc->row, bc->col, data); + s = new_supermatrix(1, 1, rows, cols, NULL, NULL, f); + } + } + + assert(s != NULL); + + if(try_rk) { + sz_current = getsize_supermatrix(s); + + r = new_rkmatrix(0, s->rows, s->cols); + r->local_eps = eps; + + ad = rk_arithmetics_adaptive(); + rk_arithmetics_set_adaptive(1); + convertsuper2_rkmatrix(s, r); + rk_arithmetics_set_adaptive(ad); + sz_rk = getsize_rkmatrix(r); + + if(sz_rk < sz_current) { + if(s->s) { + for(j=0; jblock_cols; j++) + for(i=0; iblock_rows; i++) { + del_supermatrix(s->s[i+j*s->block_rows]); + s->s[i+j*s->block_rows] = NULL; + } + free(s->s); + s->s = NULL; + } + + if(s->r) { + del_rkmatrix(s->r); + s->r = NULL; + } + + if(s->f) { + del_fullmatrix(s->f); + s->f = NULL; + } + + s->r = r; + } + else + del_rkmatrix(r); + } + + if(try_full) { + sz_current = getsize_supermatrix(s); + + sz_full = (size_t) sizeof(double) * s->rows * s->cols; + + if(sz_full < sz_current && + sz_full / s->cols == (size_t) sizeof(double) * s->rows) { + f = new_fullmatrix(s->rows, s->cols); + + convertsuper2_fullmatrix(f, s); + + if(s->s) { + for(j=0; jblock_cols; j++) + for(i=0; iblock_rows; i++) { + del_supermatrix(s->s[i+j*s->block_rows]); + s->s[i+j*s->block_rows] = NULL; + } + free(s->s); + s->s = NULL; + } + + if(s->r) { + del_rkmatrix(s->r); + s->r = NULL; + } + + if(s->f) { + del_fullmatrix(s->f); + s->f = NULL; + } + + s->f = f; + } + } + + return s; + } + + typedef struct _hcacoarsen hcacoarsen; + typedef hcacoarsen *phcacoarsen; + + struct _hcacoarsen { + psurfacebemfactory bfactory; + double hca_eps; + int kmax; + }; + + static prkmatrix + hca_rkmatrix(pccluster row, pccluster col, void *data) + { + phcacoarsen hca = (phcacoarsen) data; + psurfacebemfactory bfactory = hca->bfactory; + prkmatrix r; + + r = new_rkmatrix(0, row->size, col->size); + hca_surfacebem(r, row, col, bfactory, hca->hca_eps, hca->kmax); + + if(bfactory->disp) { + bfactory->done++; + bfactory->disp(1, bfactory->done, bfactory->total, bfactory->dispdata); + } + + return r; + } + + static pfullmatrix + hca_fullmatrix(pccluster row, pccluster col, void *data) + { + phcacoarsen hca = (phcacoarsen) data; + psurfacebemfactory bfactory = hca->bfactory; + pfullmatrix f; + + f = new_fullmatrix(row->size, col->size); + bfactory->integrate_nearfield(row->start, row->size, + col->start, col->size, + bfactory, f->e, f->rows); + + if(bfactory->disp) { + bfactory->done++; + bfactory->disp(1, bfactory->done, bfactory->total, bfactory->dispdata); + } + + return f; + } + + psupermatrix + coarsen_hca_from_blockcluster(pcblockcluster bc, + psurfacebemfactory bfactory, + double hca_eps, double coarsen_eps, + int diag, int kmax) + { + hcacoarsen hca; + psupermatrix s; + + hca.bfactory = bfactory; + hca.hca_eps = hca_eps; + hca.kmax = kmax; + + if(bfactory->disp) { + bfactory->done = 0; + bfactory->total = getnrleaves_blockcluster(bc); + bfactory->disp(0, bfactory->done, bfactory->total, bfactory->dispdata); + } + + s = coarsen_supermatrix_from_blockcluster(bc, coarsen_eps, diag, + hca_rkmatrix, hca_fullmatrix, + &hca); + + if(bfactory->disp) { + bfactory->done = bfactory->total; + bfactory->disp(2, bfactory->done, bfactory->total, bfactory->dispdata); + } + + if(bfactory->acabuf) { + free(bfactory->acabuf); + bfactory->acabuf = NULL; + bfactory->acabufsize = 0; + } + + return s; + } diff -N -r -c HLib-1.3old/Library/hcoarsening.h HLib-1.3/Library/hcoarsening.h *** HLib-1.3old/Library/hcoarsening.h 2004-11-23 10:15:58.000000000 +0100 --- HLib-1.3/Library/hcoarsening.h 2005-05-09 11:16:45.189894520 +0200 *************** *** 74,77 **** --- 74,113 ---- void compresslow_supermatrix(psupermatrix s, double eps, int k); + /* Compute a coarsened H-matrix approximation of an H-matrix + implicitly defined by the block cluster tree bc and the + callback functions buildrk and buildfull + bc: Block cluster tree defining the structure of the original matrix + eps: Error tolerance for coarsening + diag: Non-zero if we are in a diagonal block and the rank is full + buildrk: Create an rkmatrix for an admissible block + buildfull: Create a fullmatrix for an inadmissible block + data: Auxiliary data for buildrk and buildfull */ + psupermatrix + coarsen_supermatrix_from_blockcluster(pcblockcluster bc, + double eps, + int diag, + prkmatrix (*buildrk)(pccluster row, + pccluster col, + void *data), + pfullmatrix (*buildfull)(pccluster row, + pccluster col, + void *data), + void *data); + + /* Compute a coarsened H-matrix approximation of an H-matrix + implicitly defined by the block cluster tree bc and the + hybrid cross approximation of the admissible leaves + bc: Block cluster tree defining the structure of the original matrix + bfactory: Information on the integral operator + hca_eps: Error tolerance for hybrid cross approximation + coarsen_eps: Error tolerance for coarsening + diag: Non-zero if we are in a diagonal block an no coarsening is allowed + kmax: Maximal rank */ + psupermatrix + coarsen_hca_from_blockcluster(pcblockcluster bc, + psurfacebemfactory bfactory, + double hca_eps, double coarsen_eps, + int diag, int kmax); + #endif