c++ - LAPACKE 函数中对角化所需的完整矩阵或三角部分？

Question

我想计算对称矩阵的特征值，并希望使用 C++ 中英特尔 MKL 库中的 LAPACKE_dsyev 函数来计算。但我对矩阵需要传递的形式有点困惑。

来自文档https://software.intel.com/en-us/mkl-developer-reference-c-syev ，我得出的结论是我必须仅通过矩阵的上/下三角形部分。它说关于它“是一个包含对称矩阵 A 的上三角部分或下三角部分的数组”的论点。然而，似乎实际上需要将指向完整矩阵的指针传递给例程。

假设我想对角化以下矩阵:

[[-2,   0,     0.5,  0], 
[0,     0.5, -2,     0.5], 
[0.5,  -2,    0.5,  0], 
[0,     0.5,  0,    -1]], 
which has eigenvalues [ 2.56,  -2.22, -1.53, -0.81]

然后在下面的代码中，只有第一个选项给出了正确的值。

#include <iostream>
#include"mkl_lapacke.h"

using namespace std;

int main(){
    MKL_INT N = 4;
    // use full matrix
    double matrix_ex_full[16] = {-2,0,0.5,0,0,0.5,-2,0.5,0.5, -2, 0.5, 0, 0,0.5,0,-1};
    double evals_full[4];
    MKL_INT test1 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', N, matrix_ex_full,N, evals_full);
    cout << "success = " <<test1 << endl;
    for (MKL_INT i = 0;i<4;i++)
        cout << evals_full[i] << endl;
    // use upper triagonal only
    double matrix_ex_uppertri[10] = {-2, 0, 0.5, 0, 0.5, -2, 0.5, 0.5, 0, -1};
    double evals_uppertri[4];
    MKL_INT test2 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', N, matrix_ex_uppertri,N, evals_uppertri);
    cout << "success = " <<test2 << endl;
    for (MKL_INT i = 0;i<4;i++)
        cout << evals_uppertri[i] << endl;
    // (Compiled with g++ test.cpp -o main -m64 -I/share/opt/intel/compilers_and_libraries_2018.0.128/linux/mkl/include -L/share/opt/intel/compilers_and_libraries_2018.0.128/linux/mkl/lib/intel64 -Wl,--no-as-needed -lmkl_intel_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl)
}

为什么我只有在将指针传递到完整矩阵时才能获得正确的特征值？

我觉得这一定是一个微不足道的问题，但我错过了什么？如果始终必须指定完整矩阵，那么用“U”或“L”指定给出矩阵的哪个三角形对我来说是没有意义的。还是我在其他地方做错了什么？

非常感谢您的帮助!

Answer 1

根据documentation of Lapack regarding naming scheme ，函数dsyev()名称中的字母sy指的是对称矩阵，而字母d指的是 double ，ev 指的是特征值。尽管如此，格式 sy 对应于 a conventional storage作为形状与矩阵一致的二维数组。根据参数 UPLO 的值使用上三角部分或下三角部分。

要使用打包格式，请查看函数 dspev() , 其中字母 sp 对应于 packed storage of symmetric matrices 。函数LAPACKE_dspev() LAPACKE 提供了一个方便的 C 接口(interface)。

这是由g++ main.cpp -o main -llapacke -llapack -lblas -lm -Wall编译的示例代码:

#include <iostream>
#include <math.h>

extern "C" { 
#include <lapacke.h>
}

using namespace std;

int main(int argc, char *argv[])
{
    lapack_int N = 4;
    // use full matrix
    double matrix_ex_full[16] = {-2,0,0.5,0,    0,0.5,-2,0.7,    0.5, -2, 0.5, 0,     0,0.7,0,-1};
    double evals_full[4];
    lapack_int test1 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', N, matrix_ex_full,N, evals_full);
    cout << "success = " <<test1 << endl;
    for (int i = 0;i<4;i++)
        cout << evals_full[i] << endl;
    // use upper triagonal only
    double matrix_ex_uppertri[10] = {-2, 0, 0.5, 0,    0.5, -2, 0.7,    0.5, 0,   -1};
    double evals_uppertri[4];


    int test2 = LAPACKE_dspev(LAPACK_COL_MAJOR, 'N', 'L', N, matrix_ex_uppertri, evals_uppertri,NULL,N);
    cout << "success = " <<test2 << endl;
    for (int i = 0;i<4;i++)
        cout << evals_uppertri[i] << endl;
    return 0;
}

如 LAPACKE_dspev_work() 的文档所示，利用 LAPACK_COL_MAJOR 可以节省临时数组的一些额外分配。由于未计算特征向量，因此结果与 LAPACKE_dsyev(LAPACK_ROW_MAJOR,...) 的结果保持一致。