openblas系列之ssyrk

用户10637292

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

140

在上一篇中，我们深入探讨了通用矩阵乘法sgemm，今天我们将介绍它的一个重要的「特化」成员ssyrk（Single-precision Symmetric Rank-K update，单精对称矩阵秩-K更新）

关于安装及行优先，列优先的概念，参考该系列的第一篇文章。

1.函数说明

函数公式：

C = alpha * A * A^T + beta * C

2. 为什么称之为「特化」？

通过对比sgemm（C = alpha * A * B + beta * C）的公式发现，ssyrk实际上是将矩阵B替换成了A（或其转置），这种看似简单的替换，却带来了本质的变化：

对称性：由于AA^T的结果必然是一个对称方阵，所以输入矩阵C也必须是对称方阵。虽然这一约束限制了函数的通用性，但却换取了极高的效率

计算量减半：程序只需计算并更新矩阵C的上三角或下三角部分。

内存优化：由于对称性，通常只需存储矩阵的一半数据（上三角或下三角），理论上节省了近50%的存储空间。

接下来，列出fortan和cblas的函数定义

fortan形式：

//函数定义在：/opt/homebrew/opt/openblas/include/f77blas.h
//为了方便展示，将其中的一些数据类型的别名替换为真实的类型
void ssyrk_(
        const char* uplo, // 更新上三角或者下三角。上三角：U 下三角：L
        const char* trans, // 是否转置
        int* n,
        int* k,
        float* alpha,
        float* a,
        int* lda,
        float* beta,
        float* c,
        int* ldc);

cblas形式：

//函数定义在：/opt/homebrew/opt/openblas/include/cblas.h
//为了方便展示，将其中的一些数据类型的别名替换为真实的类型
void cblas_ssyrk(
  const enum CBLAS_ORDER Order, 
  const enum CBLAS_UPLO Uplo, 
  const enum CBLAS_TRANSPOSE Trans,
  const int N, 
  const int K, 
  const float alpha, 
  const float *A, 
  const int lda, 
  const float beta, 
  float *C, 
  const int ldc);

在内存布局方面，fortan原生接口默认采用列优先存储；而cblas接口可通过首个参数Order控制，支持开发者在行优先与列优先之间自由切换。其余参数在逻辑上与cblas保持一致。

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

void cblas_ssyrk(
  const enum CBLAS_ORDER Order, // 枚举值。行优先：CblasRowMajor。列优先：CblasColMajor
  const enum CBLAS_UPLO Uplo, // 枚举值。上三角：CblasUpper 下三角：CblasLower
  const enum CBLAS_TRANSPOSE Trans, // 枚举值，公式中的矩阵A是否转置。不转置：CblasNoTrans。转置：CblasTrans
  const int N, // 结果方阵的行数/列数
  const int K, // 公共维度K。
  const float alpha, // 系数alpha
  const float *A, // 矩阵A的数据
  const int lda, // 矩阵A的leading dimension。如果是列优先，即为矩阵A一列数据的行数。如果是行优先，即为矩阵A一行数据的列数
  const float beta, // 系数beta
  float *C, // 既是输入矩阵C，也是结果矩阵C。注意：存储格式需要跟第一个参数行优先或者列优先对齐。
  const int ldc);// 矩阵C的leading dimension。如果是列优先，即为矩阵C一列数据的行数。如果是行优先，即为矩阵C一行数据的列数

3.demo示例

示例1（列优先+不转置）

#include<iostream>
#include<cblas.h>
#include<vector>
using std::cout;
using std::endl;
using std::vector;
using std::string;

/**
 * 按照行展示的方式打印列优先存储的矩阵
 * 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 * A^T + beta * C
 */
int main(int argc,char* argv[]){
    int N = 2;
    int K = 3;
    float alpha = 1.0;
    float beta = 1.0;
    /* 列优先 不转置 注释中的矩阵为行展示
    A:
        45 78 89
        23 12 56
    C:
        0 0
        0 0
    */
    vector<float> A = {45,23,78,12,89,56};// 列优先存储
    vector<float> C(N*N,0);// 列优先存储
    cblas_ssyrk(CblasColMajor, // 列优先
        CblasUpper, // 更新矩阵C的上三角
        CblasNoTrans, // 不转置
        N, // 结果方阵的N，即为2
        K, // 公共维度K，因为A为2*3的矩阵，即为3
        alpha, // 系数alpha
        A.data(), // 矩阵A的数据
        N, // 矩阵A的leading dimension。因为是列优先，所以是2
        beta, // beta系数
        C.data(), // 待更新的输入矩阵C，和结果矩阵C。
        N);// 列优先存储，为2
    print_col_matrix_by_row("列优先 不转置:",C.data(),N,N);
}

// 编译链接。关于环境安装，可参考第一篇文章：https://mp.weixin.qq.com/s/1M5ixC7dnucVOHzXKL5ucA
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include ssyrk_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o ssyrk_main
// 运行
./ssyrk_main
// 运行结果
列优先 不转置::2x2
16030 6955 
0 3809

示例2（列优先+转置+输入矩阵C）

#include<iostream>
#include<cblas.h>
#include<vector>
using std::cout;
using std::endl;
using std::vector;
using std::string;

/**
 * 按照行展示的方式打印列优先存储的矩阵
 * 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 * A^T + beta * C
 */
int main(int argc,char* argv[]){
    int N = 2;
    int K = 3;
    float alpha = 1.0;
    float beta = 1.0;
    /* 列优先 转置 注释中的矩阵为行展示
    A:
        45 23
        78 12
        89 56
    C:
        1 2
        3 4
    */
    vector<float> A = {45,78,89,23,12,56};// 列优先存储
    vector<float> C = {1,3,2,4}; // 列优先存储
    cblas_ssyrk(CblasColMajor, // 列优先存储
        CblasUpper, // 只更新C矩阵的上三角
        CblasTrans, // A转置
        N, // 矩阵C的行数/列数
        K, // 公共维度K
        alpha, // alpha系数
        A.data(), // 矩阵A
        K, // 矩阵A的leading dimension。因为是列优先存储，即为矩阵A的45 78 89，总共三个
        beta, // beta系数
        C.data(), // 输入矩阵C，也是输出矩阵C。无论是输入还是输出，都为列优先存储
        N);// 矩阵C的leading dimension。为2
    print_col_matrix_by_row("列优先 转置:",C.data(),N,N);
}

// 编译链接。关于环境安装，可参考第一篇文章：https://mp.weixin.qq.com/s/1M5ixC7dnucVOHzXKL5ucA
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include ssyrk_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o ssyrk_main
// 运行
./ssyrk_main
// 运行结果
列优先 转置::2x2
16031 6957 
3 3813

// 示例2流程分析，系数alpha和beta都为1.0。不再列出结算过程
A^T :
  45 78 89
  23 12 56
(A^T)^T = A:
  45 23
  78 12
  89 56
所以公式中：A * A^T = A^T * (A^T)^T = A^T * A
的结果为2*2的矩阵：
  16030 6955 
  0 3809
然后再加上输入矩阵C：
  1 2
  3 4
因为更新上三角矩阵，所以结果为：
  16031 6957 
  3 3813

以上便是对ssyrk函数的深度解析与实战演示。推荐同学们自行探索下cblas_sgemmt函数，对比下区别与优劣。

附录：

openblas系列之sgemm

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

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

函数