openblas系列之sgesvd

A = U * Σ * V^T

// 函数定义： https://netlib.org/lapack/explore-html/d1/d7f/group__gesvd_ga23ab797a9c7feb13b3f6feb0b52673b8.html
// 为了方便展示，使用原始数据类型
int sgesvd_(
        const char* jobu, // 指定如何计算左奇异向量U  'A'返回所有M个列向量(U为M*M)  'S'返回前min(M,N)个列向量(U为M*min(M,N)) 'O'返回前min(M,N)个列向量,并覆盖在输入矩阵A上(A被U覆盖) 'N'：不计算左奇异向量
        const char* jobvt, //指定如何计算右奇异向量V^T 'A'：返回所有N个行向量（VT为N×N）  'S'：返回前min(M,N)个行向量（VT为min(M,N)×N） 'O'：返回前min(M,N)个行向量，并覆盖在输入矩阵A上（如果JOBU和JOBVT都为'O'，则A被U覆盖，且VT正常返回，不覆盖A） 'N'：不计算右奇异向量 
        int* m, //矩阵A的行数
        int* n, //矩阵A的列数
        float* a, //矩阵A
        int* lda,//列优先的A的维度数
        float* s,// 奇异值,按降序排列
        float* u,//左奇异矩阵U
        int* ldu,//U的leading-dimension
        float* vt,//右奇异矩阵V的转置
        int* ldvt,// 右奇异矩阵V^T的leading-dimension
        float* work,//工作数组
        int* lwork,//工作数组的长度
        int* info);//状态码

// 函数定义在:/opt/homebrew/opt/openblas/include/lapacke.h
// 为了方便展示，将别名类型替换为真实类型
int LAPACKE_sgesvd(
  int matrix_layout, // 行优先或者列优先。行优先：LAPACK_ROW_MAJOR  列优先：LAPACK_COL_MAJOR
  char jobu, // 指定如何计算左奇异向量U  'A'返回所有M个列向量(U为M*M)  'S'返回前min(M,N)个列向量(U为M*min(M,N)) 'O'返回前min(M,N)个列向量,并覆盖在输入矩阵A上(A被U覆盖) 'N'：不计算左奇异向量
  char jobvt, //指定如何计算右奇异向量V^T 'A'：返回所有N个行向量（VT为N×N）  'S'：返回前min(M,N)个行向量（VT为min(M,N)×N） 'O'：返回前min(M,N)个行向量，并覆盖在输入矩阵A上（如果JOBU和JOBVT都为'O'，则A被U覆盖，且VT正常返回，不覆盖A） 'N'：不计算右奇异向量 
  int m, // 矩阵A的行数
  int n, // 矩阵A的列数
  float* a, // 输入矩阵A的原始数据
  int lda, // 矩阵A的leading dimension。根据行优先或者列优先来定。如果是行优先，即为列数。如果是列优先，即为行数
  float* s, // 输出结果：奇异值数组，为Σ矩阵的对角线元素，按照降序排列。大小为min(M,N)
  float* u, // 输出结果：左奇异矩阵
  int ldu, // 左奇异矩阵的 leading dimension。根据行优先或者列优先来定。如果是行优先，即为列数。如果是列优先，即为行数
  float* vt, // 输出结果：右奇异矩阵的转置
  int ldvt,  // 右奇异矩阵转置的leading dimension。根据行优先或者列优先来定。如果是行优先，即为列数。如果是列优先，即为行数
  float* superb ); // 中间计算结果，一般不会使用。工作空间中前1到min(M,N)的计算结果。大小为：min(M,N)-1

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
// 注意：本节所用的奇异值分解函数只需包含lapacke.h头文件即可。cblas.h文件是为了重建验证使用的
#include <cblas.h>
#include <lapacke.h>
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] << "\t";
    }
    cout << endl;
  }
}
/**
    公式：A = U * Σ * V^T
*/
int main(int argc, char *agrv[]) {
  int M = 3;
  int N = 2;
  int min = std::min(M, N);  // 用于处理奇异值数组
  /**
  列优先 注释为行展示
  A:
      1 2
      3 4
      5 6
  */
  // 1. 输入矩阵A
  vector<float> A = {1, 3, 5, 2, 4, 6};  // 列优先存储
  // 2. 准备输出空间
  vector<float> S(min);  // 奇异值数组，即为Σ的对角线元素。大小为 min(M,N)
  vector<float> U(M * M);   // 左奇异向量，3*3的矩阵
  vector<float> VT(N * N);  // 右奇异向量，2*2的矩阵
  vector<float> superb(
      min - 1);  // 工作空间中前1到min(M,N)的计算结果。大小为：min(M,N)-1
  // 3. 奇异值分解
  LAPACKE_sgesvd(
      LAPACK_COL_MAJOR,  // 列优先
      'A',               // 计算全部左奇异值向量，大小为M*M
      'A',               // 计算全部右奇异值向量，大小为N*N
      M,                 // 矩阵A的行数
      N,                 // 矩阵A的列数
      A.data(),          // 矩阵A的原始数据
      M,  // 矩阵A的leading dimension.因为第一个参数指定了列优先，所以为行数3。
      S.data(),   // 奇异值数组
      U.data(),   // 左奇异向量
      M,          // 左奇异向量U的leading-dimension.
                  // 因为第一个参数指定了列优先，所以为行数M
      VT.data(),  // 右奇异向量
      N,          // 右奇异向量V转置的leading-dimension.
                  // 因为第一个参数指定了列优先，所以为行数2
      superb
          .data());  // 工作空间（中间结果）前1到min(M,N)的计算结果。大小为：min(M,N)-1
  // 4. 打印计算结果
  print_col_matrix_by_row("奇异值S:", S.data(), 1, N);
  print_col_matrix_by_row("左奇异矩阵U", U.data(), M, M);
  print_col_matrix_by_row("右奇异矩阵转置VT:", VT.data(), N, N);
  // 以上，奇异值分解就已经计算完成了。下面是重建验证的过程
  // 5.重建矩阵A
  // 5.1 将奇异值恢复为Σ矩阵
  vector<float> sigma(M * N, 0);  // M*N的矩阵
  for (int i = 0; i < min; i++) {
    sigma[i * M + i] = S[i];
  }
  print_col_matrix_by_row("sigma", sigma.data(), M, N);
  // 5.2 将U * Σ的结果存入C矩阵
  vector<float> C(M * N, 0);
  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, 1.0, U.data(),
              M, sigma.data(), M, 1.0, C.data(), M);
  // 5.3 再计算C * VT将结果存入D矩阵，然后打印D矩阵与A进行对比
  vector<float> D(M * N, 0);
  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, 1.0, C.data(),
              M, VT.data(), N, 1.0, D.data(), M);
  print_col_matrix_by_row("重建后的A矩阵", D.data(), M, N);
  return 0;
}

// 编译链接
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include sgesvd_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o sgesvd_main
// 执行
./sgesvd_main
// 执行结果
奇异值S::1x2
9.52552 0.514301
左奇异矩阵U:3x3
-0.229848 0.883461 0.408247
-0.524745 0.240782 -0.816497
-0.819642 -0.401896 0.408249
右奇异矩阵转置VT::2x2
-0.619629 -0.784895
-0.784895 0.619629
sigma:3x2
9.52552 0
0 0.514301
0 0
重建后的A矩阵:3x2
1 2
3 4
5 6