sample83.cpp

Go to the documentation of this file.

#include <iostream>
#include <iomanip>

extern "C" {
#include <atlas/cblas.h>
  //#include <cblas.h>
}

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/operation_blocked.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/timer.hpp>

namespace boost { namespace numeric { namespace ublas {

    template<class M, std::size_t IBS, std::size_t JBS, std::size_t KBS, class E1, class E2>
    BOOST_UBLAS_INLINE
    M
    block_prod_goto (const matrix_expression<E1> &e1,
                     const matrix_expression<E2> &e2,
                     row_major_tag) {
        typedef M matrix_type;
        const std::size_t i_block_size = IBS;
        const std::size_t j_block_size = JBS;
        const std::size_t k_block_size = KBS;
        typedef const E1 expression1_type;
        typedef const E2 expression2_type;
        typedef typename M::size_type size_type;
        typedef typename M::value_type value_type;

        M m (e1 ().size1 (), e2 ().size2 ());
#ifdef BOOST_UBLAS_TYPE_CHECK
        matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
        indexing_matrix_assign (scalar_assign<value_type, value_type> (), cm, prod (e1, e2), row_major_tag ());
        disable_type_check = true;
#endif
        size_type i_size = e1 ().size1 ();
        size_type j_size = e2 ().size2 ();
        size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
        m.assign (zero_matrix<value_type> (m.size1 (), m.size2 ()));
        for (size_type i_begin = 0; i_begin < i_size; i_begin += i_block_size) {
            size_type i_end = i_begin + std::min (i_size - i_begin, i_block_size);
            for (size_type k_begin = 0; k_begin < k_size; k_begin += k_block_size) {
                size_type k_end = k_begin + std::min (k_size - k_begin, k_block_size);
                // const matrix<value_type, row_major, bounded_array<value_type, i_block_size * k_block_size> > e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
                const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
                for (size_type j_begin = 0; j_begin < j_size; j_begin += j_block_size) {
                    size_type j_end = j_begin + std::min (j_size - j_begin, j_block_size);
                    matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end));
                    if (boost::is_same<typename E2::orientation_category, column_major>::value) {
                        const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end));
                        m_range.plus_assign (prod (e1_range, e2_range));
                    } else {
                        // const matrix<value_type, column_major, bounded_array<value_type, k_block_size * j_block_size> > e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
                        const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
                        m_range.plus_assign (prod (e1_range, e2_range));
                    }
                }
            }
        }
#ifdef BOOST_UBLAS_TYPE_CHECK
        disable_type_check = false;
        BOOST_UBLAS_CHECK (equals (m, cm), internal_logic ());
#endif
        return m;
    }

    template<class M, std::size_t IBS, std::size_t JBS, std::size_t KBS, class E1, class E2>
    BOOST_UBLAS_INLINE
    M
    block_prod_goto (const matrix_expression<E1> &e1,
                     const matrix_expression<E2> &e2,
                     column_major_tag) {
        typedef M matrix_type;
        const std::size_t i_block_size = IBS;
        const std::size_t j_block_size = JBS;
        const std::size_t k_block_size = KBS;
        typedef const E1 expression1_type;
        typedef const E2 expression2_type;
        typedef typename M::size_type size_type;
        typedef typename M::value_type value_type;

        M m (e1 ().size1 (), e2 ().size2 ());
#ifdef BOOST_UBLAS_TYPE_CHECK
        matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
        indexing_matrix_assign (scalar_assign<value_type, value_type> (), cm, prod (e1, e2), row_major_tag ());
        disable_type_check = true;
#endif
        size_type i_size = e1 ().size1 ();
        size_type j_size = e2 ().size2 ();
        size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
        m.assign (zero_matrix<value_type> (m.size1 (), m.size2 ()));
        for (size_type j_begin = 0; j_begin < j_size; j_begin += j_block_size) {
            size_type j_end = j_begin + std::min (j_size - j_begin, j_block_size);
            for (size_type k_begin = 0; k_begin < k_size; k_begin += k_block_size) {
                size_type k_end = k_begin + std::min (k_size - k_begin, k_block_size);
                // const matrix<value_type, column_major, bounded_array<value_type, k_block_size * j_block_size> > e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
                const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
                for (size_type i_begin = 0; i_begin < i_size; i_begin += i_block_size) {
                    size_type i_end = i_begin + std::min (i_size - i_begin, i_block_size);
                    matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end));
                    if (boost::is_same<typename E1::orientation_category, row_major>::value) {
                        const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end));
                        m_range.plus_assign (prod (e1_range, e2_range));
                    } else {
                        // const matrix<value_type, row_major, bounded_array<value_type, i_block_size * k_block_size> > e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
                        const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
                        m_range.plus_assign (prod (e1_range, e2_range));
                    }
                }
            }
        }
#ifdef BOOST_UBLAS_TYPE_CHECK
        disable_type_check = false;
        BOOST_UBLAS_CHECK (equals (m, cm), internal_logic ());
#endif
        return m;
    }

    // Dispatcher
    template<class M, std::size_t IBS, std::size_t JBS, std::size_t KBS, class E1, class E2>
    BOOST_UBLAS_INLINE
    M
    block_prod_goto (const matrix_expression<E1> &e1,
                     const matrix_expression<E2> &e2) {
        typedef typename M::orientation_category orientation_category;
        return block_prod_goto<M, IBS, JBS, KBS> (e1, e2, orientation_category ());
    }

}}}

namespace ublas = boost::numeric::ublas;
using namespace std;

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 
                                        m, n, k, 1.0, &A.data()[0], A.size2(), &B.data()[0], B.size2(), 1.0, &C.data()[0], C.size2());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size1(), &B.data()[0], B.size2(), 1.0, &C.data()[0], C.size2());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size2(), &B.data()[0], B.size1(), 1.0, &C.data()[0], C.size2());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size1(), &B.data()[0], B.size1(), 1.0, &C.data()[0], C.size2());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size2(), &B.data()[0], B.size2(), 1.0, &C.data()[0], C.size1());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size1(), &B.data()[0], B.size2(), 1.0, &C.data()[0], C.size1());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::row_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size2(), &B.data()[0], B.size1(), 1.0, &C.data()[0], C.size1());
        }

    BOOST_UBLAS_INLINE
    void
    cblas_prod (ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& C,
                            const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& A,
                const ublas::matrix<double, ublas::column_major, ublas::unbounded_array<double> >& B) {
                // C(m,n) <- aplha A(m,k) B(k,n) + beta C(m,n)
                // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
            int m = C.size1();
                int n = C.size2();
                int k = A.size2();
        BOOST_UBLAS_CHECK (C.size1() == A.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (A.size2() == B.size1(), ublas::bad_size ());
        BOOST_UBLAS_CHECK (B.size2() == C.size2(), ublas::bad_size ());
                cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                                        m, n, k, 1.0, &A.data()[0], A.size1(), &B.data()[0], B.size1(), 1.0, &C.data()[0], C.size1());
        }

void zeromatrix(ublas::matrix<double,ublas::row_major>& A) {
 for(size_t i = 0; i<A.size1(); i++)
  for(size_t j = 0; j<A.size2(); j++)
   A(i,j) = 0.0;
}
void zeromatrix(ublas::matrix<double,ublas::column_major>& A) {
 for(size_t i = 0; i<A.size1(); i++)
  for(size_t j = 0; j<A.size2(); j++)
   A(i,j) = 0.0;
}
void initmatrix(ublas::matrix<double,ublas::row_major>& A) {
 for(size_t i = 0; i<A.size1(); i++)
  for(size_t j = 0; j<A.size2(); j++)
   A(i,j) = 1.0*rand();
}
void initmatrix(ublas::matrix<double,ublas::column_major>& A) {
 for(size_t i = 0; i<A.size1(); i++)
  for(size_t j = 0; j<A.size2(); j++)
   A(i,j) = 1.0*rand();
}

#define NOCHECK
#ifdef NOCHECK
template <class M, class E>
int equals(const M& lhs, const E& rhs)
{
        return 2;
}
#endif

typedef ublas::matrix<double,ublas::row_major> RT;
typedef ublas::matrix<double,ublas::column_major> CT;

#define mul(a,b,d) \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  for (size_t a=0;a<size;a++) \
   for (size_t b=0;b<size;b++) \
    for (size_t d=0;d<size;d++) \
     X(r,c) += A(r,k)*B(k,c); \
  time = t.elapsed(); \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B))<< endl; \
}

#define mul2(xtype) \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  X.plus_assign(prod(A,B)); \
  time = t.elapsed(); \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B)) << endl; \
}

#define mul3(xtype) \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  ublas::axpy_prod(A,B,X,false); \
  time = t.elapsed(); \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B)) << endl; \
}

#define mul4(xtype) \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  ublas::opb_prod(A,B,X,false); \
  time = t.elapsed(); \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B)) << endl; \
}

#define mul5(xtype) \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  X.plus_assign(ublas::block_prod<xtype, 64>(A,B)); \
  time = t.elapsed(); \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B)) << endl; \
}

#define mul6(xtype) \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  X.plus_assign(ublas::block_prod_goto<xtype, 64, 64, 64>(A,B)); \
  time = t.elapsed(); \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B)) << endl; \
}

#define mul7 \
{ \
  zeromatrix(X); \
  initmatrix(A); \
  initmatrix(B); \
  t.restart(); \
  cblas_prod(X, A, B); \
  time = t.elapsed(); \
  /* cout << setw(8) << (2.0*(double(size1)*double(size2)*double(size3))/time)/1e+6 << flush; */ \
  cout << setw(8) << time << flush; \
  cerr << equals(X,prod(A,B)) << endl; \
}

#define mulset(xtype,atype,btype) \
{ \
 xtype X(size1,size3); \
 atype A(size1,size2); \
 btype B(size2,size3); \
 mul5(xtype); \
 mul6(xtype); \
 mul7 ; \
 cout << endl; \
}

/*
 mul2(xtype); \
 mul3(xtype); \
 mul4(xtype); \
 mul(r,k,c); \
 mul(r,c,k); \
 mul(k,r,c); \
 mul(k,c,r); \
 mul(c,r,k); \
 mul(c,k,r); \
*/

void test(const size_t size1, const size_t size2, const size_t size3) {
  boost::timer t;
        double time;
  cout << "size of metrices - " << size1 << "x" << size2 << "*" << size2 << "x" << size3 << endl << endl;
  // cout << "      prod    axpy    opb     block   atlas   rkc     rck     krc     kcr     crk     ckr" <<
  cout << "      block   goto    atlas" <<
endl;
  cout << "RRR";mulset(RT,RT,RT);
  cout << "RRC";mulset(RT,RT,CT);
  cout << "RCR";mulset(RT,CT,RT);
  cout << "RCC";mulset(RT,CT,CT);

  cout << "CRR";mulset(CT,RT,RT);
  cout << "CRC";mulset(CT,RT,CT);
  cout << "CCR";mulset(CT,CT,RT);
  cout << "CCC";mulset(CT,CT,CT);

}

int main(int argc, char *argv[]) {
  int size1 = 100;
  if (argc > 1)
    size1 = ::atoi(argv [1]);
  int size2 = 100;
  if (argc > 2)
    size2 = ::atoi(argv [2]);
  int size3 = 100;
  if (argc > 3)
    size3 = ::atoi(argv [3]);
  test(size1, size2, size3);
}


// g++ -o sample83 sample83.cpp -I /LOCAL/gwinkler/include/ -DNDEBUG -O2 -g -L /LOCAL/gwinkler/lib -l cblas -l atlas

---