我正在比较我机器上的矩阵乘法,看起来c blas非常慢.乘以1000×1000矩阵大约需要4秒钟,而
python所需的时间大约为1.5秒.我认为链接可能有问题,但我真的不知道如何解决这些问题.这是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
谢谢你的帮助！
解决方法
BLAS中的子程序是事实上的标准,并且存在实现该接口的大量优化和供应商特定的库. numpy和gsl都可以链接到各种不同的BLAS(或者在某些情况下使用它们自己的实现),但是从这个角度看,numpy和gsl几乎都是包装器 – 你获得的性能基本上只依赖于他们联系在一起的BLAS.
使用GSL,可以相对容易地与替代BLAS链接.这里有一些说明：http://www.gnu.org/software/gsl/manual/html_node/Linking-with-an-alternative-BLAS-library.html
英特尔的MKL是一种通常非常快的BLAS(至少如果你没有AMD cpu),但是很难连接起来.他们甚至有一个Web应用程序来帮助你编写链接：http://software.intel.com/en-us/articles/intel-mkl-link-line-advisor.我在OpenBLAS(http://www.openblas.net/)上运气不错,在i7-3770K cpu上获得的性能在MKL的1％或2％之内. OpenBLAS也很容易编译;它比ATLAS更令人头疼.
一旦你获得OpenBLAS,无论是从源代码编译还是从你的软件包管理器下载,如果你在* nix上,你的修改后的编译线基本上就是
g++ -o program mm.cpp -I/home/gsl/include -lm -L/home/gsl/lib -lgsl -lcblas -lopenblas
