使用cblus_sgemm接口进行CPU矩阵乘法运算

        上一篇博客讲过了GPU比对乘法运算，使用的是cublas库的cublasSegmm和cublasGemmEx接口，速度很快，但当时也提出了问题，就是GPU的内存不够大，就算是损失一定精度压缩成int8型数据，也只是会比FP32型多存储4倍数据。如果在很大型的千亿级别数据项目中，就需要很大很大的显存才可以放得下这么多数据。但GPU成本比较高，显然也存在不便之处。并且在效能上，GPU比对在数据量很大时，前期在向GPU中存入数据的过程会比较长，一旦程序关闭重启，可能又需要等待很长的时间进行拷贝。所以这时，就可以考虑另一种矩阵乘法算法——CPU矩阵乘法。

         CPU矩阵乘法接口与GPU矩阵乘法接口名字很相似，为cblas_sgemm，可以在http://www.netlib.org/blas/中下载到源码，它的定义如下：

enum CBLAS_ORDER 
{
    CblasRowMajor = 101,      //行优先
    CblasColMajor = 102       //列优先
};

enum CBLAS_TRANSPOSE 
{
    CblasNoTrans = 111,       //矩阵不转置
    CblasTrans = 112,         //矩阵转置
    CblasConjTrans = 113      //矩阵按列转置
};

void cblas_sgemm(
                const enum CBLAS_ORDER Order,        //存储的有限性，有行优先和列优先(c语言是行优先)
                const enum CBLAS_TRANSPOSE TransA,   //A矩阵转置性
                const enum CBLAS_TRANSPOSE TransB,   //B矩阵转置性
                const int M,                         //A, C 的行数 
                const int N,                         //B, C 的列数
                const int K,                         //A 的列数和 B 的行数
                const float alpha,                   //运算式的 α 值
                const float *A,                      //A矩阵
                const int lda,                       //A矩阵的列数
                const float *B,                      //B矩阵
                const int ldb,                       //B矩阵的行数
                const float beta,                    //运算式的 β 值
                float *C,                            //结果矩阵C
                const int ldc                        //C矩阵的行数
                );

        接口定义也和GPU的cubalsSegmm非常相似，它可以完成的运算表达式如下：

        C=αop(A)op(B)+βC

        其中，α和β为标量，A、B、C分别是m×k、k×n、m×n的矩阵，当α = 1, β = 0时，即为：

        C=op(A)op(B)

        也就是做矩阵乘法。α和β分别对应cublasGemmEx中的const void *alpha、const void *beta。一个例子如下所示：

#include "cblas.h"
#include <iostream>

using namespace std;
#define M 5
#define N 3
#define K 2

int main()
{
    float *m_A = (float *)malloc(M * N * sizeof(float));
    float *m_B = (float *)malloc(N * K * sizeof(float));
    float *m_C = (float *)malloc(M * K * sizeof(float));

    for (int i = 0; i < N * M; i++) 
    {
        m_A[i] = (float)(rand() % 10 + 1);
        m_B[i] = (float)(rand() % 10 + 1);
    }

    // 打印待测试的矩阵
    cout << "矩阵 A :" << endl;
    for (int i = 0; i < N * M; i++)
    {
        cout << (int)m_A[i] << " ";
        if ((i + 1) % N == 0) 
            cout << endl;
    }
    cout << endl;
    
    cout << "矩阵 B :" << endl;
    for (int i = 0; i < N * K; i++)
    {
        cout << (int)m_B[i] << " ";
        if ((i + 1) % K == 0) 
            cout << endl;
    }
    cout << endl;

    float alpha = 1.0;
    float beta = 0;
    //计算CPU矩阵乘法
    cblas_sgemm(
        CblasRowMajor,
        CblasTrans,
        CblasTrans,
        M,
        N,
        K,
        alpha,
        m_A,
        N,
        m_B,
        K,
        beta,
        m_C,
        N
    );

    cout << "计算结果的转置 ( (A*B)的转置 )：" << endl;
 
    for (int i = 0; i < M * K; i++)
    {
        cout << m_C[i] << " ";
        if ((i+1) % K == 0) 
            cout << endl;
    }

    free(m_C);
    free(m_B);
    free(m_A);
    
    return 0;
}

        运算结果为A矩阵乘以B矩阵的结果，与GPU矩阵乘法相同。但是它省去了从内存向显存拷贝和从显存向内存拷贝数据的过程，且内存空间很大可以放很多数据，单纯相乘的次数依两个矩阵大小而定，(M*N)x(N*M)的效率远大于(1*N)x(N*M)的效率。且它支持多线程的方式，也就是我们可以把大矩阵分割成n个小矩阵，每个小矩阵起一个线程，例如使用pthread_t创建线程的方式，每个线程做一个MxM的矩阵相乘，效率能够达到4~5亿/秒，但最后的集中排序速率可能会耗时更多。