HSS Related Computations

H2Pack also provides functions for HSS representations of kernel matrices.

HSS construction in H2Pack is based on the proxy point method (like construction) and resembles a method called Recursive Skeletonization (RS). Compared to RS, HSS construction in H2Pack works for more general kernel matrices but requires an additional but relatively cheap ULV decomposition step if solves are required (as they usually are if HSS is being used).

Main functionality:

• HSS matrix representation construction for a kernel matrix
• HSS matrix-vector multiplication
• HSS matrix-matrix multiplication
• ULV decomposition of HSS matrix representation
• Direct solves involving the HSS matrix representation using its ULV decomposition

Note that the parallel implementation of HSS functionality in H2Pack is not as carefully load balanced as for functionality.

Basic Framework

The typical usage of these HSS-related functions is very similar to that of -related functions in Basic Usage.

1. Initialize an H2Pack data structure and specify that HSS mode is desired
   Functions: H2P_init() and H2P_run_HSS()

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

3. Compute the proxy points. The proxy points are generated the same way as for matrix construction
   Functions: H2P_generate_proxy_point_surface() and H2P_generate_proxy_point_ID_file()

4. Construct an HSS matrix representation
   Function: H2P_build()

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

6. Construct the ULV decomposition of the HSS matrix representation
   Functions: H2P_HSS_ULV_Cholesky_factorize() and H2P_HSS_ULV_LU_factorize()
   Remarks There are two different ULV decompositions:

   • LU-based ULV decomposition, works for any HSS matrix representation
   • Cholesky-based ULV decomposition, works only for H22 matrix representations of symmetric positive definite matrices

   For more details, see the reference:

   • Xia, J., Chandrasekaran, S., Gu, M., and Li, X. S. (2010). Fast algorithms for hierarchically semiseparable matrices. Numerical Linear Algebra with Applications, 17(6), 953-976

7. Compute the ULV decomposition and solve with this decomposition
   Functions: H2P_HSS_ULV_Cholesky_solve() and H2P_HSS_ULV_LU_solve()

8. Destroy the H2Pack data structure
   Function: H2Pack_destroy()

All the functions above except those in Steps 1, 6, and 7 are exactly the same as those used for -related computations (see Basic Usage). We thus only show steps 1, 6, and 7 below.

Detailed usage

Step 1: Initialize an H2Pack data structure

Call H2P_init() to initialize a H2Pack data structure. Then set the initialized data structure in HSS mode by calling H2P_run_HSS(). An example follows.

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);
H2P_run_HSS(h2pack);

Step 6: Construct the ULV decomposition of the HSS matrix representation

Call either H2P_HSS_ULV_Cholesky_factorize() or H2P_HSS_ULV_LU_factorize() to construct the corresponding ULV decompositions. The decomposition will be stored inside the H2Pack data structure (h2pack defined earlier`) The declarations of the two functions are:

// Input parameters:
//   h2pack : H2Pack structure with constructed HSS representation
//   shift  : Possible diagonal shift parameter k and ULV decomposition of (A + k * I) is constructed
// Output parameter:
//   h2pack : H2Pack structure with ULV LU factorization
void H2P_HSS_ULV_LU_factorize(H2Pack_p h2pack, const DTYPE shift);
void H2P_HSS_ULV_Cholesky_factorize(H2Pack_p h2pack, const DTYPE shift);

Step 7: Apply the ULV decomposition for direct solve of the HSS matrix representation

Call either H2P_HSS_ULV_Cholesky_solve() or H2P_HSS_ULV_LU_solve() to solve with a given right-hand side vector. Let's write the LU-based and Cholesky-based ULV decomposition as where in the Cholesky-based ULV decomposition. Each function can perform three different operations (chosen by the operation type):

• Solve
• Solve
• Solve

where the first two partial solves can be useful in preconditioning. The declarations of the two solve functions are

// Input parameters:
//   h2pack : H2Pack structure with ULV LU/Cholesky factorization
//   op     : Operation type, 1, 2, or 3
//   b      : Size >= h2pack->krnl_mat_size, right-hand side vector
// Output parameter:
//   x : Size >= h2pack->krnl_mat_size, solution vector. 
//       If op == 1, x satisfies L_{ULV} * x = b.
//       If op == 2, x satisfies U_{ULV} * x = b.
//       If op == 3, x satisfies A_{HSS} * x = b.
void H2P_HSS_ULV_LU_solve(H2Pack_p h2pack, const int op, const DTYPE *b, DTYPE *x);
void H2P_HSS_ULV_Cholesky_solve(H2Pack_p h2pack, const int op, const DTYPE *b, DTYPE *x);

Notes

• If an HSS matrix representation is very ill-conditioned, its ULV decomposition can have low accuracy.
• In comparison with matrix construction, HSS construction can be much more expensive, especially for 3D problems.

Examples

• ./example/example_HSS.c