Given a sparse matrix \(A\) and a right hand hand side \(b\), the solution \(x\) of
\[Ax = b\]
is sought, whereby H-LU factorisation shall be used to improve the convergence of an iterative solver. No geometrical data is assumed and hence, algebraic clustering is performed for converting \(A\) into an H-matrix. Furthermore, nested dissection is applied to improve the efficiency of the H-LU factorisation.
Loading the Sparse Matrix
The program starts with the inclusion of the needed header files and the definition of the error handling macro.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <hlib-c.h>
#define CHECK_INFO { if ( info != HLIB_NO_ERROR ) \
{ char buf[1024]; hlib_error_desc( buf, 1024 ); \
printf( "\n%s\n\n", buf ); exit(1); } }
int
main ( int argc, char ** argv )
{
The next step is the definition of the sparse matrix, which shall be read from a file. 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 supports various foreign matrix formats, e.g. Matlab, SAMG, Harwell-Boeing or Matrixmarket. In this example, the matrix is stored in Harwell/Boeing format format:
hlib_matrix_t S = NULL;
HLIB_FNDECL hlib_matrix_t hlib_load_matrix(const char *filename, int *info)
The function hlib_load_matrix performs autodetection of the file format and is therefore suited for all supported formats.
Next, some properties of the matrix are queried:
int is_sym = 0;
printf( " symmetric = %d\n", is_sym );
HLIB_FNDECL int hlib_matrix_is_sym(const hlib_matrix_t A, int *info)
HLIB_FNDECL int hlib_matrix_is_herm(const hlib_matrix_t A, int *info)
HLIB_FNDECL size_t hlib_matrix_rows(const hlib_matrix_t A, int *info)
HLIB_FNDECL size_t hlib_matrix_cols(const hlib_matrix_t A, int *info)
We also want to read the right-hand side of the equation system, which is also expected to be stored in Harwell/Boeing format:
hlib_vector_t b;
HLIB_FNDECL hlib_vector_t hlib_load_vector(const char *filename, int *info)
Conversion into H-Format
Having a right hand side, the equation system above could now be solved using an iterative scheme. Unfortunately, most sparse matrices need some form of preconditioning for an acceptable convergence behaviour. The preconditioner shall be in the form of a H-LU factorisation of \(A\) and hence, of S. For this, we need a cluster tree and a block cluster tree.
Since, by assumption, no geometrical data is available, the clustering is performed purely algebraically by using the connectivity relation of the indices stored in the sparse matrix itself. Furthermore, nested dissection is to be used, which needs some special clustering.
Algebraic clustering uses graph partitioning as the underlying technique. 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 implements two graph partitioning algorithms (HLIB_ALG_BFS and HLIB_ALG_ML) and supports external algorithms, e.g. METIS and Scotch (HLIB_ALG_METIS and HLIB_ALG_SCOTCH). Here we use the BFS-based partitioning strategy together with nested dissection:
hlib_clustertree_t ct;
HLIB_FNDECL hlib_clustertree_t hlib_clt_build_alg_nd(const hlib_matrix_t S, const hlib_alg_t algtype, const unsigned int nmin, int *info)
@ HLIB_ALG_BFS
Definition: hlib-c.h:131
Having the cluster tree, the block cluster tree can be computed:
hlib_admcond_t adm;
hlib_blockclustertree_t bct;
& info ); CHECK_INFO;
HLIB_FNDECL hlib_permutation_t hlib_clt_perm_i2e(hlib_clustertree_t ct, int *info)
HLIB_FNDECL hlib_blockclustertree_t hlib_bct_build(const hlib_clustertree_t rowct, const hlib_clustertree_t colct, const hlib_admcond_t ac, int *info)
HLIB_FNDECL hlib_admcond_t hlib_admcond_alg(const hlib_adm_t crit, const double eta, const hlib_matrix_t S, const hlib_permutation_t row_perm, const hlib_permutation_t col_perm, int *info)
@ HLIB_ADM_WEAK
Definition: hlib-c.h:141
Here, the weak admissibility condition for graphs (HLIB_ADM_WEAK) is used. Since the block clusters to test are based on the internal ordering of the H-matrix but the sparse matrix is defined in the external ordering, the corresponding mappings between both is needed by the admissibility condition.
The final step is the conversion of the sparse matrix into H-format. The accuracy can be usually neglected (set to some arbitrary value), since low rank blocks are usually empty (clusters of positive distance usually have no overlapping basis functions).
hlib_matrix_t A;
HLIB_FNDECL hlib_acc_t hlib_acc_fixed_eps(const hlib_real_t eps)
HLIB_FNDECL hlib_matrix_t hlib_matrix_build_sparse(const hlib_blockclustertree_t bct, const hlib_matrix_t S, const hlib_acc_t acc, int *info)
H-LU factorisation
The factorisation is performed in place, i.e. the data of A is overwritten by the corresponding factores. This implies, that A is actually not a matrix object after the factorisation since the stored data depends on special handling. For this, linear operators are used with special versions for each factorisation algorithm to permit evaluation of the factorised matrix or its inverse.
hlib_linearoperator_t LU;
HLIB_FNDECL hlib_linearoperator_t hlib_matrix_factorise_inv(hlib_matrix_t A, const hlib_acc_t acc, int *info)
The factorisation is performed up to a block-wise accuracy of \(10^{-4}\). Depending on the given matrix this has to be changed to obtain a direct solver (single matrix-vector mult.) or a good preconditioner (small number of iterations).
Some properties of the factorisation are printed next, most notably the inversion error \(\|I-(LU)^{-1}S\|_2\) and the condition estimate \(\|LU\mathbf{1}\|_2\). However, to compare \((LU)^{-1}\) with \(S\), both need to have the same ordering of the indices. Instead of reordering \(S\), the factorised matrix may be wrapped in an object representing a permuted matrix:
hlib_linearoperator_t PLU;
LU,
& info ); CHECK_INFO;
HLIB_FNDECL hlib_permutation_t hlib_clt_perm_e2i(hlib_clustertree_t ct, int *info)
HLIB_FNDECL hlib_linearoperator_t hlib_perm_linearoperator(const hlib_permutation_t P, const hlib_linearoperator_t A, const hlib_permutation_t R, int *info)
Here, an operator was built, which first reorders w.r.t. the H-matrix, then applies the factorised matrix and finally reorders back.
Now, the inversion error maybe computed:
printf( " inversion error = %.4e\n",
HLIB_FNDECL hlib_real_t hlib_linearoperator_norm_inv_approx(const hlib_linearoperator_t A, const hlib_linearoperator_t B, int *info)
Solving the Equation System
Finally, the equation system can be solved using the inverse of the matrix factorisation as a preconditioner. The solution vector is created using the hlib_col_vector function of S, which returns a vector defined over the column index set of S.
hlib_vector_t x;
hlib_solver_t solver;
hlib_solve_info_t sinfo;
& sinfo, & info ); CHECK_INFO;
if ( solve_info.converged )
printf( " converged in %u steps with rate %.2e, |r| = %.2e\n",
solve_info.steps, solve_info.conv_rate, solve_info.res_norm );
else
printf( " not converged in %u steps\n", solve_info.steps );
HLIB_FNDECL void hlib_solver_initialise_start_value(hlib_solver_t solver, const int flag, int *info)
HLIB_FNDECL hlib_solver_t hlib_solver_auto(int *info)
HLIB_FNDECL void hlib_solver_solve(const hlib_solver_t solver, const hlib_linearoperator_t A, hlib_vector_t x, const hlib_vector_t b, const hlib_linearoperator_t W, hlib_solve_info_t *solve_info, int *info)
HLIB_FNDECL hlib_vector_t hlib_matrix_col_vector(const hlib_matrix_t A, int *info)
The iterative solver is automatically chosen depending on S and the factorisation (MINRES for symmetric matrices and GMRES for unsymmetric matrices). The start value of the iteration should be initialised by the solver instead of using the data provided by the user (no guess/data available).
The program is finished with the finalisation of 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 and the closing of the try block:
}
HLIB_FNDECL void hlib_matrix_free(hlib_matrix_t A, int *info)
HLIB_FNDECL void hlib_vector_free(hlib_vector_t v, int *info)
HLIB_FNDECL void hlib_clt_free(hlib_clustertree_t ct, int *info)
HLIB_FNDECL void hlib_bct_free(hlib_blockclustertree_t bct, int *info)
HLIB_FNDECL void hlib_done(int *info)
HLIB_FNDECL void hlib_solver_free(hlib_solver_t solver, int *info)
HLIB_FNDECL void hlib_admcond_free(const hlib_admcond_t ac, int *info)
HLIB_FNDECL void hlib_linearoperator_free(hlib_linearoperator_t A, int *info)
The Plain Program
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <hlib-c.h>
#define CHECK_INFO { if ( info != HLIB_NO_ERROR ) \
{ char buf[1024]; hlib_error_desc( buf, 1024 ); \
printf( "\n%s\n\n", buf ); exit(1); } }
int
main ( int argc, char ** argv )
{
hlib_matrix_t S = NULL;
int is_sym = 0;
hlib_vector_t b;
hlib_clustertree_t ct;
hlib_admcond_t adm;
hlib_blockclustertree_t bct;
hlib_matrix_t A;
hlib_linearoperator_t LU;
hlib_linearoperator_t PLU;
hlib_vector_t x;
hlib_solver_t solver;
hlib_solve_info_t sinfo;
printf( " symmetric = %d\n", is_sym );
& info ); CHECK_INFO;
LU,
& info ); CHECK_INFO;
printf( " inversion error = %.4e\n",
& sinfo, & info ); CHECK_INFO;
if ( solve_info.converged )
printf( " converged in %u steps with rate %.2e, |r| = %.2e\n",
solve_info.steps, solve_info.conv_rate, solve_info.res_norm );
else
printf( " not converged in %u steps\n", solve_info.steps );
}
1.9.1