Re: [Boost-Users] uBlas Lapack binds - part 2

From: Kresimir Fresl (fresl_at_[hidden])
Date: 2003-04-15 05:57:21

Hi Geert,

you wrote:

> I realise that the binding will need quite some work to get right 
> under VC6 (see previous posting). 

(This is a reply to both postings.)

Unfortunately (or, maybe, fortunately ;o), I don't have any
experience with VC6 (or VC7). Neither I have access to it,
so all that I will write will be just a guess -- I can't check it.

My first guess is that the problem is not in template specialization,
but in *partial* template specialization -- for example, in
ublas/traits.hpp there is a template class `type_traits<>' which
is *fully* specialized for float, double etc. and AFAIK VC6 doesn't
complain.

> I would welcome suggestions to do 
> the conversion directly (so without all the traits stuff) and not as 
> general. Goal would be to quickly start and then (hopefully) one day 
> be able to migrate to the traits version. Maybe the getrf function
> in Lapack would be a good starting point? Thanks ...

If I understand you correctly, by `to quickly start' and `to
migrate' you mean the following scenario:

-- write (and compile ;o) now something like:

    ublas::matrix<double, ublas::column_major> A (n, n);
    ublas::matrix<double, ublas::column_major> B (n, 1);
    std::vector<int> piv (n);
    // ...
    getrf (A, piv.begin(), piv.end());
    getrs ('N', A, piv.begin(), piv.end(), B);

-- and then, one day, with better compiler (which accepts
traits), compile it again, without to much modifications.

If you are interested only in lapack bindings (and not in blas),
then all you need is some substitute for matrix_traits<>
(that is, you don't need vector_traits<>, at least not
now, but see also Note [1] below).

You can write your own traits/traits.hpp, but without traits
classes ;o)

Namely, lapack bindings are written in terms of functions
matrix_size1(), matrix_size2(), matrix_storage() and
leading_dimension(), which are simple forwarding functions
that call corresponding static functions from matrix_traits<>
specializations -- so, in lapack/lapack.hpp there's no
direct mention of matrix_traits<>.

You can rewrite these functions to do `real work':
==========================================
template <typename MatrixT>
inline
int matrix_size1 (MatrixT const& m) { return m.size1(); }

template <typename MatrixT>
inline
int matrix_size2 (MatrixT const& m) { return m.size2(); }

template <typename MatrixT>
inline
int leading_dimension (MatrixT const& m) { return m.size1(); }

template <typename MatrixT>
inline
typename MatrixT::value_type* matrix_storage (MatrixT& m) {
     return m.data().begin();
}

template <typename MatrixT>
inline
typename MatrixT::value_type const*
matrix_storage (MatrixT const& m) {
     return m.data().begin();
}
==========================================
These should work with

     ublas::matrix<T, F, ublas::unbounded_array<T> >

but maybe not with with

     ublas::matrix<T, F, std::vector<T> >

Problem with this approach is that soon it will not work
anymore -- namely, getrf and getrs now use

     typedef typename matrix_type::value_type value_type;

but in generic implementation this should be

     typedef typename matrix_traits<matrix_type>::value_type value_type;

and, unfortunately, matrix_traits<> will reappear.
(Well, probably we can do something like

#ifdef BUGGY_COMPILER
     typedef typename matrix_type::value_type value_type;
#else
     typedef typename matrix_traits<matrix_type>::value_type value_type;
#endif

but ...)

Another possibility is to write your own traits/traits.hpp
and ublas_matrix.hpp with full specialization of matrix_traits<>:
=============================================
// in traits/traits.hpp
template <typename M>
struct matrix_traits {};

// in ublas_matrix.hpp
template<>
struct matrix_traits<
    boost::numeric::ublas::matrix<
       double,
       boost::numeric::ublas::column_major
    >
 >
{
public:
     typedef general_t matrix_structure;
     typedef double value_type;
     typedef double* pointer;
     typedef boost::numeric::ublas::matrix<
          double,
          boost::numeric::ublas::column_major
     > matrix_type;
     static pointer storage (matrix_type& m) {
          return m.data().begin();
     }
     static int size1 (matrix_type& m) { return m.size1(); }
     static int size2 (matrix_type& m) { return m.size2(); }
     // etc.
};

// same thing (copy&paste with a bit of post-editing) for
//    matrix<double> const
//    matrix<float>
//    matrix<float> const
//  etc.
=============================================

Note [1]:  as I said, currently you don't need vector_traits<>,
but with the next version of bindings you will probably need
them.  We recently discussed the interface of getrf, getrs etc.
Currently, getrf is:

   template < typename matrix_type, typename IpivIterator >
   int getrf(matrix_type& a, IpivIterator begin_ipiv, IpivIterator 
end_ipiv);

but in the next version it will be:

   template < typename matrix_type, typename int_vector_type >
   int getrf(matrix_type& a, int_vector_type& piv);

(as in the atlas/clapack bindings) and to get the starting address
of the pivot vector `piv' vector_traits<>::storage() will be needed
-- but then you can probably write full specialization of vector_traits
for std::vector<int>.

> By the way Toon, I think there is a typo in lapack.hpp:
[...]
>  int n = boost::numeric::bindings::traits::matrix_matrix_size2( a ) ;
>
> Shouldn't matrix_matrix_size2 be matrix_size2? 

Yes, it should.

Regards,

fres