openblas第一弹:openblas 使用说明和常用接口介绍

openblas 使用说明

openblas 是一个开源的矩阵计算库，包含了诸多的精度和形式的矩阵计算算法。就精度而言，包括float和double，两种数据类型的数据，其矩阵调用函数也是不一样。不同矩阵，其计算方式也是有所不同，（姑且认为向量也是一维矩阵），例如，向量与向量之间的计算，向量与矩阵之间的计算，矩阵与矩阵之间的计算。

按照不同的计算需求openblas有如下几种

1. 第一层级：向量与向量

)

对于不同的精度，使用的函数也是不一样。例如cblas_daxpy表示double类型数据，cblas_daxpy表示float类型数据

• 头文件

  "cblas.h"
  

• 描述

  cblas_?axpy系列函数执行vector-vector的操作，即：
  y := alpha*x + y;
  其中alpha是标量，x，y是向量。
  

• 参数说明

  n       Specifies the number of elements in vectors x and y.
  a       Specifies the scalar a.
  x       Array, size at least (1 + (n-1)*abs(incx)).
  incx   Specifies the increment for the elements of x.   一般设置为1 
  y       Array, size at least (1 + (n-1)*abs(incy)).
  incy   Specifies the increment for the elements of y.  一般设置为1 
  

2. 第2层级：矩阵与向量

矩阵与向量的乘法种类也有很多，常规也就是上述三种，以cblas_sgemv为例：

• 头文件：

  "cblas.h"
  

• 描述：

  cblas_?ger系列函数执行matrix-vector的操作，即：
  y := alpha*A*x + beta*y,
  其中alpha，beta是标量，x，y是向量（x(m*1)，y(n*1)），A是一个m*n的矩阵。
  

• 详细接口：

  void cblas_sgemv (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE trans, const MKL_INT m, const MKL_INT n, const float alpha, const float *a, const MKL_INT lda, const float *x, const MKL_INT incx, const float beta, float *y, const MKL_INT incy)
  

• 详细接口的具体含义：

  • Layout：布局，也就是指定矩阵的排列方式，按照行优先(CblasRowMajor)还是列优先(CblasColMajor)，需要指出的是，在函数使用中，虽然是向量与矩阵的乘法，但是实际上，传递的是两个指针向量。这样对于多维矩阵而言，指明其排列方式，十分有必要。参数（CblasRowMajor，CblasColMajor）
  • trans： 转置特定操作：
    
    • if trans=CblasNoTrans, then y := alpha*A*x + beta * y;
    • if trans=CblasTrans, then y := alpha*A’x + beta y;
    • if trans=CblasConjTrans, then y := alpha *conjg(A’)*x + beta*y.
  • m：表示矩阵A的行m
  • n：表示矩阵A的列n
  • a：表示一个矩阵，size = lda *k
    
    • For Layout = CblasColMajor, k is n，
    • For Layout = CblasRowMajor, k is m.
  • alpha：一个标量，通常为1
  • x：矩阵向量
  • incx，incy，一般为1
  • beta: 一个标量参数，通常设置为 0，之前试过设置为1，直接导致了结果不正确；在矩阵或是向量元素的值较大的时候，beta这个值，为0还是1影响不大；

• 示例程序1：

void matrix_vector(){
    float array[6] = { 1,2,3,4,5,6 };
    float x[3] = { 1,2,3 };
    int m = 2; 
    int n = 3;
    float y[2] = {};
    cblas_sgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, array, n, x, 1, beta, y, inc_y);
    //cblas_sgemv(CblasColMajor, CblasTrans, n, m, alpha, array, n, x, 1, beta, y, inc_y);
    //两个结果一致
}

//打印结果
14,32

• 示例程序2：

void matrix_vector()
{
       int i = 0;
       double x[2] = { 1.0,2.0 };
       double y[3] = { 2.0,1.0,3.0};
       double A[6] = { 0 };
       blasint rows = 2, cols = 3;
       double alpha = 10;
       blasint inc_x = 1, inc_y = 1;
       blasint lda = 2;
       //矩阵按列优先存储
       //A <== alpha*x*y' + A （y'表示y的转置）
       cblas_dger(CblasColMajor, rows, cols, alpha, x, inc_x, y, inc_y, A, lda);
       for (i = 0; i<6; i++)
              std::cout << A[i] << " ";
       std::cout << "\n";
}
//打印结果
20 40 10 20 30 60

3. 第3层级：矩阵与矩阵

矩阵与矩阵的乘积，相较于前两种较为复杂，但是用途却比较大，上面能算的，它行；上面不能算的，它也行；

• 头文件：

  "cblas.h"
  

• 描述：

  C := alpha*op(A)*op(B) + beta*C,
  其中：op（A）表示A矩阵或是转置等
  A是一个M*K矩阵
  B是一个K*N矩阵
  C是一个M*N 矩阵

• 接口描述：

  void cblas_sgemm(OPENBLAS_CONST enum CBLAS_ORDER Order, OPENBLAS_CONST enum CBLAS_TRANSPOSE TransA, OPENBLAS_CONST enum CBLAS_TRANSPOSE TransB, OPENBLAS_CONST blasint M, OPENBLAS_CONST blasint N, OPENBLAS_CONST blasint K,
   OPENBLAS_CONST float alpha, OPENBLAS_CONST float *A, OPENBLAS_CONST blasint lda, OPENBLAS_CONST float *B, OPENBLAS_CONST blasint ldb, OPENBLAS_CONST float beta, float *C, OPENBLAS_CONST blasint ldc);
  

• 详细参数描述：

  • cblas_sgemm(order,transA,transB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC);

  • 详细参数与上述的矩阵与向量乘法类似

  • 以此函数为例：

    cblas_sgemm(order,transA,transB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC);
    
    alpha =1,beta =0 的情况下，等于两个矩阵相成。
    
    第一参数 oreder 候选值 有ClasRowMajow 和ClasColMajow 这两个参数决定一维数组怎样存储在内存中,
    
    一般用ClasRowMajow
    
    参数 transA和transB ：表示矩阵A，B是否进行转置。候选参数 CblasTrans 和CblasNoTrans.
    
    参数M：表示 A或C的行数。如果A转置，则表示转置后的行数
    
    参数N：表示 B或C的列数。如果B转置，则表示转置后的列数。
    
    参数K：表示 A的列数或B的行数（A的列数=B的行数）。如果A转置，则表示转置后的列数。
    
    参数LDA：表示A的列数，与转置与否无关。
    
    参数LDB：表示B的列数，与转置与否无关。
    
    参数LDC：始终=N
    
    NOTE: 首先判定是行优先还是列优先；
    再次,依照lda，ldb分别将两个矩阵按行、列分开
    再次,判定是否进行转置操作
    输出结果；
    

• 示例程序1：

#include <vector>
#include <iostream>
#include <cblas.h>
using namespace std;
int main() {
    const int M=4;
    const int N=2;
    const int K=3;
    const float alpha=1;
    const float beta=0;
    const int lda=M;
    const int ldb=K;
    const int ldc=N;
    const float A[K*M]={1,2,3,4,5,6,7,8,9,8,7,6};
    const float B[N*K]={5,4,3,2,1,0};
    float C[M*N];

    cblas_sgemm(CblasRowMajor, CblasTrans, CblasTrans, 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]<<" ";
       }
       cout<<endl;
    }
」
输出：
52 7
58 10
64 13
70 16

• 示例程序2：

#include <vector>
#include <iostream>
#include <cblas.h>
using namespace std;
int main() {
    const int M=4;
    const int N=2;
    const int K=3;
    const float alpha=1;
    const float beta=0;
    const int lda=K;
    const int ldb=N;
    const int ldc=N;
    const float A[M*K]={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(CblasRowMajor, CblasNoTrans, CblasNoTrans, 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]<<" ";
       }   
       cout<<endl;
    }  
} 

结果：
14 8 
41 26
68 44 
67 46 

参考链接：

1：本文主要参考链接，主要是inte的MKL库的部分接口说名，openblas的这种接口没找到
MKLAPI

2：本文第三部分主要参考此博客，详细说明了多维矩阵的乘法问题多维矩阵的乘法问题

https://software.intel.com/en-us/mkl-developer-reference-c-cblas-ger
。

3:官网：http://www.openblas.net/