HLIBpro  2.9
Solving a Sparse Equation System
Remarks
This example corresponds to the C++ example "Solving a Sparse Equation System".

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.

Remarks
During the example, the definition of variables is done when needed to make the presentation more structured. This of course is not valid C code. In the plain program at the end of the presentation, all variable definitions are placed correctly at the start of code blocks.

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); } }
/*
* main function
*/
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;
S = hlib_load_matrix( "matrix.hb", & info ); CHECK_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;
is_sym = ( hlib_matrix_is_sym( S, & info ) || hlib_matrix_is_herm( S, & info ) ? 1 : 0 );
printf( " matrix has dimension %lu x %lu\n", hlib_matrix_rows( S, & info ),
hlib_matrix_cols( S, & info ) );
printf( " symmetric = %d\n", is_sym );

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;
b = hlib_load_vector( "rhs.hb", & info ); CHECK_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;
ct = hlib_clt_build_alg_nd( S, HLIB_ALG_BFS, 40, & info ); CHECK_INFO;

Having the cluster tree, the block cluster tree can be computed:

hlib_admcond_t adm;
hlib_blockclustertree_t bct;
hlib_clt_perm_i2e( ct, & info ),
hlib_clt_perm_i2e( ct, & info ),
& info ); CHECK_INFO;
bct = hlib_bct_build( ct, ct, adm, & info ); CHECK_INFO;

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;
A = hlib_matrix_build_sparse( bct, S, acc, & info ); CHECK_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_acc_t lu_acc = hlib_acc_fixed_eps( 1e-4 );
LU = hlib_matrix_factorise_inv( A, lu_acc, & info ); CHECK_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,
hlib_clt_perm_e2i( ct, & info ),
& info ); CHECK_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_linearoperator_norm_inv_approx( (hlib_linearoperator_t) S, PLU, & info ) ); CHECK_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;
x = hlib_matrix_col_vector( S, & info ); CHECK_INFO;
solver = hlib_solver_auto( & info ); CHECK_INFO;
hlib_solver_initialise_start_value( solver, 1, & info ); CHECK_INFO;
hlib_solver_solve( solver, (hlib_linearoperator_t) S, x, b, PLU,
& 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 );

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_solver_free( solver, & info ); CHECK_INFO;
hlib_vector_free( x, & info ); CHECK_INFO;
hlib_linearoperator_free( PLU, & info ); CHECK_INFO;
hlib_linearoperator_free( LU, & info ); CHECK_INFO;
hlib_vector_free( b, & info ); CHECK_INFO;
hlib_matrix_free( A, & info ); CHECK_INFO;
hlib_bct_free( bct, & info ); CHECK_INFO;
hlib_admcond_free( adm, & info ); CHECK_INFO;
hlib_clt_free( ct, & info ); CHECK_INFO;
hlib_matrix_free( S, & info ); CHECK_INFO;
hlib_done( & info ); CHECK_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); } }
/*
* main function
*/
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_acc_t lu_acc;
hlib_linearoperator_t PLU;
hlib_vector_t x;
hlib_solver_t solver;
hlib_solve_info_t sinfo;
S = hlib_load_matrix( "matrix.hb", & info ); CHECK_INFO;
is_sym = ( hlib_matrix_is_sym( S, & info ) || hlib_matrix_is_herm( S, & info ) ? 1 : 0 );
printf( " matrix has dimension %lu x %lu\n", hlib_matrix_rows( S, & info ),
hlib_matrix_cols( S, & info ) );
printf( " symmetric = %d\n", is_sym );
b = hlib_load_vector( "rhs.hb", & info ); CHECK_INFO;
ct = hlib_clt_build_alg_nd( S, HLIB_ALG_BFS, 40, & info ); CHECK_INFO;
hlib_clt_perm_i2e( ct, & info ),
hlib_clt_perm_i2e( ct, & info ),
& info ); CHECK_INFO;
bct = hlib_bct_build( ct, ct, adm, & info ); CHECK_INFO;
acc = hlib_acc_fixed_eps( 0.0 );
A = hlib_matrix_build_sparse( bct, S, acc, & info ); CHECK_INFO;
lu_acc = hlib_acc_fixed_eps( 1e-4 );
LU = hlib_matrix_factorise_inv( A, lu_acc, & info ); CHECK_INFO;
LU,
hlib_clt_perm_e2i( ct, & info ),
& info ); CHECK_INFO;
printf( " inversion error = %.4e\n",
hlib_linearoperator_norm_inv_approx( (hlib_linearoperator_t) S, PLU, & info ) ); CHECK_INFO;
x = hlib_matrix_col_vector( S, & info ); CHECK_INFO;
solver = hlib_solver_auto( & info ); CHECK_INFO;
hlib_solver_initialise_start_value( solver, 1, & info ); CHECK_INFO;
hlib_solver_solve( solver, (hlib_linearoperator_t) S, x, b, PLU,
& 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_solver_free( solver, & info ); CHECK_INFO;
hlib_vector_free( x, & info ); CHECK_INFO;
hlib_linearoperator_free( PLU, & info ); CHECK_INFO;
hlib_linearoperator_free( LU, & info ); CHECK_INFO;
hlib_vector_free( b, & info ); CHECK_INFO;
hlib_matrix_free( A, & info ); CHECK_INFO;
hlib_bct_free( bct, & info ); CHECK_INFO;
hlib_admcond_free( adm, & info ); CHECK_INFO;
hlib_clt_free( ct, & info ); CHECK_INFO;
hlib_matrix_free( S, & info ); CHECK_INFO;
hlib_done( & info ); CHECK_INFO;
}
hlib_acc_u
Definition: hlib-c.h:188
hlib_matrix_is_herm
HLIB_FNDECL int hlib_matrix_is_herm(const hlib_matrix_t A, int *info)
hlib_clt_build_alg_nd
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_matrix_build_sparse
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)
hlib_load_matrix
HLIB_FNDECL hlib_matrix_t hlib_load_matrix(const char *filename, int *info)
hlib_bct_build
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_matrix_factorise_inv
HLIB_FNDECL hlib_linearoperator_t hlib_matrix_factorise_inv(hlib_matrix_t A, const hlib_acc_t acc, int *info)
hlib_linearoperator_free
HLIB_FNDECL void hlib_linearoperator_free(hlib_linearoperator_t A, int *info)
hlib_solver_solve
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_clt_perm_e2i
HLIB_FNDECL hlib_permutation_t hlib_clt_perm_e2i(hlib_clustertree_t ct, int *info)
hlib_solver_free
HLIB_FNDECL void hlib_solver_free(hlib_solver_t solver, int *info)
HLIB_ADM_WEAK
@ HLIB_ADM_WEAK
Definition: hlib-c.h:141
hlib_solver_initialise_start_value
HLIB_FNDECL void hlib_solver_initialise_start_value(hlib_solver_t solver, const int flag, int *info)
hlib_admcond_alg
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_perm_linearoperator
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)
hlib_matrix_rows
HLIB_FNDECL size_t hlib_matrix_rows(const hlib_matrix_t A, int *info)
HLIB_ALG_BFS
@ HLIB_ALG_BFS
Definition: hlib-c.h:131
hlib_done
HLIB_FNDECL void hlib_done(int *info)
hlib_bct_free
HLIB_FNDECL void hlib_bct_free(hlib_blockclustertree_t bct, int *info)
hlib_linearoperator_norm_inv_approx
HLIB_FNDECL hlib_real_t hlib_linearoperator_norm_inv_approx(const hlib_linearoperator_t A, const hlib_linearoperator_t B, int *info)
hlib_clt_free
HLIB_FNDECL void hlib_clt_free(hlib_clustertree_t ct, int *info)
hlib_load_vector
HLIB_FNDECL hlib_vector_t hlib_load_vector(const char *filename, int *info)
hlib_admcond_free
HLIB_FNDECL void hlib_admcond_free(const hlib_admcond_t ac, int *info)
hlib_clt_perm_i2e
HLIB_FNDECL hlib_permutation_t hlib_clt_perm_i2e(hlib_clustertree_t ct, int *info)
hlib_matrix_cols
HLIB_FNDECL size_t hlib_matrix_cols(const hlib_matrix_t A, int *info)
hlib_matrix_col_vector
HLIB_FNDECL hlib_vector_t hlib_matrix_col_vector(const hlib_matrix_t A, int *info)
hlib_vector_free
HLIB_FNDECL void hlib_vector_free(hlib_vector_t v, int *info)
hlib_solver_auto
HLIB_FNDECL hlib_solver_t hlib_solver_auto(int *info)
hlib_matrix_free
HLIB_FNDECL void hlib_matrix_free(hlib_matrix_t A, int *info)
hlib_matrix_is_sym
HLIB_FNDECL int hlib_matrix_is_sym(const hlib_matrix_t A, int *info)
hlib_acc_fixed_eps
HLIB_FNDECL hlib_acc_t hlib_acc_fixed_eps(const hlib_real_t eps)