# blaze / Vector and Matrix Customization

## Custom Data Members

So far the Blaze library does not provide a lot of flexibility to customize the data members of existing vector types and matrix types. However, to some extend it is possible to customize vectors and matrices by inheritance. The following example gives an impression on how to create a simple variation of `CompressedMatrix`, which automatically takes care of acquiring and releasing custom memory.

```template< typename Type                    // Data type of the matrix
, bool SO = defaultStorageOrder >  // Storage order
class MyCustomMatrix
: public CustomMatrix< Type, unaligned, unpadded, SO >
{
public:
explicit inline MyCustomMatrix( size_t m, size_t n )
, array_( new Type[m*n] )
{
this->reset( array_.get(), m, n );
}

private:
std::unique_ptr<Type[]> array_;
};
```

Please note that this is a simplified example with the intent to show the general approach. The number of constructors, the memory acquisition, and the kind of memory management can of course be adapted to specific requirements. Also, please note that since none of the Blaze vectors and matrices have virtual destructors polymorphic destruction cannot be used.

## Custom Operations

There are two approaches to extend Blaze with custom operations. First, the `map()` functions provide the possibility to execute componentwise custom operations on vectors and matrices. Second, it is possible to add customized free functions.

#### The map() Functions

Via the unary and binary `map()` functions it is possible to execute componentwise custom operations on vectors and matrices. The unary `map()` function can be used to apply a custom operation on each single element of a dense vector or matrix or each non-zero element of a sparse vector or matrix. For instance, the following example demonstrates a custom square root computation on a dense matrix:

```blaze::DynamicMatrix<double> A, B;

B = map( A, []( double d ) { return std::sqrt( d ); } );
```

The binary `map()` function can be used to apply an operation pairwise to the elements of two dense vectors or two dense matrices. The following example demonstrates the merging of two matrices of double precision values into a matrix of double precision complex numbers:

```blaze::DynamicMatrix<double> real{ { 2.1, -4.2 }, { 1.0,  0.6 } };
blaze::DynamicMatrix<double> imag{ { 0.3,  1.4 }, { 2.9, -3.4 } };

blaze::DynamicMatrix< complex<double> > cplx;

// Creating the matrix
//    ( (-2.1,  0.3) (-4.2, -1.4) )
//    ( ( 1.0,  2.9) ( 0.6, -3.4) )
cplx = map( real, imag, []( double r, double i ){ return complex( r, i ); } );
```

These examples demonstrate the most convenient way of defining a unary custom operation by passing a lambda to the `map()` function. Alternatively, it is possible to pass a custom functor:

```struct Sqrt
{
double operator()( double a ) const
{
return std::sqrt( a );
}
};

B = map( A, Sqrt() );
```

In order for the functor to work in a call to `map()` it must define a function call operator, which accepts arguments of the type of the according vector or matrix elements.

Although the operation is automatically parallelized depending on the size of the vector or matrix, no automatic vectorization is possible. In order to enable vectorization, a `load()` function can be added to the functor, which handles the vectorized computation. Depending on the data type this function is passed one of the following Blaze SIMD data types:

• SIMD data types for fundamental data types
• `blaze::SIMDint8`: Packed SIMD type for 8-bit signed integral data types
• `blaze::SIMDuint8`: Packed SIMD type for 8-bit unsigned integral data types
• `blaze::SIMDint16`: Packed SIMD type for 16-bit signed integral data types
• `blaze::SIMDuint16`: Packed SIMD type for 16-bit unsigned integral data types
• `blaze::SIMDint32`: Packed SIMD type for 32-bit signed integral data types
• `blaze::SIMDuint32`: Packed SIMD type for 32-bit unsigned integral data types
• `blaze::SIMDint64`: Packed SIMD type for 64-bit signed integral data types
• `blaze::SIMDuint64`: Packed SIMD type for 64-bit unsigned integral data types
• `blaze::SIMDfloat`: Packed SIMD type for single precision floating point data
• `blaze::SIMDdouble`: Packed SIMD type for double precision floating point data
• SIMD data types for complex data types
• `blaze::SIMDcint8`: Packed SIMD type for complex 8-bit signed integral data types
• `blaze::SIMDcuint8`: Packed SIMD type for complex 8-bit unsigned integral data types
• `blaze::SIMDcint16`: Packed SIMD type for complex 16-bit signed integral data types
• `blaze::SIMDcuint16`: Packed SIMD type for complex 16-bit unsigned integral data types
• `blaze::SIMDcint32`: Packed SIMD type for complex 32-bit signed integral data types
• `blaze::SIMDcuint32`: Packed SIMD type for complex 32-bit unsigned integral data types
• `blaze::SIMDcint64`: Packed SIMD type for complex 64-bit signed integral data types
• `blaze::SIMDcuint64`: Packed SIMD type for complex 64-bit unsigned integral data types
• `blaze::SIMDcfloat`: Packed SIMD type for complex single precision floating point data
• `blaze::SIMDcdouble`: Packed SIMD type for complex double precision floating point data

All SIMD types provide the `value` data member for a direct access to the underlying intrinsic data element. In the following example, this intrinsic element is passed to the AVX function `_mm256_sqrt_pd()`:

```struct Sqrt
{
double operator()( double a ) const
{
return std::sqrt( a );
}

SIMDdouble load( simd_double_t a ) const
{
return _mm256_sqrt_pd( a.value );
}
};
```

In this example, whenever vectorization is generally applicable, the `load()` function is called instead of the function call operator for as long as the number of remaining elements is larger-or-equal to the width of the packed SIMD type. In all other cases (which also includes peel-off and remainder loops) the scalar operation is used.

Please note that this example has two drawbacks: First, it will only compile in case the intrinsic `_mm256_sqrt_pd()` function is available (i.e. when AVX is active). Second, the availability of AVX is not taken into account. The first drawback can be alleviated by making the `load()` function a function template. The second drawback can be dealt with by adding a `simdEnabled()` function template to the functor:

```struct Sqrt
{
double operator()( double a ) const
{
return std::sqrt( a );
}

template< typename T >
T load( const T& a ) const
{
return _mm256_sqrt_pd( a.value );
}

template< typename T >
static constexpr bool simdEnabled() {
#if defined(__AVX__)
return true;
#else
return false;
#endif
}
};
```

The `simdEnabled()` function must be a `static`, `constexpr` function and must return whether or not vectorization is available for the given data type `T`. In case the function returns `true`, the `load()` function is used for a vectorized evaluation, in case the function returns `false`, `load()` is not called.

Note that this is a simplified example that is only working when used for dense vectors and matrices with double precision floating point elements. The following code shows the complete implementation of the according functor that is used within the Blaze library. The Blaze `Sqrt` functor is working for all data types that are providing a square root operation:

```namespace blaze {

struct Sqrt
{
template< typename T >
BLAZE_ALWAYS_INLINE auto operator()( const T& a ) const
{
return sqrt( a );
}

template< typename T >
static constexpr bool simdEnabled() { return HasSIMDSqrt<T>::value; }

template< typename T >
BLAZE_ALWAYS_INLINE auto load( const T& a ) const
{
BLAZE_CONSTRAINT_MUST_BE_SIMD_PACK( T );
return sqrt( a );
}
};

} // namespace blaze
```

The same approach can be taken for binary custom operations. The following code demonstrates the `Min` functor of the Blaze library, which is working for all data types that provide a `min()` operation:

```struct Min
{
explicit inline Min()
{}

template< typename T1, typename T2 >
BLAZE_ALWAYS_INLINE decltype(auto) operator()( const T1& a, const T2& b ) const
{
return min( a, b );
}

template< typename T1, typename T2 >
static constexpr bool simdEnabled() { return HasSIMDMin<T1,T2>::value; }

template< typename T1, typename T2 >
BLAZE_ALWAYS_INLINE decltype(auto) load( const T1& a, const T2& b ) const
{
BLAZE_CONSTRAINT_MUST_BE_SIMD_PACK( T1 );
BLAZE_CONSTRAINT_MUST_BE_SIMD_PACK( T2 );
return min( a, b );
}
};
```

For more information on the available Blaze SIMD data types and functions, please see the SIMD module in the complete Blaze documentation.

#### Free Functions

In order to extend Blaze with new functionality it is possible to add free functions. Free functions can be used either as wrappers around calls to the `map()` function or to implement general, non-componentwise operations. The following two examples will demonstrate both ideas.

The first example shows the `setToZero()` function, which resets a sparse matrix to zero without affecting the sparsity pattern. It is implemented as a convenience wrapper around the `map()` function:

```template< typename MT  // Type of the sparse matrix
, bool SO >    // Storage order
void setToZero( blaze::SparseMatrix<MT,SO>& mat )
{
(~mat) = blaze::map( ~mat, []( int ){ return 0; } );
}
```

The `blaze::SparseMatrix` class template is the base class for all kinds of sparse matrices and provides an abstraction from the actual type `MT` of the sparse matrix. However, due to the Curiously Recurring Template Pattern (CRTP) it also enables a conversion back to the actual type. This downcast is performed via the tilde operator (i.e. `operator~()`). The template parameter `SO` represents the storage order (`blaze::rowMajor` or `blaze::columnMajor`) of the matrix.

The second example shows the `countZeros()` function, which counts the number of values, which are exactly zero, in a dense, row-major matrix:

```template< typename MT >
size_t countZeros( blaze::DenseMatrix<MT,rowMajor>& mat )
{
const size_t M( (~mat).rows() );
const size_t N( (~mat).columns() );
size_t count( 0UL );

for( size_t i=0UL; i<M; ++i ) {
for( size_t j=0UL; j<N; ++j ) {
if( blaze::isDefault<strict>( (~mat)(i,j) ) )
++count;
}
}

return count;
}
```

The `blaze::DenseMatrix` class template is the base class for all kinds of dense matrices. Again, it is possible to perform the conversion to the actual type via the tilde operator.

The following two listings show the declarations of all vector and matrix base classes, which can be used for custom free functions:

```template< typename VT  // Concrete type of the dense or sparse vector
, bool TF >    // Transpose flag (blaze::columnVector or blaze::rowVector)
class Vector;

template< typename VT  // Concrete type of the dense vector
, bool TF >    // Transpose flag (blaze::columnVector or blaze::rowVector)
class DenseVector;

template< typename VT  // Concrete type of the sparse vector
, bool TF >    // Transpose flag (blaze::columnVector or blaze::rowVector)
class SparseVector;
```
```template< typename MT  // Concrete type of the dense or sparse matrix
, bool SO >    // Storage order (blaze::rowMajor or blaze::columnMajor)
class Matrix;

template< typename MT  // Concrete type of the dense matrix
, bool SO >    // Storage order (blaze::rowMajor or blaze::columnMajor)
class DenseMatrix;

template< typename MT  // Concrete type of the sparse matrix
, bool SO >    // Storage order (blaze::rowMajor or blaze::columnMajor)
class SparseMatrix;
```

## Custom Data Types

The Blaze library tries hard to make the use of custom data types as convenient, easy and intuitive as possible. However, unfortunately it is not possible to meet the requirements of all possible data types. Thus it might be necessary to provide Blaze with some additional information about the data type. The following sections give an overview of the necessary steps to enable the use of the hypothetical custom data type `custom::double_t` for vector and matrix operations. For example:

```blaze::DynamicVector<custom::double_t> a, b, c;
// ... Resizing and initialization
c = a + b;
```

The Blaze library assumes that the `custom::double_t` data type provides `operator+()` for additions, `operator-()` for subtractions, `operator*()` for multiplications and `operator/()` for divisions. If any of these functions is missing it is necessary to implement the operator to perform the according operation. For this example we assume that the custom data type provides the four following functions instead of operators:

```namespace custom {

double_t add ( const double_t& a, const double_t b );
double_t sub ( const double_t& a, const double_t b );
double_t mult( const double_t& a, const double_t b );
double_t div ( const double_t& a, const double_t b );

} // namespace custom
```

The following implementations will satisfy the requirements of the Blaze library:

```inline custom::double_t operator+( const custom::double_t& a, const custom::double_t& b )
{
}

inline custom::double_t operator-( const custom::double_t& a, const custom::double_t& b )
{
return sub( a, b );
}

inline custom::double_t operator*( const custom::double_t& a, const custom::double_t& b )
{
return mult( a, b );
}

inline custom::double_t operator/( const custom::double_t& a, const custom::double_t& b )
{
return div( a, b );
}
```

Blaze will use all the information provided with these functions (for instance the return type) to properly handle the operations. In the rare case that the return type cannot be automatically determined from the operator it might be additionally necessary to provide a specialization of the following four Blaze class templates:

```namespace blaze {

template<>
using Type = custom::double_t;
};

template<>
struct SubTrait<custom::double_t,custom::double_t> {
using Type = custom::double_t;
};

template<>
struct MultTrait<custom::double_t,custom::double_t> {
using Type = custom::double_t;
};

template<>
struct DivTrait<custom::double_t,custom::double_t> {
using Type = custom::double_t;
};

} // namespace blaze
```

The same steps are necessary if several custom data types need to be combined (as for instance `custom::double_t` and `custom::float_t`). Note that in this case both permutations need to be taken into account:

```custom::double_t operator+( const custom::double_t& a, const custom::float_t& b );
custom::double_t operator+( const custom::float_t& a, const custom::double_t& b );
// ...
```

Please note that only built-in data types apply for vectorization and thus custom data types cannot achieve maximum performance!

Previous: Configuration Files ---- Next: Error Reporting Customization

Updated