openblas系列之ssyev

A = Q * Λ * Qᵀ

// 找到了原始fortran语言的函数定义。https://netlib.org/lapack/explore-html/d8/d1c/group__heev_gab3474da7a9f45653c4b7bc7d523ebc61.html
// 为了方便展示，使用原始变量类型
int ssyev_(
    const char* jobz, // N只计算特征值. V计算特征值和特征向量
    const char* uplo, // 更新上三角或者下三角(U/L)
    int* n, // 矩阵维度
    float* a, // 矩阵A
    int* lda, // 矩阵的 leading dimension
    float* w, // 输出的特征值数组
    float* work, // 输出的特征向量数组
    int* lwork, // 特征向量数组长度
    int* info);// 返回的状态码 0成功完成 大于0有非法值 小于0算法未能收敛

// 函数定义在:/opt/homebrew/opt/openblas/include/lapacke.h
// 为了方便展示，将别名类型替换为真实类型
// A = Q * Λ * Qᵀ
int LAPACKE_ssyev(
    int matrix_layout, // 行优先或者列优先。行优先：LAPACK_ROW_MAJOR  列优先：LAPACK_COL_MAJOR
    char jobz, // 是否计算特征向量。 N只计算特征值. V计算特征值和特征向量
    char uplo, // 上三角矩阵或者下三角矩阵，只需提供上三角或者下三角即可，其他元素为0。上三角矩阵：U，下三角矩阵：L
    int n, // 结果矩阵的阶数
    float* a, // 输入矩阵A,计算完成后，如果需要计算特征向量，将会被特征向量矩阵Q原地覆盖
    int lda, // 矩阵A的leading dimension。 根据行优先或者列优先来定。如果是行优先，即为列数。如果是列优先，即为行数
    float* w); // 计算出的特征值Λ，是个数组，只保存特征对角线上的特征值，大小为对角线的长度

#include <iostream>
#include <string>
#include <vector>
// 注意：本节的函数定义是在lapacke.h包内，引入cblas.h是为了重建验证
#include <lapacke.h>
#include <cblas.h>

using std::string;
using std::cout;
using std::endl;
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<<"\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 = Q * Λ * Q^T
*/
int main(int argc,char* argv[]){
    int N = 2;
    /** 列优先，注释为行展示
    A:
        16030 6955
        0 3809
    */
    // 1. 构造数据
    vector<float> A = {16030,0,6955,3809}; // 列优先存储。输入矩阵，计算完成后被替换为矩阵Q,即为求得的所有特征向量。
    vector<float> W(N,0); // 求得的特征值，由于中间的特征矩阵Λ是个对角矩阵，这里只用数组存储了对角线元素，节省空间。
    // 2. 计算特征值和特征向量
    LAPACKE_ssyev(
        LAPACK_COL_MAJOR, // 列优先存储
        'V', // V：既要计算特征值，也计算特征向量
        'U', // U：上三角矩阵
        N, // 结果矩阵Q的阶数
        A.data(), // 矩阵A的数据
        N, // 矩阵A的leading dimension,由于这里是列优先存储，即为行数2
        W.data());// 求得的特征值
    // 3. 打印结果
    cout<<"所有的特征值Λ:";
    for (auto& item : W) {
        cout<<item<<" ";
    }
    cout<<endl;
    print_col_matrix_by_row("特征向量Q:", A.data(), N, N);
    // 以上就已经完成了求解的过程，下面是利用以前的所学的内容进行重建验证。
    // 4.重建矩阵A
    // 4.1 先把特征值数组还原为对角矩阵
    vector<float> M = {W[0],0,0,W[1]};// 列优先存储
    // 4.1 进行两次矩阵乘法,先计算Q*Λ将结果存储C中
    vector<float> C(4,0);
    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, N, N, 1.0, A.data(), N, M.data(), N, 0.0,C.data(),N);
    // 4.2 再计算C*Q^T将结果存储D中，打印结果矩阵D
    vector<float> D(4,0);
    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasTrans, N, N, N, 1.0, C.data(), N, A.data(), N, 0.0, D.data(), N);
    print_col_matrix_by_row("重建矩阵A:", D.data(), N, N);
    return 0;
}

# 编译链接
g++ -std=c++17 -I/opt/homebrew/opt/openblas/include ssyev_main.cpp -L/opt/homebrew/opt/openblas/lib -lopenblas -o ssyev_main
# 执行
./ssyev_main
# 执行结果
所有的特征值Λ:661.518 19177.5 
特征向量Q::2x2
0.412295 -0.91105 
-0.91105 -0.412295 
重建矩阵A::2x2
16030 6955 
6955 3809