Quickstart

Htool-DDM can be used for different use cases:

Dependencies

Htool-DDM is a C++ header-only library, it requires:

  • C++ compiler with standard 14.

  • BLAS, to perform algebraic dense operations.

And

  • For DDM solvers: a MPI implementation, HPDDM and LAPACK.

  • For shared memory parallelism when using in-house \(\mathcal{H}\)-matrix algebra: OpenMP

  • for SVD (re)compression: LAPACK

  • CMake to generate a build-system for the examples and tests.

Installation

It is sufficient to include the include folder of this repository in your library. If you prefer, the following command copies the include folder your OS-specific include directory to make it more widely available: in the root of this repository on your system,

cmake -B build
cmake --build build --config Release --target install -- -j $(nproc)

Note

The install prefix can be modified using cmake ..  -DCMAKE_INSTALL_PREFIX:PATH=your/install/path instead of the first line.

Once installed, you can also use Htool-DDM as a cmake imported target as follows in your CMakeLists.txt:

find_package(Htool REQUIRED)
find_package(MPI REQUIRED)
find_package(BLAS REQUIRED)
# other optional packages

add_executable(your_target your_file.cpp)
target_link_libraries(your_target Htool::htool)

Geometric clustering

The hierarchical clustering of a geometry is contained in htool::Cluster and can be created using a htool::ClusterTreeBuilder. For example,

int number_points = 10000;
int spatial_dimension = 3;
int number_of_partitions = 2;
int number_of_children = 2;
std::vector<double> coordinates(number_points*spatial_dimension);
// coordinates = {...}
htool::ClusterTreeBuilder<double> cluster_tree_builder;
htool::Cluster<double> cluster = cluster_tree_builder.create_cluster_tree(number_points, spatial_dimension, coordinates.data(), number_of_children, number_of_partitions);

where

  • number_of_partitions defines the number of children at a level of the cluster tree called “partition”, which can be used to distribute data in a MPI context.

  • number_of_children defines the number of children for the other nodes that are not leaves.

Hierarchical compression

To use the in-house hierarchical compression htool::HMatrix, Htool-DDM will need two inputs:

  1. The underlying geometry of the kernel to be compressed for the Geometric clustering.

  2. A function to generate any coefficient of the matrix to be compressed.

Coefficient generator

A coefficient generator is defined by following the interface htool::VirtualGenerator, and in particular htool::VirtualGenerator::copy_submatrix(). For example:

class UserOperator: public htool::VirtualGenerator<double>{
    virtual void copy_submatrix(int M, int N, const int *rows, const int *cols, T *ptr) const override {
        for (int i = 0; i < M; i++) {
            for (int j = 0; j < N; j++) {
                ptr[i + M * j] = ...// A(rows[i], cols[j]) where A is the kernel to be compressed;
            }
        }
    }
};

Build a hierarchical matrix

To build a htool::HMatrix, a htool::HMatrixBuilder can be used:

int number_points = 10000;
int spatial_dimension = 3;
double epsilon = 1e-3;
double eta = 10;
char symmetry = 'N';
char uplo = 'N';
std::vector<double> coordinates(number_points*spatial_dimension);
// coordinates = {...}
htool::HMatrixBuilder<double> hmatrix_builder(number_points, spatial_dimension, coordinates.data());
htool::HMatrix<double> hmatrix = hmatrix_builder.build(UserOperator(),htool::HMatrixTreeBuilder<double>(epsilon, eta, symmetry, uplo));

Note

The geometric clustering is done within the constructor of htool::HMatrixBuilder. You can still access the resulting target and source clusters as public members of hmatrix_builder.

Use a hierarchical matrix

Basic linear algebra is provided for htool::HMatrix. Shared-memory parallelism is also supported via execution policy traits.

Table 1 Supported linear algebra

Operations

C++ function

Supported execution policy

\(\mathcal{H}\)-matrix vector product

htool::add_hmatrix_vector_product()

std::execution::seq, std::execution::par

\(\mathcal{H}\)-matrix matrix product

htool::add_hmatrix_matrix_product()

std::execution::seq, std::execution::par

\(\mathcal{H}\)-LU factorisation

htool::lu_factorization()

None

\(\mathcal{H}\)-LU solve

htool::lu_solve()

None

\(\mathcal{H}\)-Cholesky factorisation

htool::cholesky_factorization()

None

\(\mathcal{H}\)-Cholesky solve

htool::cholesky_solve()

None

Note

We try to have a similar API to BLAS/LAPACK or <linalg>. In particular, parallel version of linear algebra functions are available using execution policy traits. If none is given, it defaults to sequential operation.

DDM Solver

Distributed operator

To assemble a distributed operator, first a geometric clustering needs to be applied. The partition defined by the target geometric cluster tree is then used to define a row-wise distributed operator (see Row-wise distributed operator). htool::DefaultApproximationBuilder builds a row-wise distributed operator htool::DistributedOperator where each block of rows is compressed using a htool::HMatrix, which can be accessed as a public member.

htool::HMatrixTreeBuilder<double> hmatrix_tree_builder(epsilon, eta, symmetry, uplo)
htool::DefaultApproximationBuilder<double> distributed_operator_builder(UserOperator(), target_cluster, source_cluster, hmatrix_tree_builder,MPI_COMM_WORLD);
DistributedOperator<double> &distributed_operator = distributed_operator_builder.distributed_operator;

A htool::DistributedOperator object provides distributed products with matrices and vectors.

Linear solver

The iterative linear solve will be done via a htool::DDM object, which can be built using htool::DDMSolverBuilder. We give here an example for a simple block-Jacobi solver, in particular we use htool::DefaultApproximationBuilder::block_diagonal_hmatrix, which is a pointer to the \(\mathcal{H}\) matrix for the block diagonal associated with the local subdomain.

// Create DDM object
htool::DDMSolverBuilder<double> ddm_solver_build(distributed_operator,distributed_operator_builder.block_diagonal_hmatrix);
htool::DDM<double> solver = ddm_solver_build.solver;

// Prepare solver
HPDDM::Option &opt = *HPDDM::Option::get();
opt.parse("-hpddm_schwarz_method asm ");
ddm_with_overlap.facto_one_level();

// Solve
int mu = 1 // number of right-hand side
std::vector<double> x = ...// unknowns
std::vector<double> b = ...// right-hand sides
ddm_with_overlap.solve(b.data(), x.data(), mu);

Note

We rely on the external library HPDDM for the implementation of efficient iterative solvers. We refer to its cheatsheet listing the various possible options for the solver.