Provide information on how to work with other lapack routines (and other libraries), and how to extend blaze functionality yourself

[Jonas Glesaaen]

Would it be possible to get some more information on how to correctly use the blaze vectors and matrices with other blas and lapack functions? I read through the source code for the wrappers you have already written, but it would be good to have some information on how one can extend this oneself.

Very specifically, I need to compute some eigenvalues of various matrices, and need to know how to call the lapack-functions correctly. Reading the source I see that you call member-functions of the result after calling the storage destructors, which I'd very much like to learn the reason for.

As a general enquiry it would also be good to know how to interact with vector and matrix routines from other libraries, but I assume that might be difficult due to padding. A word or two would be nice anyway, as at the moment the only safe way I see of doing that is by copying the contents of a blaze container into a temporary raw array, and call the function on this matrix insted (possibly having to copy the data back again).

Finally, is there any way to provide information on how to write functions that take blaze-expressions? Or is this simply too complicated and cumbersome? I was writing my own trace function yesterday, and though that there is definitely some performance to be gained by working with the matrix expressions themselves when computing traces of large expressions.

Cheers, Jonas

[Klaus Iglberger]

Hi Jonas!

Thanks for pointing out these three issues. Please allow me to answer them separately:

Calling BLAS/LAPACK routines

Assuming that you have a column-major dense matrix, all you need are the member functions rows(), columns(), spacing(), and data(). rows() and columns() give you the current number of rows and columns, respectively, spacing() gives you the number of elements between the start of two columns and data() gives you a pointer to the internal storage. Most BLAS and LAPACK functions expect these parameters. For instance, the following example shows a call to the LAPACK dgetrf() and dgetri() functions to invert a matrix:

constexpr size_t N( 100UL );

blaze::DynamicMatrix<double,blaze::columnMajor> A( N, N );
// ... Initialize the matrix

int m   ( numeric_cast<int>( A.rows()    ) );  // == N
int n   ( numeric_cast<int>( A.columns() ) );  // == N
int lda ( numeric_cast<int>( A.spacing() ) );  // >= N
int info( 0 );
int lwork( n*lda );

const std::unique_ptr<int[]> ipiv( new int[N] );
const std::unique_ptr<double[]> work( new double[N] );

getrf( m, n, A.data(), lda, ipiv.get(), &info );
getri( n, A.data(), lda, ipiv.get(), work.get(), lwork, &info );

Based on this example you should be able to call any BLAS or LAPACK function for a column-major matrix. The complexity increases in case you have a row-major matrix, since LAPACK historically only supports column-major matrices. In case you have a row-major matrix, you have to adapt the parameters accordingly and you might have to accept that the result is the transpose of the result that you would get if you would use a column-major matrix. Alternatively, you can also use the according C wrappers for LAPACK (clapack, LAPACKE, etc.), which simplifies the task. Please note that due to the lda parameter padding is not an issue.

Compatibility of Blaze to Other Libraries

Blaze currently does not provide any guarantees on compatibility to libraries other than BLAS and LAPACK. This is especially true for sparse matrix libraries. Currently, there are also no conversion functions to vector or matrix formats of other libraries (especially including other C++ math libraries).

In case another library expects all matrix elements to lie contiguously in memory, the correct way is to copy the elements between data structures. Alternatively, you can disable padding via the usePadding switch in the <blaze/config/Optimizations.h> header. Then all matrix elements are guaranteed to be contiguous. However, since this has a significant negative performance impact on Blaze this approach is not recommended.

For the time being we don't have plans to add guarantees on compatibility. This might change, however, as we add features to Blaze, but this will always happen in conjunction with a new feature.

Extending Blaze

So far there is no guideline or wiki page on how to write functions for Blaze expressions. At this point, we also don't plan to add this. However, the following C++14 implementation of a trace() function should give you an impression on how to write a general function for all types of matrices (including expressions). Please note that this function works with all possible input matrices, but is not (!) optimized for performance. An optimized implementation will be added to Blaze in the context of issue #61.

template< typename MT, bool SO >
auto trace( const Matrix<MT,SO>& mat )
{
   if( !isSquare( ~mat ) ) {
      BLAZE_THROW_INVALID_ARGUMENT( "Invalid input matrix for trace computation" );
   }

   ElementType_<MT> tmp{};

   for( size_t i=0UL; i<(~mat).rows(); ++i ) {
      tmp += (~mat)(i,i);
   }

   return tmp;
}

The Matrix template is the base class for all matrices (including expressions). The template parameters MT and SO represent the actual type of the matrix and its storage order (rowMajor or columnMajor). Via the ~ operator you get a type-safe downcast to the actual type of the matrix.

If you have further questions on any of these topics we suggest that you contact us directly via email. We will try to answer as quickly as possible.

Since we realize that some examples on how to call BLAS/LAPACK functions with Blaze data structures would be helpful, we will keep the ticket open until we have provided additional documentation for this.

Best regards,

Klaus!

[Jonas Glesaaen]

Thank you very much for the detailed answer! It gave me enough pointers so that I should be able to implement most of what I need in Blaze.

Kind regards,

Jonas

[Klaus Iglberger]

The documentation for the Blaze LAPACK wrappers has been extended with additional examples. The documentation is immediately available via cloning the Blaze repository and will be officially released with Blaze 3.1.