openblas系列之strsm

用户10637292

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

1060

前面两篇我们详细探讨了矩阵乘法，这一节我们将介绍strsm，用于求解三角矩阵方程。本节内容的验证环节，需要用到第一篇的sgemm（矩阵乘法），如有不了解的可自行参阅。

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

1.函数说明

1.1 函数公式：

A * X = alpha * B 或者 X * A = alpha * B

A：为上三角矩阵或者下三角矩阵。可以在X的左侧或右侧。

X：求解的矩阵。

B：输入矩阵，当计算完成后，会原地替换为输出矩阵。也就是X矩阵的值。

整个函数形式比较简单。接下来，列出fortran和cblas的函数定义

1.2 fortran形式：

//函数定义在：/opt/homebrew/opt/openblas/include/f77blas.h
//为了方便展示，将其中一些数据类型的别名替换为真实的数据类型
void strsm_(
        const char* side, // 字符串. A在X的左边或者右边. 左边:L 对应的公式 A * X = alpha * B  右边: R 对应的公式 X * A = alpha * B
        const char* uplo, // 字符串. A矩阵是上三角矩阵或者下三角矩阵. 上三角: U  下三角: L
        const char* transa, // 字符串. 矩阵A是否转置. 转置: T  不转置: N
        const char* diag, // 字符串. 是否为单位三角矩阵(对角线全部为1). 是: U  否: N
        int* m, // 矩阵B的行数
        int* n, // 矩阵B的列数
        float* alpha, // 矩阵B的系数
        float* a, // 矩阵A的数据
        int* lda,// 矩阵A的主维度. 列优先存储,即为矩阵A的行数
        float* b, // 矩阵B的数据
        int* ldb); // 矩阵B的主维度. 列优先存储,即为矩阵B的行数

1.3 cblas形式：

//函数定义在：/opt/homebrew/opt/openblas/include/cblas.h
//为了方便展示，将其中一些数据类型的别名替换为真实的数据类型
void cblas_strsm(
  const enum CBLAS_ORDER Order, // 枚举。矩阵使用行优先或者列优先存储。 行优先：CblasRowMajor  列优先：CblasColMajor
  const enum CBLAS_SIDE Side, // 枚举。矩阵A在X的左边或者右边。 左边：CblasLeft  右边：CblasRight
  const enum CBLAS_UPLO Uplo, // 枚举。矩阵A是上三角矩阵或下三角矩阵。 上三角矩阵：CblasUpper 下三角矩阵：CblasLower
  const enum CBLAS_TRANSPOSE TransA,// 枚举。矩阵A是否转置。 转置：CblasTrans  不转置：CblasNoTrans
  const enum CBLAS_DIAG Diag, // 枚举。矩阵A是否为单位三角矩阵。 单位三角矩阵：CblasUnit 非单位三角矩阵：CblasNonUnit
  const int M, // 矩阵B的行数
  const int N, // 矩阵B的列数
  const float alpha, // 系数alpha
  const float *A, // 矩阵A的数据。存储格式需要按照第一个参数order指定的行优先或者列优先
  const int lda, // 矩阵A的leading dimension. 如果是行优先，就是矩阵A的列数，如果是列优先，就是矩阵A的行数
  float *B, // 矩阵B的数据。存储格式也需要和order指定的保持一致
  const int ldb);// 矩阵B的leading dimension. 如果是行优先，就是矩阵B的列数，如果是列优先，就是矩阵B的行数

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

2.demo示例

2.1 示例1（列优先+非三角单位矩阵+A在左侧）

#include <string>
#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<<"\n";
    for (size_t i = 0; i < rows; i++){
        for (size_t j = 0; j < cols; j++){
            cout<<matrix[j*rows+i]<<" ";
        }
        cout<<endl;
    }
}
/**
    公式: A * X = alpha * B
*/
int main(int argc,char* argv[]){
    int M = 3;
    int N = 2;
    float alpha = 1.0;
    /** 列优先，注释为行展示
    A:
      11 28 33
      0 45  82
      0 0  77
    B:
      100 20
      56  63
      764 91
    */
    // 1. 构建数据
    vector<float> A = {11,0,0,28,45,0,33,82,77}; // 列优先存储
    vector<float> B = {100,56,764,20,63,91}; // 列优先存储
    // 2. 求解矩阵X
    cblas_strsm(CblasColMajor, // 列优先
        CblasLeft, // A矩阵在求解矩阵X的左边
        CblasUpper,  // 上三角矩阵
        CblasNoTrans, // 不转置
        CblasNonUnit, // 对角线的元素不为1，所以是非单位三角矩阵。
        M, // 结果矩阵X的行数
        N, // 结果矩阵X的列数
        alpha, // 系数alpha
        A.data(),M, // 矩阵A. 因为是列优先存储，所以它的leading dimension即为行数3
        B.data(), M);// 输入矩阵/输出矩阵B. 求解的X会覆盖在B矩阵上。因为是列优先存储，所以leading dimension为行数3
    // 3. 打印结果矩阵X
    print_col_matrix_by_row("列优先。求解的矩阵X", B.data(),M, N);
    cout<<"---------------------------------------"<<endl;
    // 4. 下面使用第一篇讲过的矩阵乘法A*X来重建矩阵B
    vector<float> C(6,0);
    float beta = 1.0;
    // 4.1 矩阵乘法A*X
    cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,
     M,N, M,alpha, 
     A.data(), M, 
     B.data(), M, beta, 
     C.data(), M);
    // 4.2 打印重建矩阵，观察是否与输入矩阵B相等
    print_col_matrix_by_row("列优先，重建矩阵B", C.data(), M, N);
    return 0;
}

// 编译链接
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include strsm_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o strsm_main
// 执行
./strsm_main
// 执行结果
列优先。求解的矩阵X:3x2
22.1794 0.190817 
-16.8358 -0.753535 
9.92208 1.18182 
---------------------------------------
列优先，重建矩阵B:3x2
100 20 
56 63 
764 91

2.2 示例2（列优先+三角单位矩阵+A在右侧）

#include <string>
#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<<"\n";
    for (size_t i = 0; i < rows; i++){
        for (size_t j = 0; j < cols; j++){
            cout<<matrix[j*rows+i]<<" ";
        }
        cout<<endl;
    }
}

/**
    公式 X * A = alpha * B
*/
int main(int argc,char* argv[]){
    int M = 3;
    int N = 2;
    float alpha = 1.0;
    /** 列优先，注释为行展示
    A:
        1 28
        0 1 
    B:
        100 20
        56  63
        764 91
    */
    // 1. 构建数据
    vector<float> A = {1,0,28,1}; // 列优先存储
    vector<float> B = {100,56,764,20,63,91}; // 列优先存储
    // 2. 求解矩阵X
    cblas_strsm(CblasColMajor, // 列优先
        CblasRight, // A矩阵在求解矩阵X的右边
        CblasUpper,  // A为上三角矩阵
        CblasNoTrans, // 不转置
        CblasUnit, // 对角线的元素为1，所以是单位三角矩阵。
        M, // 结果矩阵X的行数
        N, // 结果矩阵X的列数
        alpha, // 系数alpha
        A.data(),N, // 矩阵A. 因为是列优先存储，所以它的leading dimension即为行数2
        B.data(), M);// 输入矩阵/输出矩阵B. 求解的X会覆盖在B矩阵上。因为是列优先存储，所以leading dimension为行数3
    // 3. 打印结果矩阵X
    print_col_matrix_by_row("列优先。求解的矩阵X", B.data(),M, N);
    cout<<"---------------------------------------"<<endl;
    // 4. 下面使用第一篇讲过的矩阵乘法X*A来重建矩阵B，因为X存储在矩阵B中，所以下面的乘法实际为B*A
    vector<float> C(6,0);
    float beta = 1.0;
    // 4.1 矩阵乘法X*A
    cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, M,N, N,alpha, B.data(), M, A.data(), N, beta, C.data(), M);
    // 4.2 打印重建矩阵，观察是否与输入矩阵B相等
    print_col_matrix_by_row("列优先，重建矩阵B", C.data(), M, N);
    return 0;
}

// 编译链接
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include strsm_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o strsm_main
// 执行
./strsm_main
// 执行结果
列优先。求解的矩阵X:3x2
100 -2780 
56 -1505 
764 -21301 
---------------------------------------
列优先，重建矩阵B:3x2
100 20 
56 63 
764 91

需要特别注意的是，在strsm运算中，矩阵A位于方程左侧 (A * X = B) 还是右侧 (X * A = B)，会直接决定矩阵X与B的维度定义，这是正确调用该函数的关键。作为众多高阶矩阵分解（如LU、Cholesky）不可或缺的核心子程序，strsm在数值计算中扮演着基石般的角色，其重要性与高效性值得我们深入理解和掌握。

附录：

openblas系列之sgemm

openblas系列之ssyrk

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

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

存储