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 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 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 #include "traits.h" namespace boost { namespace numeric { namespace bindings { namespace traits { // Need to find a way to extract the orientation class (row / column) template // 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::storage (m.data()); } static const value_type *storage (const matrix_type& m) { typedef typename matrix_type::array_type array_type; return vector_traits::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 #include #include <vector> // Intel MKL Lapack #include #include "ublas_matrix.h" // To make life easiertypedef ublas::matrix 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(mtraits::storage(a)); double *bStorage = mtraits::storage(b); int lda = mtraits::leading_dimension(a); int ldb = mtraits::leading_dimension(b); int *ipivStorage = const_cast(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 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