Basic Usage

Input

To construct and apply an matrix representation of a kernel matrix K(X,X), the following inputs are needed:

• Point set X

  • pt_dim: dimension of the space that the points lie in (i.e., 1D, 2D, or 3D)
  • n_point: number of points
  • coord: coordinates of the points; an array of dimensions pt_dim * n_point

• Kernel function K(x,y)

  • krnl_eval: a function that takes two sets of points X₀, Y₀ and outputs kernel block K(X₀,Y₀)
  • krnl_param: optional parameters for the kernel function
  • krnl_dim: the dimension of kernel function, i.e., scalar kernels like the Gaussian have krnl_dim = 1 and 3D tensor kernels like Stokes have krnl_dim = 3

  See Using and Writing Kernel Functions for how to define and use krnl_eval and krnl_param.

• Parameters for the low-rank approximation in matrix construction

  • rel_tol with compression mode QR_REL_NRM: relative error threshold for the low-rank approximations
  • abs_tol with compression mode QR_ABS_NRM: absolute error threshold for the low-rank approximations
  • rank with compression mode QR_RANK: a fixed rank for the low-rank approximations

  It is recommended to use rel_tol with QR_REL_NRM which generally gives matrix-vector multiplication with a relative error of approximately rel_tol.

Basic Framework

The typical usage of H2Pack involves the following ordered steps:

1. Initialize an H2Pack data structure
   Function: H2P_init()

2. Hierarchically partition the point set coord, which bisects each dimension and is associated with a 2-/4-/8-ary partition tree for 1D/2D/3D points
   Function: H2P_partition_points()

3. Construct the proxy points that are used for low-rank interpolative decomposition (ID) in matrix construction.
   Functions: H2P_generate_proxy_point_surface() and H2P_generate_proxy_point_ID_file()
   The proxy points are used for efficiently compressing matrix blocks in matrix construction. H2Pack provides two proxy point selection schemes. Please read Proxy Points and their Reuse for more details.

4. Construct an matrix representation
   Function: H2P_build()
   In an representation, there are many kernel blocks that can be either precomputed and stored or dynamically computed. This corresponds to two different modes in H2Pack: ahead-of-time (AOT) mode or just-in-time (JIT) mode. Please read Two Running Modes for H2Pack for more details.

5. Multiply the matrix representation by matrices or vectors
   Functions: H2P_matmul() and H2P_matvec()

6. Destroy the H2Pack data structure
   Function: H2P_destroy()

Detailed usage of H2Pack

Step 0: Add H2Pack header file and compile

To use H2Pack in a C/C++ program, add #include "H2Pack.h" to your source files and add -I<H2Pack path>/include to the compiler flags (see the example directory for examples). You can also include H2Pack_2D_kernels.h or H2Pack_3D_kernels.h to use the optimized kernels provided by H2Pack (see Using and Writing Kernel Functions for more information).

Step 1: Initialize an H2Pack data structure

Define a pointer to an H2Pack data structure and call H2P_init() to initialize the H2Pack data structure. For example, use the following code to initialize an H2Pack data structure for a scalar kernel function for 3D points and set the relative error threshold to be .

int pt_dim = 3, krnl_dim = 1;
DTYPE rel_tol = 1e-6;
kernel_eval_fptr krnl_eval = Coulomb_3D;
void *krnl_param = NULL;
H2Pack_p h2pack;
H2P_init(&h2pack, pt_dim, krnl_dim, QR_REL_NRM, &rel_tol);

Note:DTYPE is a macro defined in include/H2Pack_config.h. DTYPE controls the floating point data type used in H2Pack, and by default is defined as double.

The declaration of H2P_init() is:

// Input parameters:
//   pt_dim        : Dimension of point coordinate
//   krnl_dim      : Dimension of kernel function output
//   QR_stop_type  : Partial QR stop criteria: QR_RANK, QR_REL_NRM, or QR_ABS_NRM
//   QR_stop_param : Pointer to QR stop parameter
// Output parameter:
//   h2pack_ : Initialized H2Pack data structure
void H2P_init(
    H2Pack_p *h2pack_, const int pt_dim, const int krnl_dim, 
    const int QR_stop_type, void *QR_stop_param
);

Step 2: Hierarchically partition the point set

First call H2P_calc_enclosing_box() to compute the root box that encloses all the points. If a file of precomputed proxy points is provided, H2P_calc_enclosing_box will adjust the root box in order to use the precomputed proxy points.

// Input parameters:
//   pt_dim  : Point dimension
//   n_point : Number of points
//   coord   : Size pt_dim-by-npt, each column is the coordinate of a point
//   fname   : Proxy point file name, can be NULL
// Output parameter:
//   enbox_ : Box that encloses all points in this node.
//            enbox[0 : pt_dim-1] are the corner with the smallest
//            x/y/z/... coordinates. enbox[pt_dim : 2*pt_dim-1] are 
//            the sizes of this box.
void H2P_calc_enclosing_box(const int pt_dim, const int n_point, const DTYPE *coord,
    const char *fname, DTYPE **enbox_);

Then call H2P_partition_points() to construct a hierarchical partitioning of the points.

// Input parameters:
//   n_point         : Number of points
//   coord           : Size pt_dim * n_point, each column is the coordinate of a point
//   max_leaf_points : Maximum point in a leaf node's box. If <= 0, will use 200 for
//                     2D points and 400 for 3D
//   max_leaf_size   : Maximum edge length of a leaf node's box. If == 0, max_leaf_points
//                     will be the only restriction.
//   h2pack          : H2Pack data structure initialized using H2P_init()
// Output parameter:
//   h2pack : H2Pack structure with point partitioning info
void H2P_partition_points(
    H2Pack_p h2pack, const int n_point, const DTYPE *coord, 
    int max_leaf_points, DTYPE max_leaf_size
);

Usually, you can use max_leaf_points = 0 and max_leaf_size = 0 for H2P_partition_points(). Here is an example code:

// The next two lines use the default values in H2Pack for maximum number
// of points in the leaf node and maximum edge length of leaf box. You can 
// change these values. 
int max_leaf_points = 0;
DTYPE max_leaf_size = 0.0;
char *pp_fname = "./PP_3DCoulomb_1e-6.dat";
H2P_calc_enclosing_box(pt_dim, n_point, coord, pp_fname, &h2pack->root_enbox);
H2P_partition_points(h2pack, n_point, coord, max_leaf_points, max_leaf_size);

After the hierarchical partitioning, a set of reordered point coordinates is stored in h2pack->coord. Array h2pack->coord_idx stores the mappings from h2pack->coord to the input point set coord: the i-th point in h2pack->coord is the h2pack->coord_idx[i]-th point in coord. In matrix-vector and matrix-matrix multiplications, the input and output vectors and matrices are automatically reordered using h2pack->coord_idx. External programs only need to use the ordering of the original coordinates coord.

Step 3: Select proxy points for matrix construction

H2Pack provides two methods to select proxy points for matrix construction.

The first method is a numerical selection scheme that works for general kernel functions. An example code:

// For storing proxy points on each level
H2P_dense_mat_p *pp;
char *pp_fname = "./PP_3DCoulomb_1e-6.dat";  // You can choose any file name you like
H2P_generate_proxy_point_ID_file(h2pack, krnl_param, krnl_eval, pp_fname, &pp);

The declaration of H2P_generate_proxy_point_ID_file() is:

// Input parameters:
//   h2pack     : Initialized H2Pack structure
//   krnl_param : Pointer to kernel function parameter array
//   krnl_eval  : Pointer to kernel matrix evaluation function
//   fname      : Proxy point file name, if == NULL or cannot find that file, compute all proxy points
// Output parameter:
//   pp_  : Array of proxy points for each level
void H2P_generate_proxy_point_ID_file(
    H2Pack_p h2pack, const void *krnl_param, kernel_eval_fptr krnl_eval, 
    const char *fname, H2P_dense_mat_p **pp_
);

This function will first try to load the proxy points from a file. If a file name is given but the file does not exist, then this function will compute the proxy points and write them into the file. For more details, please read Proxy Points and their Reuse).

Note: Using saved proxy points that were generated for a different kernel function or for a different relative error threshold could lead to inaccurate matrix construction.

The second method is to select proxy surface points that work for kernel functions from potential theory, such as the 2D/3D Laplace kernel, 3D Stokes kernel, and 3D RPY kernel. An example code:

// For storing proxy points on each level
H2P_dense_mat_p *pp;    
// The edge length of the root box enclosing all the points
DTYPE max_L = h2pack->root_enbox[pt_dim];
// Parameters used to decide the number of proxy surface points
// This is a heuristic method that works based on numerical tests
int num_pp, num_pp_dim = ceil(-log10(rel_tol));
if (num_pp_dim < 4 ) num_pp_dim = 4;
if (num_pp_dim > 10) num_pp_dim = 10;
if (pt_dim == 2) num_pp = 2 * pt_dim * num_pp_dim;
if (pt_dim == 3) num_pp = 2 * pt_dim * num_pp_dim * num_pp_dim;
H2P_generate_proxy_point_surface(
    pt_dim, pt_dim, num_pp, h2pack->max_level, 
    h2pack->min_adm_level, max_L, &pp
);

The declaration of H2P_generate_proxy_point_surface() is:

// Input parameters:
//   pt_dim      : Dimension of point coordinate
//   xpt_dim     : Dimension of extended point coordinate 
//                 (for points with radius info xpt_dim == pt_dim+1, otherwise xpt_dim == pt_dim)
//   min_npts    : Minimum number of proxy points on the box surface
//   max_level   : Maximum level (included) of a H2 tree, (root level == 0)
//   start_level : Minimum level that needs proxy points
//   max_L       : The size of the root node's enclosing box
// Output parameter:
//   pp_  : Array of proxy points for each level
void H2P_generate_proxy_point_surface(
    const int pt_dim, const int xpt_dim, const int min_npts, const int max_level, 
    const int start_level, DTYPE max_L, H2P_dense_mat_p **pp_
);

Note: If any of tol_norm and max_L changes, you should call H2P_generate_proxy_point_surface() to generate a new set of proxy points. Using saved proxy points in this case could lead to inaccurate matrix construction.

Note: To destroy a set of proxy points, use the following code:

for (int i = 0; i <= h2pack->max_level; i++)
    H2P_dense_mat_destroy(&pp[i]);

Step 4: Construct an matrix representation

Call H2P_build() to construct the matrix representation. The declaration of H2P_build() is:

// Input parameters:
//   h2pack          : H2Pack structure with point partitioning info
//   pp              : Array of proxy points for each level
//   BD_JIT          : 0 or 1, if B and D matrices are computed just-in-time in matvec
//   krnl_param      : Pointer to kernel function parameter array
//   krnl_eval       : Pointer to kernel matrix evaluation function
//   krnl_bimv       : Pointer to kernel matrix bi-matvec function
//   krnl_bimv_flops : FLOPs needed in kernel bi-matvec
// Output parameter:
//   h2pack : H2Pack structure with H2 representation matrices
void H2P_build(
    H2Pack_p h2pack, H2P_dense_mat_p *pp, const int BD_JIT, void *krnl_param, 
    kernel_eval_fptr krnl_eval, kernel_bimv_fptr krnl_bimv, const int krnl_bimv_flops
);

BD_JIT controls whether H2Pack uses ahead-of-time (AOT) mode or just-in-time (JIT) mode. If BD_JIT = 0, H2Pack computes and stores all kernel block components in the matrix representation. If BD_JIT = 1, H2Pack skips the computation of these blocks when constructing the matrix representation and computes these blocks dynamically when needed in the multiplication.

See Bi-Kernel Matvec (BKM) Functions for information on how to set krnl_bimv and krnl_bimv_flops for using the bi-kernel matvec optimization. If you wish to avoid using this optimization, you can set krnl_bimv = NULL and krnl_bimv_flops = 0.

Step 5: Multiply the matrix representation by vectors and matrices

Call H2P_matvec() to multiply the constructed matrix representation by a vector. The declaration of H2P_matvec() is:

// Input parameters:
//   h2pack : H2Pack structure with H2 representation
//   x      : Size h2pack->krnl_mat_size, input dense vector
// Output parameter:
//   y      : Size h2pack->krnl_mat_size, output dense vector
void H2P_matvec(H2Pack_p h2pack, const DTYPE *x, DTYPE *y);

To multiply the constructed matrix representation by multiple vectors (a dense matrix), you can call H2P_matvec multiple times (if these vectors are stored in column-major format) or call H2P_matmul(). The declaration of H2P_matmul() is:

// Input parameters:
//   h2pack : H2Pack structure with H2 representation
//   layout : CblasRowMajor/CblasColMajor if x & y are stored in row/column-major style
//   n_vec  : Number of column vectors in mat_x
//   mat_x  : Size >= h2pack->krnl_mat_size * ldx, input dense matrix, the leading 
//            h2pack->krnl_mat_size-by-n_vec part of mat_x will be used
//   ldx    : Leading dimension of mat_x, should >= n_vec if layout == CblasRowMajor,
//            should >= h2pack->krnl_mat_size if layout == CblasColMajor
//   ldy    : Leading dimension of mat_y, the same requirement of ldx
// Output parameter:
//   mat_y  : Size is the same as mat_x, output dense matrix, mat_y := A_{H2} * mat_x
void H2P_matmul(
    H2Pack_p h2pack, const CBLAS_LAYOUT layout, const int n_vec, 
    const DTYPE *mat_x, const int ldx, DTYPE *mat_y, const int ldy
);

For a performance comparison between calling H2P_matvec multiple times and calling H2P_matmul a single time, see Comparative tests on H2-matvec and H2-matmul.

Step 6: Destroy an H2Pack data structure

Call H2Pack_destroy(h2pack) to destroy the H2Pack data structure pointed by h2pack. (Before destroying it, you can call H2P_print_statistic(h2pack) to print diagnostic information contained in the H2Pack structure pointed by h2pack.) Here is the example code:

for (int i = 2; i <= h2pack->max_level; i++)
    H2P_dense_mat_destroy(&pp[i]);
H2P_destroy(&h2pack);  // Call it after destroying all proxy points

Examples

• ./example/example_H2.c
• ./example/example_H2_tensor.c