HLIBpro  2.9.1
BLAS/LAPACK Interface

The foundation of HLIBpro is provided by the BLAS and LAPACK software packages. Especially the advanced functions of LAPACK, e.g. singular value decomposition, play a crucial role in the H-arithmetics. To simplify the usage of these functions and to reduce the potential for errors, an extra layer was implemented between the H-arithmetics and the BLAS/LAPACK functions. The BLAS module provides this layer.

All data classes in the BLAS namespace share special properties with respect to the memory management of the underlying data:

  • Each class is implemented as a template with the base datatype as the template parameter, e.g. double or complex.
  • Data is not necessarily contiguous inside a memory block, e.g. you may address each second element of a memory block. This is especially needed to address sub matrices, where each column starts with a certain offset with respect to the previous column.
  • You may not change the size of a vector or matrix, only assigning a vector or matrix of a different size.
  • In contrast to H-matrices and vectors, indices for BLAS vectors and matrices always start at 0.

Since the matrix and vector classes were implemented with BLAS/LAPACK in mind, some restrictions apply:

  • All data is stored column-wise.
  • No mix of datatypes is possible, e.g. multiplying double and complex matrices.

For the following code, an

using namespace BLAS;

is assumed.


Memory Management

Due to performance issues on multi-core systems the memory management of the BLAS classes was changed in HLIBpro. Please read the following carefully to prevent programming errors.

Both vector and matrix classes in the BLAS namespace store a pointer to some memory block. This memory may be either owned or referenced by the corresponding object. If the pointer only references some memory block owned by another object, no reference counting is done. Hence, if the owning object is freed, the reference will be invalid.

A major cause of such problems are temporary objects, e.g. some auxilliary vectors or matrices inside a function, which are assigned to some return value:

void compute_x ( Vector< double > & x ) {
Vector t;
...
x = t;
}

Here, x will by default only reference t. To make the new owner of the data in t, one has to use the move semantics in C++:

void compute_x ( Vector< double > & x ) {
Vector t;
...
x = std::move( t );
}

Now, t will be empty and x owns the data of t.

If the result is the return value of the function, the compiler will automatically choose the move semantics:

Vector< double > compute_x () {
Vector t;
...
return t;
}

Ranges

A range is an index set in the BLAS namespace.

Range r1( 0, 10 );

Here, an index set containing \(\left\{0,1,\ldots,10\right\}\) is defined. Please note, that the index interval is closed on both ends, e.g. the first and the last index is given.

In addition to a standard index set, you may also define a stride, e.g. a distance between elements in the index set:

Range r2( 10, 20, 2 );

This defines the index set \(\left\{10,12,14,\ldots,20\right\}\)

For convenience, a constructor for an index set containing just one element is also available:

Range r3( 21 );

yields \(\left\{21\right\}\).

Furthermore, the special range Range::all indicates, that the whole range shall be used:

Matrix< double > M( 10, 20 );
Matrix< double > M2( M, Range( 5, 10 ), Range::all );

Here, all columns of M are used in M2 but only the rows 5 to 10.


Vectors

Construction

Usually, a BLAS vector is intialised with the corresponding size as the constructor argument:

Vector< double > x( 10 );

creates a double vector of size 10. Leaving out the size argument, intialises the vector with size 0:

Vector< double > y;

Providing another vector in the constructor, creates a new vector pointing to the same data as the argument vector:

Vector< double > z( x );

Changing z also changes x!

With the usage of range objects, parts of vectors can be represented as another vector:

Vector< double > x1( 10 );
Vector< double > x2( x1, Range( 5, 10 ) );
Vector< double > x3( x2, Range( 0, 5, 2 ) );

Here, x2 addresses the second half of the coefficients of x1 whereas x3 contains each second coefficient of x2 and by that, the coefficients 5, 7 and 9 of x1.

By default, the policy for the copy constructors is to use a reference to the given data. You may change that, such that a copy of the data is created, e.g. a new memory block is allocated:

Vector< double > x3( x2, Range( 0, 5, 2 ), copy_value );

Now, x3 does not address any coefficients in x1 or x2.

Accessing Coefficients and Vector Properties

Single coefficients in vectors are accessed using the () operator:

Vector< double > x( 10 );
for ( idx_t i = 0; i < 10; ++i )
x( i ) = i;
for ( idx_t i = 0; i < 10; ++i )
std::cout << x( i ) << std::endl;

The length of a vector is accessible via the length() function:

for ( idx_t i = 0; i < x.length(); ++i )
x( i ) = i;

Algebra

Most of the BLAS vector functions are available in HLIBpro and some, which have no representation in the BLAS package.

An examples, involving filling, scaling vectors, forming linear combinations and computing dot products:

Vector< double > x( 10 );
Vector< double > y( 10 );
Vector< double > z( 10 );
fill( 2, x ); // x is fill with 2
fill( 3.5, y ); // y is fill with 3.5
copy( x, z ); // copy coefficients of x to z
add( 0.5, y, z ); // set z ≔ z + 0.5 x
std::cout << dot( x, z ) << std::endl; // compute x^H·z
void add(const T_value alpha, const TMatrix *A, const T_value beta, TMatrix *C, const TTruncAcc &acc)
C ≔ α·A + β·C.
std::enable_if_t< is_vector< T1 >::value &&is_vector< T2 >::value &&is_same_type< typename T1::value_t, typename T2::value_t >::value, typename T1::value_t > dot(const T1 &x, const T2 &y)
compute <x,y> = x^H · y
Definition: Algebra.hh:279
std::enable_if_t< is_vector< T2 >::value &&std::is_same< T1, typename T2::value_t >::value, void > fill(const T1 f, T2 &x)
fill vector with constant
Definition: Algebra.hh:71
std::enable_if_t< is_vector< T1 >::value &&is_vector< T2 >::value &&is_same_type< typename T1::value_t, typename T2::value_t >::value, void > copy(const T1 &x, T2 &y)
copy x into y
Definition: Algebra.hh:145

As mentioned before, you can not mix vectors of different value type, e.g.

Vector< double > x( 10 );
Vector< float > y( 10 );
copy( x, y ); // ERROR

since BLAS/LAPACK do not support such mixing

For a complete overview of the algebra functions refer to the BLAS module.


Matrices

Construction

Constructing matrices follows the same scheme as vectors, with the only difference being the additional dimension:

Matrix< double > M1( 10, 20 );
Matrix< complex > M2( 20, 20 );
Matrix< double > M3( M1 );
Matrix< double > M4;

Addressing subblocks of a given matrix is done by providing Range objects:

Matrix< double > M1( 20, 20 );
Matrix< double > M2( M1, Range( 10, 19 ), Range( 0, 19 ) );
Matrix< double > M3( M2, Range( 0, 9 ), Range( 0, 19, 2 ) );

Here, M2 references the last 10 rows of M1. M3 also shares these coefficients but further restricts the elements to each second column of M1. Hence

M1( 14, 10 ) = 2.0;
M2( 4, 10 ) = 3.0;
M3( 4, 5 ) = 4.0;

will affect the same memory position.

Remarks
This also introduces the coefficient access for matrices, again being the standard index operator.

The standard copy policy for matrices, i.e. copy_reference, can be changed by the last, optional parameter to the above constructor:

Matrix< double > M4( M2, Range( 0, 9 ), Range( 0, 19, 2 ), copy_value );

This will create a 10×10 matrix with coefficients copied from matrix M2.

Algebra

As for vectors, the arithmetic functionality is implemented as functions and not as special C++ operators for objects of type Matrix.

Matrix views

𝖧𝖫𝖨𝖡𝗉𝗋𝗈 provides different <emph>matrix view</emph> for accessing matrices, e.g. transposed and hermitian view. A view object will not change or copy the data of a matrix but changes for instance the way the data is access, e.g. interchanging the indices. For a given matrix

Matrix< double > M( 10, 10 );

The transposed view of M is defined by

transpose( M );
std::enable_if_t< is_matrix< T1 >::value, void > transpose(T1 &A)
transpose matrix A: A → A^T ASSUMPTION: A is square matrix
Definition: Algebra.hh:569

and the adjoint view by

adjoint( M );
matop_t adjoint(const matop_t op)
return adjoint of given matrix operation
Definition: matrix_view.hh:65

These views may be used for arithmetic functions, albeit not all functions, to simplify the usage and to enhance the performance of these functions.

Matrix-Vector Arithmetic

The most common form of matrix-vector arithmetic is the matrix vector multiplication in the general form

\[ y := \alpha A x + \beta y \]

:

Matrix< double > M( 20, 10 );
Vector< double > x( 10 );
Vector< double > y( 20 );
mulvec( 2.0, M, x, 1.0, y );
std::enable_if_t< is_matrix< T2 >::value &&is_vector< T3 >::value &&is_vector< T4 >::value &&is_same_type< T1, typename T2::value_t >::value &&is_same_type< T1, typename T3::value_t >::value &&is_same_type< T1, typename T4::value_t >::value, void > mulvec(const T1 alpha, const T2 &A, const T3 &x, const T1 beta, T4 &y)
compute y ≔ β·y + α·A·x
Definition: Algebra.hh:821

Using matrix views allows the multiplication with the transposed or adjoint matrix:

mulvec( 2.0, transposed(M), x, 1.0, y );
mulvec( 2.0, adjoint(M), x, 1.0, y );
matop_t transposed(const matop_t op)
return transposed of given matrix operation
Definition: matrix_view.hh:47

If M is a triangular matrix, the multiplication is performed in place, i.e. the single argument vector is modified:

\[ x := A \cdot x \]

The exact shape of the matrix is defined by additional arguments of type tri_type_t (lower_triangular or upper_triangular) and diag_type_t (unit_diag or general_diag). The second argument describes the diagonal which is either identical to 1 or not:

mulvec_tri( upper_triangular, general_diag, transposed(M), x );
std::enable_if_t< is_matrix< T1 >::value &&is_vector< T2 >::value &&is_same_type< typename T1::value_t, typename T2::value_t >::value, void > mulvec_tri(const tri_type_t shape, const diag_type_t diag, const T1 &A, T2 &x)
compute x ≔ M·x, where M is upper or lower triangular with unit or non-unit diagonal
Definition: Algebra.hh:950
Remarks
The shape of the matrix, e.g. lower or upper triangular, always refers to the original matrix and not the provided view!

Also available is the rank-1 update of a matrix

\[ A := A + \alpha a \cdot b^H \]

:

Matrix< double > M( 20, 30 );
Vector< double > a( 20 );
Vector< double > b( 30 );
add_r1( 0.5, a, b, M );
std::enable_if_t< is_vector< T2 >::value &&is_vector< T3 >::value &&is_matrix< T4 >::value &&is_same_type< T1, typename T2::value_t >::value &&is_same_type< T1, typename T3::value_t >::value &&is_same_type< T1, typename T4::value_t >::value, void > add_r1(const T1 alpha, const T2 &x, const T3 &y, T4 &A)
compute A ≔ A + α·x·y^H
Definition: Algebra.hh:709

For triangular matrices \(A\), solving equation systems \(Ax=b\) is supported via:

Matrix< double > A( 20, 20 );
Vector< double > x( 20 );
solve_tri( lower_triangular, unit_diag, adjoint(A), x )
std::enable_if_t< is_matrix< T1 >::value &&is_vector< T2 >::value &&is_same_type< typename T1::value_t, typename T2::value_t >::value, void > solve_tri(const tri_type_t uplo, const diag_type_t diag, const T1 &A, T2 &b)
solve A·x = b with known A and b; x overwrites b
Definition: Algebra.hh:1559

Here, upon input x contains \(b\)b and is overwritten with the result \(x\).

Matrix-Matrix Arithmetic

The BLAS namespace defines common functions for single matrices and between matrices. Single matrix functions are

Matrix< complex > M( 20, 20 );
fill( 3.14, M ); // fill M with 3.14
scale( 2.72, M ); // scale M by 2.72
conj( M ); // replace coefficients with conjugate
transpose( M ); // inplace transpose or
conj_transpose( M ); // conjugate transpose of M
std::enable_if_t< is_vector< T2 >::value &&is_same_type< T1, typename T2::value_t >::value, void > scale(const T1 f, T2 &x)
scale vector by constant
Definition: Algebra.hh:122
std::enable_if_t< is_vector< T1 >::value, void > conj(T1 &x)
fill vector with repeated function evaluation, x_i = f()
Definition: Algebra.hh:102
std::enable_if_t< is_matrix< T1 >::value, void > conj_transpose(T1 &A)
conjugate transpose matrix A: A → A^H ASSUMPTION: A is square matrix
Definition: Algebra.hh:592
Remarks
PLease note, that the last two functions will only work for square matrices!

Different matrix norms are available:

std::cout << normF( M ) << std::endl; // print Frobenius norm
std::cout << norm2( M ) << std::endl; // print spectral norm
std::cout << cond( M ) << std::endl; // print condition number
std::enable_if_t< is_matrix< T1 >::value, typename real_type< typename T1::value_t >::type_t > cond(const T1 &M)
return condition of M
std::enable_if_t< is_matrix< T1 >::value, typename real_type< typename T1::value_t >::type_t > normF(const T1 &M)
return Frobenius norm of M
Definition: Algebra.hh:1315
std::enable_if< is_vector< T1 >::value, typename real_type< typename T1::value_t >::type_t >::type norm2(const T1 &x)
compute ∥x∥₂
Definition: Algebra.hh:408

Other functions, involving more matrices are:

Matrix< double > M1( 10, 20 );
Matrix< double > M2( 10, 20 );
Matrix< double > M3( 10, 10 );
Matrix< double > M4( 10, 10 );
prod( 1.0, M1, transposed( M2 ), 0.0, M3 );
add( 2.0, M3, M4 );
std::enable_if_t< is_matrix< T2 >::value &&is_matrix< T3 >::value &&is_matrix< T4 >::value &&is_same_type< T1, typename T2::value_t >::value &&is_same_type< T1, typename T3::value_t >::value &&is_same_type< T1, typename T4::value_t >::value, void > prod(const T1 alpha, const T2 &A, const T3 &B, const T1 beta, T4 &C)
compute C ≔ β·C + α·A·B
Definition: Algebra.hh:1031

Inversion and factorisation functions are:

invert( M ); // compute M^-1
lu( M ); // compute LU factorisation
ldlt( M ); // compute LDL^T factorisation
ldlh( M ); // compute LDL^H factorisation
llt( M ); // compute LL^T factorisation
llh( M ); // compute LL^H factorisation
qr( M, R ); // compute QR factorisation
void invert(TMatrix *A, const TTruncAcc &acc, const inv_options_t &opts=inv_options_t())
compute
void lu(TMatrix *A, const TTruncAcc &acc, const fac_options_t &options=fac_options_t())
compute LU factorisation
std::enable_if_t< is_matrix< T1 >::value, void > llh(T1 &A)
compute L·L^H factorisation of given hermitian n×n matrix A
std::enable_if_t< is_matrix< T1 >::value, void > ldlh(T1 &A)
compute L·D·L^H factorisation of given hermitian n×n matrix A
std::enable_if_t< is_matrix< T1 >::value, void > ldlt(T1 &A)
compute L·D·L^T factorisation of given symmetric n×n matrix A
std::enable_if_t< is_matrix< T1 >::value, void > llt(T1 &A)
compute L·L^T factorisation of given symmetric n×n matrix A
std::enable_if_t< is_matrix< T1 >::value, void > qr(T1 &A, Matrix< typename T1::value_t > &R)
Compute QR factorisation of the matrix A.
Remarks
All factorisation functions are without pivoting!

It is also possible to solve matrix equations of the form

\[ A X = B \]

with known \(A\) and \(b\):

Matrix< double > A( 10, 10 );
Matrix< double > X( 10, 10 );
solve( A, X );
std::enable_if_t< is_matrix< T1 >::value &&is_vector< T2 >::value &&is_same_type< typename T1::value_t, typename T2::value_t >::value, void > solve(T1 &A, T2 &b)
solve A·x = b with known A and b; x overwrites b (A is overwritten upon exit!)
Definition: Algebra.hh:990

X will initially contain the matrix \(B\). Furthermore, the content of the matrix A will be destroyed during computation.

If \(A\) is a triangular matrix, the equation may take the more general forms:

\[ A X = \alpha B \quad \textrm{or} \quad X A = \alpha B \]

The corresponding side from which to apply \(A\) to \(X\) is defined by a parameter (form_left or form_right):

Matrix< double > A( 10, 10 );
Matrix< double > X( 10, 5 );
solve_tri( from_left, upper_triangular, unit_diag, 1.0, A, X );

Eigenvalues and singular values and vectors of a matrix can also be computed:

Vector< double > eig_val;
Matrix< double > eig_vec;
eigen( M, eigval, eigvec ); // compute all eigen values/vectors of M
Vector< double > S;
Matrix< double > V;
svd( M, S, V ); // compute SVD M = USV^H (U overwrites M)
sv( M, S ); // compute singular values of M
std::enable_if_t< is_matrix< T1 >::value, void > eigen(T1 &M, Vector< typename T1::value_t > &eig_val, Matrix< typename T1::value_t > &eig_vec)
compute eigenvalues and eigenvectors of matrix M
std::enable_if_t< is_matrix< T1 >::value, void > svd(T1 &A, Vector< typename real_type< typename T1::value_t >::type_t > &S, Matrix< typename T1::value_t > &V)
compute SVD decomposition of the nxm matrix A with n×min(n,m) matrix U, min(n,m)×min(n,...
std::enable_if_t< is_matrix< T1 >::value, void > sv(T1 &A, Vector< typename real_type< typename T1::value_t >::type_t > &S)
compute SVD decomposition of the nxm matrix A but return only the singular values S ∈ ℝ^min(n,...

If M is symmetric, selected eigen values may be computed:

eigen( M, B::Range( 0, 1 ), eigval, eigvec ); // compute first 2 eigen values/vectors