As an introduction; my goal is to link together the uBlas
matrix/vector classes with the Intel MKL library (I looked at ATLAS,
but the Lapack coverage was insufficient). I have tried the bindings
code in the Boost Sandbox (on which Toon and Kresimir are working),
but VC6 falls over.
I managed to work around the partial template specialisation problem
(see my previous postings), thanks to Kresimir's and Toon's
feedback. There is however one remaining question I need some help
with.
First the question; after that I post all the code which I've hacked
together to assist other interested users using VC6.
Due to my ignorance, I am obliged to pass the orientation as a
template parameter to the matrix_traits class. What I would like to
do is to obtain this information from the matrix class directly. The
issue is that matrix_type::functor_type does not exist for all
matrix classes. For example, hermitian, triangular and symmetrix
matrix use matrix_type::functor2_type to store the orientation of
the matrix. With a good compiler (i.e. *NOT* VC6), one could
specialise matrix_traits for these matrix classes (at least that is
what I think - but this isn't done yet in the ublas_matrix.hpp file
in the Sandbox), here however I'm stuck. As a hack, I allow the
orientation to be passed as a template parameter. Is there a better
way to do this (the matrix class M knows the orientation already)?
// M is type of matrix, RC is the row/column orientation
template <typename M, // type of matrix
typename RC = boost::numeric::ublas::row_major>
struct matrix_traits
{
public:
typedef RC row_column_major;
typedef M matrix_type;
...
static int stride1 (matrix_type& m) {
//return matrix_type::functor_type::one1 (m.size1(), m.size2
());
return row_column_major::one1 (m.size1(), m.size2());
}
...
};
Here is the code I wrote:
---------------------//----------------------
traits.h:
// (-*- C++ -*- header)
/*
*
* Copyright (c) Kresimir Fresl and Toon Knapen 2002
*
* Permission to copy, modify, use and distribute this software
* for any non-commercial or commercial purpose is granted provided
* that this license appear on all copies of the software source
code.
*
* Authors assume no responsibility whatsoever for its use and makes
* no guarantees about its quality, correctness or reliability.
*
* First author acknowledges the support of the Faculty of Civil
* Engineering, University of Zagreb, Croatia.
*
* --> Modified/simplified for VC6
*/
#ifndef BOOST_NUMERIC_BINDINGS_TRAITS_TRAITS_HPP
#define BOOST_NUMERIC_BINDINGS_TRAITS_TRAITS_HPP
#include <cstddef>
//namespace boost { namespace numeric { namespace bindings {
namespace traits {
// vector_traits<> base:
// .. implementation -- do not use it directly
template <typename V, typename T = typename V::value_type>
struct default_vector_traits {
typedef T value_type;
typedef T* pointer;
// assumption: iterator==pointer
// .. e.g. ublas::(un)bounded_array
static pointer storage (V& v) { return v.begin(); }
static const T* storage (const V& v) { return v.begin(); }
static int size (V& v) { return v.size(); }
static int size (const V& v) { return v.size(); }
static int stride (V& v) { return 1; }
static int stride (const V& v) { return 1; }
};
// vector_traits<> generic version:
template <typename V>
struct vector_traits : default_vector_traits<V> {};
// matrix structure tags:
struct general_t {};
struct symmetric_t {};
struct symmetric_packed_t {};
struct hermitian_t {};
struct hermitian_packed_t {};
// triangular, banded etc.
struct unknown_structure_t {};
/*
// matrix_traits<> generic version:
template <typename M>
struct matrix_traits {
// typedefs:
// matrix_structure
// value_type
// pointer
// ordering_type
// static functions:
// pointer storage()
// int storage_size()
// int size1()
// int size2()
// int leading_dimension() // not all matrix types
// symmetric/hermitian typedefs:
// uplo_type
};
*/
// storage ordering tags:
struct row_major_t {};
struct column_major_t {};
// upper/lower tags:
struct upper_t {};
struct lower_t {};
// free accessor functions:
template <typename VectorT>
inline
typename vector_traits<VectorT>::pointer vector_storage (VectorT&
v) {
return vector_traits<VectorT>::storage (v);
}
template <typename VectorT>
inline
int vector_size (VectorT& v) {
return vector_traits<VectorT>::size (v);
}
template <typename VectorT>
inline
int vector_stride (VectorT& v) {
return vector_traits<VectorT>::stride (v);
}
// All the matrix_traits<MatrixT> code has been removed here ...
}}}}
#endif // BOOST_NUMERIC_BINDINGS_TRAITS_TRAITS_HPP
---------------------//----------------------
ublas_matrix.h:
// (-*- C++ -*- header)
/*
*
* Copyright (c) Kresimir Fresl and Toon Knapen 2002
*
* Permission to copy, modify, use and distribute this software
* for any non-commercial or commercial purpose is granted provided
* that this license appear on all copies of the software source
code.
*
* Authors assume no responsibility whatsoever for its use and makes
* no guarantees about its quality, correctness or reliability.
*
* First author acknowledges the support of the Faculty of Civil
* Engineering, University of Zagreb, Croatia.
*
* --> Modified/simplified for VC6
*/
#ifndef BOOST_NUMERIC_BINDINGS_TRAITS_UBLAS_MATRIX_H
#define BOOST_NUMERIC_BINDINGS_TRAITS_UBLAS_MATRIX_H
#include <boost/numeric/ublas/matrix.hpp>
#include "traits.h"
namespace boost { namespace numeric { namespace bindings { namespace
traits {
// Need to find a way to extract the orientation class (row /
column)
template <typename M, //
type of matrix
typename RC = boost::numeric::ublas::row_major> // row
or column major
struct matrix_traits
{
public:
typedef RC row_column_major;
typedef general_t matrix_structure;
typedef M::value_type value_type;
typedef M::value_type *pointer;
typedef M matrix_type;
static pointer storage (matrix_type& m) {
typedef typename matrix_type::array_type array_type;
return vector_traits<array_type>::storage (m.data());
}
static const value_type *storage (const matrix_type& m) {
typedef typename matrix_type::array_type array_type;
return vector_traits<array_type>::storage(m.data());
}
static int size1 (matrix_type& m) { return m.size1(); }
static int size1 (const matrix_type& m) { return m.size1(); }
static int size2 (matrix_type& m) { return m.size2(); }
static int size2 (const matrix_type& m) { return m.size2(); }
static int storage_size (matrix_type& m) {
return size1 (m) * size2 (m);
}
static int storage_size (const matrix_type& m) {
return size1 (m) * size2 (m);
}
// stride1 == distance (m (i, j), m (i+1, j))
static int stride1 (matrix_type& m) {
//return matrix_type::functor_type::one1 (m.size1(), m.size2
());
return row_column_major::one1 (m.size1(), m.size2());
}
static int stride1 (const matrix_type& m) {
//return matrix_type::functor_type::one1 (m.size1(), m.size2
());
return row_column_major::one1 (m.size1(), m.size2());
}
// stride2 == distance (m (i, j), m (i, j+1))
static int stride2 (matrix_type& m) {
//return matrix_type::functor_type::one2 (m.size1(), m.size2
());
return row_column_major::one2 (m.size1(), m.size2());
}
static int stride2 (const matrix_type& m) {
//return matrix_type::functor_type::one2 (m.size1(), m.size2
());
return row_column_major::one2 (m.size1(), m.size2());
}
static int leading_dimension (matrix_type& m) {
//return matrix_type::functor_type::size2 (m.size1(), m.size2
());
return row_column_major::size2 (m.size1(), m.size2());
}
static int leading_dimension (const matrix_type& m) {
//return matrix_type::functor_type::size2 (m.size1(), m.size2
());
return row_column_major::size2 (m.size1(), m.size2());
}
};
}}}}
#endif
---------------------//----------------------
boost_mkl_lapack.cpp:
namespace ublas = boost::numeric::ublas;
#ifdef BOOST_MSVC
#pragma warning (disable: 4355)
#pragma warning (disable: 4503)
#pragma warning (disable: 4786)
#endif
#define BOOST_UBLAS_USE_EXCEPTIONS
#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <vector>
// Intel MKL Lapack
#include <mkl/include/mkl_lapack64.h>
#include "ublas_matrix.h"
// To make life easiertypedef ublas::matrix<double,ublas::row_major>
Matrix;
typedef std::vector<int> IntVector;
template<class M>
int getrf(M& a, IntVector& ipiv)
{
typedef matrix_traits<M> mtraits;
typedef vector_traits<IntVector> vtraits;
int m = a.size1();
int n = a.size2();
assert(ipiv.size() == (m < n ? m : n));
double *aStorage = mtraits::storage(a);
int lda = mtraits::leading_dimension(a);
int *ipivStorage = vtraits::storage(ipiv);
int info; // return status
dgetrf(&m, &n, aStorage, &lda, ipivStorage, &info);
return info;
}
template<class M>
int getrs(const M& a, const IntVector& ipiv, M& b, char trans
= 'N' )
{
typedef matrix_traits<M> mtraits;
typedef vector_traits<IntVector> vtraits;
int m = a.size1();
int n = a.size1();
assert(ipiv.size() == (m < n ? m : n));
// Note: we need to cast away the const for a and ipiv
double *aStorage = const_cast<double*>(mtraits::storage(a));
double *bStorage = mtraits::storage(b);
int lda = mtraits::leading_dimension(a);
int ldb = mtraits::leading_dimension(b);
int *ipivStorage = const_cast<int*>(vtraits::storage(ipiv));
// if b is column major, take size2 (i.e. columns)
int nrhs = mtraits::size1(b);
int info; // return status
dgetrs(&trans, &n, &nrhs, aStorage, &lda, ipivStorage, bStorage,
&ldb, &info);
return info;
}
template <typename M>
void init_symm (M& m, char uplo = 'f') {
size_t n
= matrix_traits<M>::size1 (m);
for (size_t i = 0; i < n; ++i)
{
m(i, i) = n;
for (size_t j = i + 1; j < n; ++j) {
if (uplo == 'u' || uplo == 'U')
m(i, j) = n - (j - i);
else if (uplo == 'l' || uplo == 'L')
m(j, i) = n - (j - i);
else
m(i, j) = m(j, i) = n - (j - i);
}
}
}
// printing:
template <typename M>
void print_m (const M& m, char const* ch = 0) {
if (ch)
std::cout << ch << ":\n";
size_t sz1
= matrix_traits<M>::size1 (m);
size_t sz2
= matrix_traits<M>::size2 (m);
for (std::size_t i = 0 ; i < sz1 ; ++i) {
for (std::size_t j = 0 ; j < sz2 ; ++j)
std::cout << m(i, j) << " ";
std::cout << std::endl;
}
// std::cout << std::endl;
}
// printing:
template <typename V>
void print_v (V const& v, char const* ch = 0) {
if (ch)
std::cout << ch << ": ";
size_t sz
= vector_traits<V>::size (v);
for (std::size_t i = 0; i < sz; ++i)
std::cout << v[i] << " ";
std::cout << std::endl;
}
using std::size_t;
void main(void)
{
size_t n = 3;
Matrix a(n, n), acopy(n, n);
IntVector ipiv(n);
init_symm(a);
acopy = a;
print_m (a, "A");
// Now factor
int status = getrf(a, ipiv);
print_v(ipiv, "ipiv");
print_m (a, "A");
// The solution (and the right hand side) are row oriented
size_t nrhs = 2;
Matrix b(nrhs, n), x(nrhs, n);
int i;
for(i = 0; i < n; i++)
{
x(0,i) = (double)i + 1;
x(1,i) = (double)i + 4.0;
}
b = ublas::trans(ublas::prod(acopy, ublas::trans(x)));
print_m (b, "b original");
print_m (x, "x original");
status = getrs (a, ipiv, b); // solve from factorization
print_m (b, "X");
std::cout << std::endl;
// Now try the same for a symmetrix matrix
ublas::symmetric_matrix<double, ublas::lower> ml(n, n);
init_symm(ml);
print_m(ml, "symmetric matrix");
status = getrf(ml, ipiv);
print_v(ipiv, "ipiv");
print_m (a, "A");
}
You will need your project settings to be such that you link with
the Intel lib files and the DLL should be in your path (see Intel
MKL documentation). Your link tab in project/settings should contain
(the lib directory is set using tools/options, same for the include):
mkl_c.lib mkl_p3.lib mkl_lapack.lib
You may want to do this in a different namespace than the boost
namespace (I actually commented out all the namespace stuff above,
but reintroduced it in this posting - if it does not compile, than
that's probably the reason). More binds to Lapack still need to be
done, but this is a start already.
Geert
