Re: [boost] ublas & clapack

From: Joerg Walter (jhr.walter_at_[hidden])
Date: 2002-04-10 14:07:23

Hi Toon,

you wrote:

> On Wednesday 10 April 2002 08:51, you wrote:
> > P.S.: To give an impression, how the test integration works, I'll
include a
> > sample implementation and a cutout of the modified CLAPACK's blaswrap.h.
>
> and I suppose gmm is implemented using ublas routines.

Yep, here it is:

 template<class M1, class T1, class T2, class M2, class M3>
 M1 &
 gmm (char transpose_2, char transpose_3, M1 &m1, const T1 &t1, const T2
&t2, const M2 &m2, const M3 &m3) {
     if (transpose_2 == 'n') {
         if (transpose_3 == 'n') {
             return m1 = t1 * m1 + t2 * prod (m2, m3);
         } else if (transpose_3 == 't') {
             return m1 = t1 * m1 + t2 * prod (m2, trans (m3));
         } else if (transpose_3 == 'c') {
             return m1 = t1 * m1 + t2 * prod (m2, herm (m3));
         }
     } else if (transpose_2 == 't') {
         if (transpose_3 == 'n') {
             return m1 = t1 * m1 + t2 * prod (trans (m2), m3);
         } else if (transpose_3 == 't') {
             return m1 = t1 * m1 + t2 * prod (trans (m2), trans (m3));
         } else if (transpose_3 == 'c') {
             return m1 = t1 * m1 + t2 * prod (trans (m2), herm (m3));
         }
     } else if (transpose_2 == 'c') {
         if (transpose_3 == 'n') {
             return m1 = t1 * m1 + t2 * prod (herm (m2), m3);
         } else if (transpose_3 == 't') {
             return m1 = t1 * m1 + t2 * prod (herm (m2), trans (m3));
         } else if (transpose_3 == 'c') {
             return m1 = t1 * m1 + t2 * prod (herm (m2), herm (m3));
         }
     }
 }

> So this also compares
> the speed of ublas with BLAS

Yes.

> (I'v always been thinking in the other
> direction, that you made binding for BLAS to be able to call BLAS using
ublas
> matrices)

But then I couldn't use the CLAPACK test suite to verify functionality ;-)

> BTW, the matrix_adaptor type here is, I suppose, a template-type ?

Yes, a simple wrapper:

    template<class T>
    class matrix_adaptor {
    public:
 typedef T value_type;
 typedef T *pointer;

     matrix_adaptor (std::size_t size1, std::size_t size2, pointer data,
std::size_t leading_dimension):
     aa_ (leading_dimension * size2, data),
     m_ (leading_dimension, size2, aa_),
     mr_ (m_, range (0, size1), range (0, size2)) {
     assert (leading_dimension >= size1);
 }

 const matrix_range<matrix<value_type, column_major,
array_adaptor<value_type> > > &operator () () const {
     return mr_;
 }
 matrix_range<matrix<value_type, column_major, array_adaptor<value_type> > >
&operator () () {
     return mr_;
 }

    private:
 array_adaptor<value_type> aa_;
        matrix<value_type, column_major, array_adaptor<value_type> > m_;
 matrix_range<matrix<value_type, column_major, array_adaptor<value_type> > >
mr_;
    };

Best regards

Joerg