openblas系列之sgemm

用户10637292

发布于 2026-07-13 21:08:18

200

随着NLP领域的发展，AI真正的走入普通大众的视野。而线性代数运算作为人工智能的核心基础，其重要性不言而喻。无论是深度学习的tensorflow、pytorch，亦或是向量计算的faiss，都可以看到它的身影。openblas以其跨平台兼容，高性能已经成为事实上的标准实现。本系列探讨其内计算函数的使用，以期能在阅读源码，理解算法原理方面为读者提供些许助力。今天我们从最常用的函数sgemm单精度矩阵乘法开始。

1. 环境准备

以mac作为演示环境，demo程序尽量做到小而美、可复现、快速验证。

1.1 安装openblas

brew install openblas

安装完成后。

头文件目录：/opt/homebrew/opt/openblas/include

lib文件目录：/opt/homebrew/opt/openblas/lib

目录可能会有所不同，观察安装日志即可。

1.2 编译命令

g++ -std=c++17 -I/opt/homebrew/opt/openblas/include sgemm_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o sgemm_main

说明：

-I/opt/homebrew/opt/openblas/include 指定openblas头文件路径
sgemm_main.cpp 编写的demo示例
-L/opt/homebrew/opt/openblas/lib -lopenblas 指定openblas库
-o sgemm_main 指定输出的文件名

2. 深度解析

我们先来解释一下前置的知识。

2.1 矩阵存储：行优先和列优先

blas最初是由fortan语言编写的，而fortan语言默认是按照列优先的。openblas的实现既支持行优先也支持列优先。

行优先：

列优先：

2.2 函数说明：

函数支持的公式：

C = alpha * A * B + beta * C

最初的科学计算都是使用fortan语言编写的，随着C，C++等高级语言的发展，CBLAS在fortan核心之上包了一层C语言的壳，而openblas兼容了两种形式的接口。

fortan形式：

//函数定义在：/opt/homebrew/opt/openblas/include/f77blas.h
//为了方便展示，将其中的一些数据类型的别名替换为真实的类型
int sgemm_(const char* transa,
        const char* transb, 
        int* m, 
        int* n,
        int* k, 
        const float* alpha, 
        const float* a, 
        int* lda,
        const float* b, 
        int* ldb, 
        float* beta, 
        float* c, 
        int* ldc);

cblas形式：

//函数定义在：/opt/homebrew/opt/openblas/include/cblas.h
//为了方便展示，将其中的一些数据类型的别名替换为真实的类型
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);

其中，fortan语言形式少一个参数，只支持列优先的存储。而cblas接口支持行优先和列优先，可以根据第一个参数Order指定。其他的参数两种形式全部一致。

下面，我们详细讲解下cblas形式的函数：

// 公式：C = alpha * A * B + beta * C
void cblas_sgemm(
 const enum CBLAS_ORDER Order, // 枚举值。行优先：CblasRowMajor。列优先：CblasColMajor
 const enum CBLAS_TRANSPOSE TransA, // 枚举值，公式中的矩阵A是否转置。不转置：CblasNoTrans。转置：CblasTrans
 const enum CBLAS_TRANSPOSE TransB, // 枚举值，公式中的矩阵B是否转置。不转置：CblasNoTrans。转置：CblasTrans
 const int M, // 结果矩阵C的行数，也即矩阵A的行数
 const int N, // 结果矩阵C的列数，也即矩阵B的列数
 const int K, // 矩阵B的行数
 const float alpha, // A和B矩阵相乘的控制系数alpha。
 const float *A, // 矩阵A的原始数据。注意：存储格式需要跟第一个参数行优先或者列优先对齐。
 const int lda, // 矩阵A的leading dimension。如果是行优先，就是矩阵A一行数据的数量。如果是列优先，就是矩阵A一列数据的数量。
 const float *B, // 矩阵B的原始数据。注意：存储格式需要跟第一个参数行优先或者列优先对齐。
 const int ldb, // 矩阵B的leading dimension。如果是行优先，就是矩阵B一行数据的数量。如果是列优先，就是矩阵B一列数据的数量。
 const float beta, // 输入矩阵C的系数beta
 float *C, // 既是输入矩阵C，也是结果矩阵C。注意：存储格式需要跟第一个参数行优先或者列优先对齐。
 const int ldc);// 矩阵C的leading dimension。如果是行优先，就是矩阵C一行数据的数量。如果是列优先，就是矩阵C一列数据的数量。

假设矩阵A为3*2的矩阵，B为2*3的矩阵，那么输入矩阵C必须为3*3的矩阵，结果矩阵C也为3*3的矩阵。讲解参数有点抽象，我们从具体的例子开始出发吧。

2.3 demo示例（列优先）

#include <vector>
//引入c接口的blas头文件
#include <cblas.h>
#include <iostream>

using std::cout;
using std::endl;
using std::string;
using std::vector;

/**
 * 按照行展示的方式打印列优先存储的矩阵
 * name: 名字
 * matrix: 矩阵
 * rows：行数
 * cols: 列数
 */
void print_col_matrix_by_row(string name,float *matrix,int rows,int cols){
    cout<<name<<":"<<rows<<"x"<<cols<<endl;
    for (size_t i = 0; i < rows; i++){
        for (size_t j = 0; j < cols; j++){
            cout<<matrix[j*rows+i]<<" ";
        }
        cout<<endl;
    }
}

/**
 * 公式：
 * C = alpha * A * B + beta * C
 */
int main(int argc,char *argv[]){
    int M = 2;
    int N = 2;
    int K = 3;
    float alpha = 1.0;
    float beta = 1.0;
    // 行展示 不转置
    /*
    A:
        45 78 89
        23 12 56
    B:
        34 56
        28 67
        64 15
    C:
        1 2
        3 4
    */
    vector<float> A = {45,23,78,12,89,56};// 列优先存储
    vector<float> B = {34,28,64,56,67,15}; // 列优先存储
    vector<float> C = {1,3,2,4}; // 列优先存储
    cblas_sgemm(CblasColMajor,// 指定按照列优先存储
        CblasNoTrans, CblasNoTrans,// A,B矩阵均不转置
        M, // 结果矩阵C的行数
        N, // 结果矩阵C的列数
        K, // 矩阵B的行数
        alpha, // 矩阵相乘的系数
        A.data(), M, // 因为是按照列优先存储的，矩阵A的一列为2
        B.data(), K, // 因为是按照列优先存储的，矩阵B的一列为3
        beta, 
        C.data(), M);// 因为是按照列优先存储的，矩阵C的一列为2
    print_col_matrix_by_row("矩阵乘(列优先 不转置)",C.data(),M,N);
}

// 编译链接
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include sgemm_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o sgemm_main
// 执行
./sgemm_main

因为alpha为1.0，所以矩阵A*B的结果为：

// 45*34+78*28+89*64  45*56+78*67+89*15
// 23*34+12*28+56*64  23*56+12*67+56*15
9410 9081
4702 2932

因为beta为1.0，上述A*B的结果需要加上输入矩阵C，结果（会存储在矩阵C中）为：

9411 9083
4705 2936

当然beta也可以为0，或者输入矩阵C全部为0，那么结果矩阵C就是矩阵A*B的结果。

以上，就是关于sgemm函数的解析。重点关注行优先和列优先的存储方式以及leading dimension应该怎么填。本节使用列优先的方式进行演示，希望同学们能够在使用行优先或者加入转置的方式进行练习，其中的leading dimension会有很大的不同（对于理解函数至关重要）。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-01-04，如有侵权请联系 cloudcommunity@tencent.com 删除

存储

本文分享自松鼠Agent 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度