HLIBpro
1.2

Table of contents 
Given a sparse matrix and a right hand hand side , the solution of
is sought, whereby ๐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 into an ๐matrix. Furthermore, nested dissection is applied to improve the efficiency of the ๐LU factorisation.
The program starts with the inclusion of the needed header files and the initialisation of ๐๐ซ๐จ๐ก๐๐๐.
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, HarwellBoeing or Matrixmarket. In this example, the matrix is stored in Matlab format:
The class TAutoMatrixIO performs autodetection of the file format and is therefore suited for all supported formats.
By using autoptr, the corresponding object is automatically deleted when the current block is left, thereby reducing the danger of memory leaks. This technique is extensively used in ๐๐ซ๐จ๐ก๐๐๐.
Since any kind of matrix may be stored in the given file, e.g. a dense matrix, but in the following sparse matrices are expected, the type of the matrix is tested:
The macro IS_TYPE
returns true
, if the corresponding object is of the given type. This is an example of the runtime type information system (RTTI) in ๐๐ซ๐จ๐ก๐๐๐ and works similar to dynamic casts in C++. The advantage of the ๐๐ซ๐จ๐ก๐๐๐ RTTI is, that further information can be gathered, e.g. the type name through the function typestr()
, independent of the corresponding C++ compiler support.
Finally, since we know by now that M
contains a sparse matrix, we cast the pointer for easier use in the following. The macro ptrcast
is just an abbreviation of a C++ cast and was introduced in ๐๐ซ๐จ๐ก๐๐๐ for ease of use and for debugging purposes (ptrcast
is different depending on compilation options).
We now make use of the TSparseMatrix pointer S
and may output some information about it:
This includes the dimension of the matrix, i.e. the number of rows and columns, as well as the number of non zero elements. Furthermore, the value type, real of complex, and matrix properties are printed. Finally, the storage size and the Frobenius norm is computed.
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 ๐LU factorisation of 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 (TBFSAlgPartStart and TMLAlgPartStrat) and supports external algorithms, e.g. METIS and Scotch. Here we use the multilevel partitioning strategy implemented by TMLAlgPartStrat.
Nested dissection clustering is based on a standard bissection clustering algorithm, which therefore has to be defined explicitly in the form of TAlgCTBuilder, which on the other hand uses the graph partitioning strategy. Finally, TNDAlgCTBuilder computes the cluster tree for the index set defined by S:
Having the cluster tree, the block cluster tree can be computed:
Here, the weak admissibility condition for graphs, implemented by TWeakAlgAdmCond, is used. Since the block clusters to test are based on the internal ordering of the ๐matrix but the sparse matrix is defined in the external ordering, the corresponding mappings between both is needed by the admissibility condition.
The printed sparsity constant gives some hint about the complexity of the ๐matrix constructed over the block cluster tree, the smaller, the better.
The final conversion of the sparse matrix into ๐format is performed by TSparseMBuilder. Here, again the mappings between internal and external ordering are needed. The accuracy can be usually neglected, since low rank blocks are usually empty (clusters of positive distance usually have no overlapping basis functions).
Matrix construction may use multiple threads, so the first argument to the build
function is corresponding number to use. The mach_info
object holds some information about the machine, e.g. number of available processors.
After matrix construction, the size of the ๐matrix and its norm is printed together with the (relative) norm of the difference , indicating any approximation error.
The factorisation of A
may be computed differently, depending on the format of A
, e.g. if it is symmetric or not. In the latter case the standard factorisation is used, implemented in TLU, while in the symmetric case the (or for hermitian matrices) factorisation implemented in TLDL is computed.
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, inverse factorisation matrices are used which are available for each factorisation algorithm.
Some properties of the factorisation are printed next, most notably the inversion error and the condition estimate :
Finally, the equation system can be solved using the inverse of the matrix factorisation as a preconditioner. The right hand side is assumed to be available in Matlab format and is read in using TAutoVectorIO, which again performs file format auto detection. The solution vector is created using the col_vector
function of S
, which returns a vector defined over the column index set of S
.
The TAutoVectorIO class can also be used to write the solution vector to a file, again in Matlab format.
The program is finished with the finalisation of ๐๐ซ๐จ๐ก๐๐๐ and the closing of the try block: