openblas系列之qr分解

A = Q * R

// 为了方便展示，将名类型替换为真实类型
// 第一个函数sgeqrf,用于求解R矩阵，以及生成Q矩阵的中间因子
// https://netlib.org/lapack/explore-html//d0/da1/group__geqrf_ga01418ec8a343078109cde4febb6cda7d.html
int sgeqrf_(
        int* m,//行数
        int* n,//列数
        float* a,//输入: 待分解矩阵A的m*n矩阵. 输出:上半部分存R,下半部分存用于构造Q的householder向量
        int* lda,//矩阵A的主维度
        float* tau,//输出 长度为min(m,n)的数组,存储householder变换的标量因子
        float* work,//工作空间,临时缓冲区,存储中间计算结果
        int* lwork,//work数组的大小,若为-1,则函数会进行工作空间查询,返回最优大小
        int* info);//输出,状态标志
// 第二个函数sorgqr,用于生成矩阵Q
// https://netlib.org/lapack/explore-html/d4/dfc/group__ungqr_gacb2375e4540a08ea3c59dcfea6701ca2.html
int sorgqr_(
        int* m,// Q的行数
        int* n,// 想要生成Q的列数.通常N=M或N=K(生成经济型矩阵)
        int* k,// 变换的阶数,通常等于sgeqrf_中的min(M,N)
        float* a,// 输入: sgeqrf_处理后的矩阵A(包含householder向量)  输出:生成的显式Q矩阵(会覆盖输入内容)
        int* lda,// 矩阵A的主维数
        float* tau,// sgeqrf_生成的标量因子数组
        float* work,// 工作空间,临时缓冲区
        int* lwork,// work数组大小
        int* info);// 状态码. 0成功

// 为了方便展示，将名类型替换为真实类型
// 函数路径：/opt/homebrew/opt/openblas/include/lapacke.h
// 第一个函数，用于求解R矩阵，以及生成Q矩阵的中间因子
int LAPACKE_sgeqrf( int matrix_layout, // 内存布局，行优先：LAPACK_ROW_MAJOR  列优先：LAPACK_COL_MAJOR
    int m, // 输入矩阵A的行数
    int n,// 输入矩阵A的列数
    float* a, // 输入矩阵A，同时也是输出矩阵。上三角部分为R矩阵，下半部分存用于构造Q的householder向量
    int lda, // 输入矩阵A的leading dimension.如果为行优先，即为列数。如果为列优先，即为行数
    float* tau ); // 输出值。长度为min(m,n)的数组,存储householder变换的标量因子
 // 第二各函数，用于生成Q矩阵
 int LAPACKE_sorgqr( int matrix_layout, // 内存布局，行优先：LAPACK_ROW_MAJOR  列优先：LAPACK_COL_MAJOR  
     int m, // Q的行数
     int n, // Q的列数
     int k, // 变换的阶数，通常为min(m,n)
     float* a, // 输入: LAPACKE_sgeqrf处理后的矩阵A(包含householder向量)  输出:生成的显式Q矩阵(会覆盖输入内容)
     int lda, // 矩阵Q的leading dimension。行优先即为列数，列优先即为行数
     const float* tau);// 输入值，LAPACKE_sgeqrf生成的中间因子

#include <algorithm>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
// 注意：本节所用的QR分解只需包含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;
  }
}
/**
    QR分解公式：A = Q*R
*/
int main(int argc, char* argv[]) {
  int M = 3;
  int N = 2;
  /**列优先，按照行展示
  A:
   12 -51
   6 167
   -4 24
  */
  // 1.构造输入数据
  vector<float> A = {12, 6, -4, -51, 167, 24};  // 内存布局使用列优先
  vector<float> tau(std::min(M, N), 0);  // 中间变换因子，长度为min(M,N)的数组
  // 2. 执行第一步分解，拿到上三角矩阵R
  LAPACKE_sgeqrf(
      LAPACK_COL_MAJOR,  // 列优先存储
      M, N,              // 输入矩阵A的行和列
      A.data(),  // 输入矩阵A的原始数据，输出的R矩阵会存放在A的上三角矩阵
      M,  // 输入矩阵A的leading-dimension, 因为为列优先，即为行数3
      tau.data());  // 中间变换因子，用于继续求解Q矩阵
  // 3.
  // 打印R矩阵。A矩阵分为两部分，上三角矩阵为R,下三角矩阵为生成Q矩阵所需的数据。
  print_col_matrix_by_row("R", A.data(), M, N);
  // 4. 以下为生成Q矩阵的过程
  vector<float> Q(M * M);  // 生成完全Q矩阵即为m*m的大小。
  std::memcpy(Q.data(), A.data(),
              M * N * sizeof(float));  // sgeqrf处理后的矩阵A
  LAPACKE_sorgqr(LAPACK_COL_MAJOR,     // 内存布局使用列优先
                 M, M,                 // 矩阵Q的行数和列数
                 M,  // 矩阵Q的leading-dimension，因为是列优先，所以是行数3
                 Q.data(),  // 输入矩阵，同时也是输出矩阵
                 M,  // 输入矩阵的leading-dimension,因为是列优先，所以是行数3
                 tau.data());  // 中间变换因子，为上一步sgeqrf求得的结果
  // 5. 打印结果Q矩阵
  print_col_matrix_by_row("Q", Q.data(), M, M);
  // 6. 以上即为QR分解的全过程，下面为重建验证的部分
  vector<float> C(M * N);
  // 7. 使用Q*R重建矩阵A,观察是否等于输入矩阵A
  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, 1.0, Q.data(),
              M, A.data(), M, 1.0, C.data(), M);
  print_col_matrix_by_row("重建A", C.data(), M, N);
  return 0;
}

// 编译链接
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include qr_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o qr_main
// 执行
./qr_main
//执行结果
R:3x2
-14     -21
0.230769        -175
-0.153846       0.0555556
Q:3x3
-0.857143       0.394286        0.331429
-0.428571       -0.902857       -0.0342857
0.285714        -0.171429       0.942857
重建A:3x2
12.04   -50.9816
5.79692 166.998
-4.18462        24.0524