Basic Matrix Algebra

Table of Contents

• Matrix-Vector Multiplication
• Matrix Addition
• Matrix Multiplication

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: \( A \),
• apply_transposed: \( A^T \),
• apply_adjoint: \( 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< value_t > >  x( ... ); // define some vector
unique_ptr< TVector< value_t > >  y( ... ); // define some vector
unique_ptr< TMatrix< value_t > >  A( ... ); // define some matrix
 
A->mul_vec( 1.0, x.get(), 0.0, y.get(), apply_normal );

Additionally, the function mul_vec is available. To use it, the header file algebra/mul_vec.hh has to be included.

Matrix Addition

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< value_t > >  A( ... ); // define some matrix
unique_ptr< TMatrix< value_t > >  C( ... ); // define some matrix
auto                              acc = fixed_prec( 1e-6 );
 
add(  2.0, A.get(), -1.0, C.get(), acc );

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:

add_identity( A.get(), 2.0 );

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< value_t > >  A( ... ); // define some matrix
unique_ptr< TMatrix< value_t > >  B( ... ); // define some matrix
unique_ptr< TMatrix< value_t > >  C( ... ); // define some matrix
auto                              acc = fixed_prec( 1e-6 );
 
multiply(  2.0, apply_normal, A.get(), apply_transposed, B.get(), -1.0, C.get(), acc );

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.