c – 布拉斯似乎很慢

#include <stdio.h>
#include <iostream>
#include <time.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sf_bessel.h>
using namespace std;


double diffclock(clock_t clock1,clock_t clock2) { double diffticks=clock1-clock2; double diffms=(diffticks*1000)/CLOCKS_PER_SEC; return diffms; }




int
main (void)
{
  double* a=new double[1000*1000];

  double* b=new double[1000*1000];

  double* c=new double[1000*1000];

  for (int i=0;i<1000*1000;i++){
  a[i]=i;
  b[i]=i/5+i*i/100;}

  gsl_matrix_view A = gsl_matrix_view_array(a,1000,1000);
  gsl_matrix_view B = gsl_matrix_view_array(b,1000);
  gsl_matrix_view C = gsl_matrix_view_array(c,1000);

  /* Compute C = A B */
  cout<<"start"<<endl;
  clock_t begin=clock();

  gsl_blas_dgemm (CblasNoTrans,CblasNoTrans,1.0,&A.matrix,&B.matrix,0.0,&C.matrix);
  clock_t end=clock();
  cout<<double(diffclock(end,begin))<<endl;

  return 0;  
}

我正在编译使用
// g -o program mm.cpp -I / home / gsl / include -lm -L / home / gsl / lib -lgsl -lgslcblas

python代码是

import time
import numpy as np



n=1000
a=np.zeros((n,n))
b=np.zeros((n,n))
for i in range(0,n):
    for j in range(0,n):
        a[i,j]=i*n+j
        b[i,j]=(i*n+j)/5+(n*i+j)**2/5
print "start"
start=time.time()
c=np.dot(a,b)
end=time.time()
print end-start

谢谢你的帮助！