LAPACKE 函数中用于对角化的全矩阵或三角部分？

Question

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

根据文档 https://software.intel.com/en-us/mkl-developer-reference-c-syev，我得出的结论是我只需要传递矩阵的 upper/lower 三角形部分。它在那里说了它 "is an array containing either upper or lower triangular part of the symmetric matrix 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

根据this example，LAPACKE_dsyev似乎需要在完整存储方案中存储上三角或下三角部分。甚至文档也说，参数 a 应该是一个大小为 max(1, lda * n) 的数组，在您的情况下为 16。

当我切换到这个矩阵定义时：

double matrix_ex_uppertri[16] =
   { -2.0,  0.0,  0.5,  0.0,
      0.0,  0.5, -2.0,  0.5,
      0.0,  0.0,  0.5,  0.0,
      0.0,  0.0,  0.0, -1.0 };

然后，我得到了正确的特征值。

需要承认文档在这里非常具有误导性。不明白为什么他们不明确指定存储方案。

Answer 2

根据documentation of Lapack regarding naming scheme, letters sy in the name of function dsyev() refer to a symmetric matrix, while letters d refers to double precision and ev refers to eigenvalues. Nevertheless, The format sy corresponds to a conventional storage作为一个二维数组，其形状与矩阵一致。根据参数 UPLO 的值使用上三角部分或下三角部分。

要使用打包格式，请查看 LAPACKE 中的函数 dspev(), where the letters sp corresponds to packed storage of symmetric matrices. The function LAPACKE_dspev()，它提供了一个方便的 C 接口。

这里是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,...).

的结果保持一致

LAPACKE 函数中用于对角化的全矩阵或三角部分？

Full matrix or triagonal part necessary in LAPACKE function for diagonalization?

c++

lapack

intel-mkl