HLIBpro
3.0

When constructing Hmatrices consisting of several differently defined subblocks, e.g., the block diagonal
\[A = \begin{pmatrix} A_0 & & \\ & A_1 & \\ & & A_2 \end{pmatrix} \]
out of given matrices \(A_0,A_1,A_2\), 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 supports two different methods.
The function
implements a direct way for assembling a new block matrix out of a given set of (sub) matrices.
Here, block_rows and block_cols define the number of block rows and block columns of the new matrix while submatrices stores pointers to the corresponding subblocks in column wise ordering. If the parameter copy is true
, a copy of each subblock is stored in the new matrix instead.
A basic example on the usage of assemble_block
for constructing a \(2 \times 2\) matrix is
For the function to work properly, the cluster trees of all subblocks have to be compatible, e.g., all subblocks in the same block column must be defined using the same column cluster tree. The same holds for the row cluster tree of all subblocks in a block row.
Upoin construction of \(A\), the index sets of all subblocks are adjusted to be numbered consecutively, e.g., A_00>row_is().last()+1 == A_01>row_is().first()
.
Subblocks may also be NULL
if the corresponding block does not exists, e.g., as in the above block diagonal example.
If all subblocks are of type THMatrix, the new matrix is also constructed as an Hmatrix. In particular, the new matrix will store matching permutations based upon the permutations of the given subblocks, which is important for importing/exporting data, e.g., RHS vectors. Furthermore, in this case all subblocks will loose the permutation information and will be reduced to standard TBlockMatrix objects.
The function assemble_block
is also available in the C interface:
The above example for a \(2 \times 2\) block matrix would look as
While assemble_block
provides a simple interface for direct block matrix assembly, it is limited in the supported options, e.g., symmetric matrices or reusing subblocks. Furthermore, the preconditions for the function to work properly are quite strict and hence, easily lead to errors if not take care off.
As an alternative, one may also set up a special matrix builder class, which handles the construction of each subblock of the block matrix directly. However, this process is more involved and starts with the partitioning of the index set to ensure the correct block structure in the resulting matrix.
Assuming a \(p \times p\) block matrix we need to define the indices belonging to each sub block of the corresponding cluster tree:
Here, the array first_part
holds the partition number of each each, while get_partition_of_index
returns this partition number for the given index. Since this is application dependend, it has to be implemented by the user.
Afterwards, the cluster tree builder TGeomPartCTBuilder may be used to use partitioning stored in first_part
to form the first level of the resulting cluster tree. For this, TGeomPartCTBuilder also needs a base builder for further partitioning:
The coordinate set coord
is assumed to be already defined by the user above.
The construction of the block cluster tree is as usual:
Once the block cluster tree is available, the final matrix can be constructed using a specialised matrix builder class. In the following a simple example of such a special builder is given, assuming a unified base matrix builder for all subblocks. Alternatively, special matrix builder for each subblock may be supplied, or matrix blocks may be constructed based on previously built matrix blocks.
The basic class definition for our matrix builder is:
The important function is build
, which actually constructs the final Hmatrix and needs to be overloaded to implement the special functionality needed for the blocked matrix.
For build
we start with a basic test and gather the actual block structure of our matrix:
Afterwards, the construction of the block matrix is performed one block after the other using the provided base matrix builder object:
nullptr
but either to use a zero rank TRkMatrix
or the special matrix class TZeroMatrix instead.Finally, the permutation of the matrix are set, which is independent on the structure of the matrix and is always performed during Hmatrix construction:
All other functions are abstract in TMatBuilder and hence, need to be implemented in TBlockMatBuilder. Since all actual matrix block constructions is delegated to _base_builder
, no matrix leafs are build in TBlockMatBuilder. Furthermore, ghost matrices are only needed in the distributed case and can be omitted in our example. For simplicity, we just throw an exception if one of these functions is called:
Finally, the remaining construction and the function matrix_format
are implemented as:
Here, the general case of an unsymmetric matrix is assumed.
Having defined TBlockMatBuilder
, the code for the actual matrix construction in the main function looks as:
Since the actual type of the base_builder
is application dependent, it is not defined here.