- Generated by
1.8.20
|
HLIBpro
2.8.1
|
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.
The program starts with the inclusion of the needed header files and the definition of the error handling macro.
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:
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:
We also want to read the right-hand side of the equation system, which is also expected to be stored in Harwell/Boeing 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:
Having the cluster tree, the block cluster tree can be computed:
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).
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.
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:
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:
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.
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:
1.8.20