MKL库矩阵乘法

MKL库矩阵乘法

Posted on  2008年09月3日 by cici1020

我再次来学术一把。
我今天因为不想看3000多页的大文档，就在GOOGLE狂翻用Intel的Math Kernel Library的库如何作矩阵运算。最后还是看那个3000多页的大文档找到的解决方法。
前人栽树后人乘凉……
关键词：Intel MKL MKL库 矩阵乘积 矩阵乘法 样例 例子 代码 标程 sample code
这样以后人家都能搜到这里了:D

MKL的大文档里，有各种matrix*matrix的说明，不过都是fortran的。
而MKL自己也在sample文件夹下提供了sample,我是照着SAMPLE看，不太清楚CBLAS_ORDER的定义才查到下文的。
感觉下文写的很清楚:)

底下ZZ自http://blog.csdn.net/cleverysm/archive/2007/12/09/1925549.aspx

CBLAS的安装与使用
 
烤鱼片(@ eii.dlmu)
cleverysm@163.com
 
 
CBLAS是BLAS的C语言接口。BLAS的全称是Basic Linear Algebra Subprograms，中文大概可以叫做基础线性代数子程序。主要是用于向量和矩阵计算的高性能数学库。本身BLAS是用Fortran写的，为了方便C/C++程序的使用，就有了BLAS的C接口库CBLAS。BLAS的主页是 http://www.netlib.org/blas/，CBLAS的下载地址也可以在这个页面上找到。
 
CBLAS安装需要先装BLAS，从主页上下载 blas.tgz，解压，根据系统修改make.inc和Makefile，make，就会生成一个blas_LINUX.a文件。然后，下载cblas.tgz，解压，在目录下将Makefile.*文件改名或者做一个链接文件为Makefile.in文件，比如在linux下就是ln -s Makefile.LINUX  Makefile.in，根据具体情况修改Makefile.in文件，主要是BLAS的库文件路径BLLIB和CBLAS的安装目录CBDIR，make help就可以打印出可以使用的make命令，要生成全部文件就是用make all。在$(CBDIR)目录下的$(CBLIBDIR)将生成CBLAS的库文件$(CBLIB)，cblas_LINUX.a。
 
在CBLAS的安装目录$(CBDIR)下的src目录中有个cblas.h是包括的CBLAS的函数和常量的头文件，使用CBLAS的时候就需要这个头文件，同时还需要BLAS的库文件$(BLLIB )和CBLAS的库文件$(CBLIB)。
 
CBLAS/BLAS分为3个level，level1是用于向量的计算，level2是用于向量和矩阵之间的计算，level3是矩阵之间的计算。比如计算矩阵的乘法就是属于level3，这里就用矩阵乘法来学习使用CBLAS。
 
计算矩阵乘法的函数之一是 cblas_sgemm，使用单精度实数，另外还有对应双精度实数，单精度复数和双精度复数的函数。在此以 cblas_sgemm为例。
 
函数定义为：
void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
 
                 const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
 
                 const int K, const float alpha, const float *A,
 
                 const int lda, const float *B, const int ldb,
 
                 const float beta, float *C, const int ldc)
关于此函数的详细定义可以在 http://www.netlib.org/blas/sgemm.f找到，只不过是fortran语言的,这个C语言版的略有差别。
此函数计算的是 C = alpha*op( A )*op( B ) + beta*C,
 
const enum CBLAS_ORDER Order，这是指的数据的存储形式，在CBLAS的函数中无论一维还是二维数据都是用一维数组存储，这就要涉及是行主序还是列主序，在C语言中数组是用行主序，fortran中是列主序。我还是习惯于是用行主序，所以这个参数是用CblasRowMajor，如果是列主序的话就是CblasColMajor。
 
const enum CBLAS_TRANSPOSE TransA和 const enum CBLAS_TRANSPOSE TransB，这两个参数影响的是op( A )和op( B)，其中TransA = CblasNoTrans, op( A ) = A，TransA = CblasTrans, op( A ) = A’，TransA = CblasConjTrans, op( A ) = A’。 TransB类似。
 
const int M，矩阵A的行，矩阵C的行
const int N，矩阵B的列，矩阵C的列
const int K，矩阵A的列，矩阵B的行
const float alpha， const float beta，计算公式中的两个参数值，如果只是计算C=A*B，则alpha=1,beta=0
 
const float *A， const float *B， const float *C，矩阵ABC的数据
 
const int lda， const int ldb， const int ldc，在BLAS的文档里，这三个参数分别为ABC的行数，但是实际使用发现，在CBLAS里应该是列数。
 
我在这里计算两个简单矩阵的乘法。
A:
1,2,3
4,5,6
7,8,9
8,7,6
 
B:
5,4
3,2
1,0
 
程序代码：
//因为程序是C++，而CBLAS是C语言写的，所以在此处用extern关键字
extern"C"
{
    #include   
} 
#include
using namespace std;
int main(void) {
    
    const enum CBLAS_ORDER Order=CblasRowMajor;
    const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
    const enum CBLAS_TRANSPOSE TransB=CblasNoTrans;
    const int M=4;//A的行数，C的行数
    const int N=2;//B的列数，C的列数
    const int K=3;//A的列数，B的行数
    const float alpha=1;
    const float beta=0;
    const int lda=K;//A的列
    const int ldb=N;//B的列
    const int ldc=N;//C的列
    const float A[K*M]={1,2,3,4,5,6,7,8,9,8,7,6};
    const float B[K*N]={5,4,3,2,1,0};
    float C[M*N]; 
    cblas_sgemm(Order, TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
    for(int i=0;i<M;i++)
    {
       for(int j=0;j<N;j++)
       {
           cout<<C[i*N+j]<<"t";
       }
       cout<<endl;
    }
   return EXIT_SUCCESS;

}

from:

https://cici1020.wordpress.com/2008/09/03/mkl%E5%BA%93%E7%9F%A9%E9%98%B5%E4%B9%98%E6%B3%95/