HLIBpro  2.9.1
Basic Matrix Algebra

All H-matrix arithmetic, with the exception of matrix-vector multiplication, is approximativ with respect to a predefined accuracy, e.g., fixed rank or fixed blockwise precision (see Accuracy Control for details).

# Matrix-Vector Multiplication

Matrix-vector multiplication in 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 is provided in the form of a vector update:

$y := \alpha \cdot M \cdot x + \beta y$

Here, $$M$$ is either $$A$$, $$A^T$$ or $$A^H$$ for a given matrix $$A$$. The individual form of $$A$$ is specified by a value of type matop_t:

• apply_normal/MATOP_NORM: $$A$$,
• apply_transposed/MATOP_TRANS: $$A^T$$,
• apply_adjoint/MATOP_ADJ: $$A^H$$.
Remarks
Since 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 is based on BLAS/LAPACK, using just the conjugate of a matrix, e.g. $$\bar A$$, is not supported.

Each matrix classes provides the method mul_vec for real and cmul_vec for complex valued factors $$\alpha,\beta$$ (see TMatrix), which perform the corresponding operations:

unique_ptr< TVector > x( ... ); // define some vector
unique_ptr< TVector > y( ... ); // define some vector
unique_ptr< TMatrix > A( ... ); // define some matrix
A->mul_vec( 1.0, x.get(), 0.0, y.get(), apply_normal );
A->mul_vec( complex( -2.0, 0.5 ), x.get(), complex( 1.0, 4.0 ), y.get(), apply_adjoint );

Additionally, the function mul_vec is available, which also handles the case of a distributed matrix (on distributed memory). To use it, the header file algebra/mul_vec.hh has to be included.

After including algebra/mat_add.hh into the corresponding source file, the function add is available for matrix addition. It implements the matrix update

$C := \alpha \cdot A + \beta \cdot C$

with matrices $$A, C$$. Here, $$\alpha$$ and $$\beta$$ are real valued. For the complex valued case, the corresponding function is called cadd:

unique_ptr< TMatrix > A( ... ); // define some matrix
unique_ptr< TMatrix > C( ... ); // define some matrix
TTruncAcc acc( 1e-6 );
-1.0, C.get(), acc );
cadd( complex( -0.5, 2.0 ), A.get(),
complex( 0.0, -1.0 ), C.get(), acc );
void add(const T_value alpha, const TMatrix *A, const T_value beta, TMatrix *C, const TTruncAcc &acc)
C ≔ α·A + β·C.

A crucial requirement for the matrix addition as well as for all matrix operations, is the compatibility of the corresponding index sets and cluster trees, i.e. if $$A \in K^{I \times J}$$ then also $$C \in K^{I \times J}$$ must hold and the cluster trees for $$I$$ and $$J$$ have to be equal for both matrices. Only the block cluster trees may differ, e.g. a different admissibility condition can be used for both matrices.

A special for of matrix addition is implemented in add_identity (cadd_identity), in which

$A := A + \lambda I$

is implemented:

void add_identity(TMatrix *A, const T_value lambda)
compute A ≔ A + λ·I

Currently, this is only implemented if $$A$$ has dense diagonal blocks, e.g. which is th case for standard admissibility. Furthermore, since the operation only affects the diagonal part, it is an exact addition without truncation.

# Matrix Multiplication

General matrix multiplication is provided by the module mat_mul, e.g. you have to include algebra/mat_mul.hh . It is again implemented in the form of a matrix update:

$C := \alpha \cdot A \cdot B + \beta \cdot C$

Both, $$A$$ and $$B$$ may furthermore be transposed or conjugate transposed during multiplication, e.g. $$A$$ may be replaced by $$A^T$$ or $$A^H$$ and analogously for $$B$$.

The actual multiplication is implemented in the functions multiply (for real valued $$\alpha,\beta$$) and cmultiply (for complex valued $$\alpha,\beta$$).

unique_ptr< TMatrix > A( ... ); // define some matrix
unique_ptr< TMatrix > B( ... ); // define some matrix
unique_ptr< TMatrix > C( ... ); // define some matrix
TTruncAcc acc( 1e-6 );
multiply( 2.0, apply_normal, A.get(), apply_transposed, B.get(),
-1.0, C.get(), acc );
void multiply(const T_value alpha, const matop_t op_A, const TMatrix *A, const matop_t op_B, const TMatrix *B, const T_value beta, TMatrix *C, const TTruncAcc &acc, TProgressBar *progress=NULL)
compute C ≔ β·C + α·op(A)·op(B)

Since matrix multiplication is computationally expensive, a progress meter may be provided to show the current state of the operation:

TConsoleProgressBar progress( cout );
multiply( 2.0, apply_normal, A.get(), apply_transposed, B.get(),
-1.0, C.get(), acc,
& progress );

As for the matrix addition, the index sets and cluster trees of the involved matrices have to be compatible, otherwise an exception is thrown during multiplication.