00001 #include <iostream>
00002 #include <iomanip>
00003
00004 #include <boost/numeric/ublas/duff.hpp>
00005
00006 #include <boost/numeric/ublas/matrix.hpp>
00007 #include <boost/numeric/ublas/operation_blocked.hpp>
00008 #include <boost/numeric/ublas/io.hpp>
00009 #include <boost/timer.hpp>
00010
00011 namespace ublas = boost::numeric::ublas;
00012 using namespace std;
00013
00014 void initmatrix(ublas::matrix<double>& A) {
00015 for(size_t i = 0; i<A.size1(); i++)
00016 for(size_t j = 0; j<A.size2(); j++)
00017 A(i,j) = 1.0*rand();
00018 }
00019 void initmatrix(ublas::matrix<double,ublas::column_major>& A) {
00020 for(size_t i = 0; i<A.size1(); i++)
00021 for(size_t j = 0; j<A.size2(); j++)
00022 A(i,j) = 1.0*rand();
00023 }
00024
00025 template <class M>
00026 void zeromatrix(M& X)
00027 {
00028 X.assign( ublas::zero_matrix<double>(X.size1(), X.size2()) );
00029 }
00030
00031 typedef ublas::matrix<double,ublas::column_major> CT;
00032 typedef ublas::matrix<double,ublas::row_major> RT;
00033
00034 #define mul(a,b,m) \
00035 { \
00036 initmatrix(A); \
00037 initmatrix(B); \
00038 t.restart(); \
00039 for (size_t a=0;a<size;a++) \
00040 for (size_t b=0;b<size;b++) \
00041 size_t m = 0; \
00042 DD (size, 4, e, ((X(r,c) += A(r,k)*B(k,c)), ++m) ); \
00043 time = t.elapsed(); \
00044 cout << setw(8) << time << flush; \
00045 }
00046
00047 #define mul2(xtype) \
00048 { \
00049 zeromatrix(X); \
00050 initmatrix(A); \
00051 initmatrix(B); \
00052 t.restart(); \
00053 X.plus_assign(prod(A,B)); \
00054 time = t.elapsed(); \
00055 cout << setw(8) << time << flush; \
00056 cerr << equals(X,prod(A,B)) << endl; \
00057 }
00058
00059 #define mul3(xtype) \
00060 { \
00061 zeromatrix(X); \
00062 initmatrix(A); \
00063 initmatrix(B); \
00064 t.restart(); \
00065 X.plus_assign(ublas::block_prod<xtype>(A,B,64)); \
00066 time = t.elapsed(); \
00067 cout << setw(8) << time << flush; \
00068 cerr << equals(X,prod(A,B)) << endl; \
00069 }
00070
00071 #define mulset(xtype,atype,btype) \
00072 { \
00073 xtype X(size,size); \
00074 atype A(size,size); \
00075 btype B(size,size); \
00076 mul2(xtype); \
00077 mul3(xtype); \
00078 mul(r,k,c); \
00079 mul(r,c,k); \
00080 mul(k,r,c); \
00081 mul(k,c,r); \
00082 mul(c,r,k); \
00083 mul(c,k,r); \
00084 cout << endl; \
00085 }
00086
00087 void test2(const size_t size) {
00088 boost::timer t;
00089 double time;
00090 cout << "size of metrices - " << size << "x" << size << endl << endl;
00091 cout << " ublas block rkc rck krc kcr crk ckr" <<
00092 endl;
00093 cout << "RRR";mulset(RT,RT,RT);
00094 cout << "RRC";mulset(RT,RT,CT);
00095 cout << "RCR";mulset(RT,CT,RT);
00096 cout << "RCC";mulset(RT,CT,CT);
00097 cout << "CRR";mulset(CT,RT,RT);
00098 cout << "CRC";mulset(CT,RT,CT);
00099 cout << "CCR";mulset(CT,CT,RT);
00100 cout << "CCC";mulset(CT,CT,CT);
00101 }
00102
00103 int main(int argc, char* argv[]) {
00104 test2(100);
00105 }
