HSS Related Computations

H2Pack provides experimental functions for HSS-related operation for kernel matrices.

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
   Involved function: 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
   Involved function: H2P_partition_points()

3. Construct the proxy points that are used for low-rank interpolative decomposition (ID) in matrix construction.
   Involved function: 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
   Involved function: H2P_build()

5. Multiply the HSS matrix representation by vectors and vectors
   Involved function: H2P_matmul() and H2P_matvec()

6. Construct the ULV decomposition of the HSS matrix representation
   Involved function: 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
   Involved function: H2P_HSS_ULV_Cholesky_solve() and H2P_HSS_ULV_LU_solve()

8. Destroy the H2Pack data structure
   Involved function: H2Pack_destroy()

All the functions above except in Steps 1, 6, and 7 are exactly the same one used for -related computations.

Detailed usage

Step 1: Initialize an H2Pack data structure

Define a pointer to a H2Pack data structure and 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 is as 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() and 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 definition of the two functions is

// 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() and H2P_HSS_ULV_LU_solve() to solve a given right-hand-side vector. First formally write LU-based and Cholesky-based ULV decomposition as where in the Cholesky-based ULV decomposition. There are three operations for the direct solve:

• Solve
• Solve
• Solve

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

// 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 extremely ill-conditioned, its ULV decomposition can have low accuracy as well.
• In comparison with H2 matrix construction, HSS construction can be much more expensive, especially for 3D problems.

Examples

• ./example/example_HSS.c