Hi Jesse,

Thank you for your help. I downloaded the binding and linked it with blas and I saw very significant improvement. I realized that last time my program did not link with blas. Anyway, here is a copy of my test:


>> g++ test1.cpp -o test1 -Wall -O3 -DNDEBUG -lblas

>> ./test1 1000000 10 20
ublas  default implementation in blas
Time : 0.44897
default implementation in atlas
Time : 0.343825



----------------------------------------------
#include <iostream>
//#include <valarray>

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

#include <boost/numeric/bindings/atlas/cblas3.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
namespace atlas = boost::numeric::bindings::atlas;

# include <sys/time.h>

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


double gettime() {
    struct timeval t;
    gettimeofday(&t, NULL);
    return (double)t.tv_sec+t.tv_usec*0.000001;
}

void initmatrix(ublas::matrix<double>& 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();
}


int main(int argc, char *argv[])
{
    long int D,r,n ;

    //std::clock_t start;
    //double diff;
    double t1,t2 ;


    //D = 1000 ;
    D = atoi(argv[1]) ;
    //r = 1000 ;
    r = atoi(argv[2]) ;
    //n = 1000 ;
    n = atoi(argv[3]) ;

   {
   cout << "ublas  ";

   ublas::matrix<double> A(D,r);
   ublas::matrix<double,ublas::column_major> B(r,n);
   ublas::matrix<double> X(D,n);
   initmatrix(A);
   initmatrix(B);

   cout << "default implementation in blas" << endl ;
   t1 = gettime() ;
   //X = ublas::prod< ublas::matrix<double> >(A,B);
   noalias(X) = ublas::prod(A,B) ;
   t2 = gettime() ;
   cout << "Time : " <<  t2 - t1 << endl;

  
   }

   {
        ublas::matrix<double, ublas::column_major> A(D,r),B(r,n),C(D,n);
     initmatrix(A);
        initmatrix(B);

    cout << "default implementation in atlas" << endl ;
    t1 = gettime() ;
       atlas::gemm(A,B,C);
    t2 = gettime() ;
    cout << "Time : " <<  t2 - t1 << endl;

   
       
   } 


}