Class hmatrix

As part of castor, the hmatrix class provides a fully abstract interface for the manipulation of hierarchical matrices 1. The class itself and the related functions can be found in Class hmatrix. It uses another class Class bintree dedicated to binary space partitioning.

This is work-in-progress! The descriptions of the different features and examples will be added to this section progressively.

1

19. Börm, L. Grasedyck and W. Hackbusch, Hierarchical Matrices, 2015 (technical report)

template<typename T>
class castor::hmatrix

Hierarchical matrices.

Array with element of type T, constructed using X and Y nodes of type S, stored in hierarchical format using rank default compression.

Values are stored in matrix<T>, associated to hierarchical indices in matrix<std::size_t> of the same dimensions of the values. Recursion is constructed using std::vector<hmatrix<T>> object.

hmatrix<T> can be constructed using dense matrix, sparse matrix or lambda function defining block matrix values.

Public Functions

inline hmatrix()

hmatrix(std::size_t m, std::size_t n, double tol, T v = 0)

Single value constructor, fill leaf with tensor product. No tree is needed.

hmatrix(matrix<T> const &A, matrix<T> const &B, double tol)

Low-rank constructor, fill leaf with tensor product. No tree is needed.

template<typename S>
hmatrix(matrix<S> const &X, matrix<S> const &Y, double tol, matrix<T> const &M)

Dense matrix contructor using nodes, fill close interactions with dense matrix and far with compressed. Parallelized leaves constructor.

template<typename S>
hmatrix(matrix<S> const &X, matrix<S> const &Y, double tol, smatrix<T> const &Ms)

Sparse matrix contructor using nodes, fill non-zeros leaves with full matrix and empty leaves with a zero tensor product. Recursive constructor.

template<typename S>
hmatrix(matrix<S> const &X, matrix<S> const &Y, double tol, std::function<matrix<T>(matrix<std::size_t>, matrix<std::size_t>)> const &fct)

Lambda function contructor using nodes, fill close interactions with dense matrix and far with compressed. Parallelized leaves constructor.

template<typename S>
hmatrix(bintree<S> const &Xb, bintree<S> const &Yb, double tol)

Empty contructor using tree, define far leaves as compressed and close leaves as dense, fill with nothing.

template<typename S>
hmatrix(bintree<S> const &Xb, bintree<S> const &Yb, double tol, smatrix<T> const &Ms)

Sparse matrix contructor using tree, fill non-zeros leaves with full matrix and empty leaves with low-rank representation of a zeros block.

template<typename S>
hmatrix(bintree<S> const &Xb, bintree<S> const &Yb, smatrix<T> const &Ml, hmatrix<T> const &Mh, smatrix<T> const &Mr)

Projector constructor : Ml * Mh * Mr.

void full2lowrank()

Convert full leaf to low-rank using ACA compression.

void fusion()

Low-rank fusion if possible if and only if hierarchical leaf has four low-rank blocks.

void hfull(matrix<T> &M, matrix<std::size_t> I = {}, matrix<std::size_t> J = {}) const

In-place conversion to dense matrix : Ah -> M.

void hfull2(matrix<T> &M, matrix<std::size_t> idx, matrix<std::size_t> jdx, matrix<std::size_t> I = {}, matrix<std::size_t> J = {}) const

In-place conversion to dense sub-matrix : Ah(I,J) -> M.

void hinv()

In-place hierarchical matrix inversion : Ah -> Ah^(-1).

void hllowsolve(hmatrix<T> const &Lh)

In-place hmatrix-hmatrix lower-triangular left resolution : Ah -> Lh^(-1) * Ah.

void hllowsolve(matrix<T> &B, matrix<std::size_t> I = {}) const

In-place hmatrix-matrix lower-triangular left resolution : B -> Lh^(-1) * B.

void hlmtimes(matrix<T> const &B, matrix<T> &C, matrix<std::size_t> I = {}, matrix<std::size_t> J = {}) const

In-place hmatrix-matrix product : C -> C + Ah * B.

void hlu(hmatrix<T> &U)

In-place hmatrix-hmatrix lower-upper factorization : Ah -> mtimes(Ah,Uh).

void hlupsolve(matrix<T> &B, matrix<std::size_t> I = {}) const

In-place hmatrix-matrix upper-triangular left resolution : B -> Uh^(-1) * B.

void hrmtimes(matrix<T> const &A, matrix<T> &C, matrix<std::size_t> I = {}, matrix<std::size_t> J = {}) const

In-place matrix-hmatrix product : C -> C + A * Bh.

void hrupsolve(hmatrix<T> const &Uh)

In-place hmatrix-hmatrix upper-triangular right resolution : Bh -> Bh * Uh^(-1).

void hrupsolve(matrix<T> &B, matrix<std::size_t> J = {}) const

In-place matrix-hmatrix upper-triangular right resolution : B -> B * Uh^(-1).

void htranspose()

In-place hmatrix transposition : Bh -> Bh^t.

void leafptr(std::vector<hmatrix<T>*> &ptr, std::vector<matrix<std::size_t>> &idx, std::vector<matrix<std::size_t>> &jdx, matrix<std::size_t> I = {}, matrix<std::size_t> J = {})

Extract leaves pointers to compute externaly leaves data.

void recompress(double tol)

Recompression of dense and low-rank leaves with fusion.

void spydata(smatrix<logical> &Sf, smatrix<logical> &Sc, matrix<std::size_t> I = {}, matrix<std::size_t> J = {}) const

Extract sparse matrix Sf representing full leaves and Sc compressed ones.

void stat(matrix<std::size_t> &info) const

Extract hierarchical storage statistics, with : Leaf : #hierarchical #compressed #full Rate : #values/numel * 100 tol : accuracy

void tgeabhm(T alpha, matrix<T> const &A, matrix<T> const &B, T beta, matrix<std::size_t> I = {}, matrix<std::size_t> J = {})

In-place hmatrix product with low-rank matrix representation : Ch -> alpha*(A*B) + beta*Ch.

void tgehmhm(T alpha, hmatrix<T> const &Ah, hmatrix<T> const &Bh, T beta)

In-place hmatrix product : Ch -> alpha*(Ah*Bh) + beta*Ch.

hmatrix<T> &operator+=(T b)

In-place addition of hmatrix or scalar : Ah -> Ah + b, Ah -> Ah + Bh.

hmatrix<T> &operator*=(T b)

In-place hmatrix scalar multiplication : Ah -> Ah*b.

hmatrix<T> &operator+=(hmatrix<T> const &Bh)

matrix<T> operator()(matrix<std::size_t> const &I, matrix<std::size_t> const &J) const

Access to sub matrix with dense conversion.

inline hmatrix<T> const &getChildren(int i) const

inline matrix<std::size_t> const &getColumn(int i) const

inline matrix<T> const &getData(int i) const

inline matrix<std::size_t> const &getRow(int i) const

inline std::size_t getSize(int i) const

inline double getTolerance() const

inline int getType() const

inline std::vector<hmatrix<T>> &setChildren()

inline hmatrix<T> &setChildren(int i)

inline matrix<std::size_t> &setColumn(int i)

inline matrix<T> &setData(int i)

inline matrix<std::size_t> &setRow(int i)

inline matrix<std::size_t> &setSize()

inline double &setTolerance()

inline int &setType()