HSS Related Computations

H2Pack also provides functions for HSS representations of kernel matrices.

HSS construction 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.

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 the HSS mode
   Functions: H2P_init() and H2P_run_HSS()

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()
   Remarks: These are the same set of proxy points for matrix construction

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

5. Multiply the HSS matrix representation by vectors and 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 representations
   • Cholesky-based ULV decomposition, only works for symmetric positive definite HSS matrix representations.

   For more details, see the reference:

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

7. Apply the ULV decomposition for direct solve of the HSS matrix representation
   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. 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