Полная матрица или трехугольная часть, необходимая в функции LAPACKE для диагонализации? - PullRequest
Новая
       85

Полная матрица или трехугольная часть, необходимая в функции LAPACKE для диагонализации?

3 голосов
/ 01 октября 2019

Я хотел бы вычислить собственные значения симметричной матрицы и хотел бы использовать для этого функцию LAPACKE_dsyev из библиотеки Intel MKL в C ++. Но я немного озадачен относительно формы, в которой необходимо передать матрицу.

Из документации 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', какой треугольник матрицы дан. Или я что-то не так делаю в другом месте?

Большое спасибо за вашу помощь!

Ответы [ 2 ]

3 голосов
/ 02 октября 2019

Согласно документации Лапака относительно схемы именования , буквы sy в названии функции dsyev() относятся к симметричной матрице, а буквы d относятся к двойной точности и evотносится к собственным значениям. Тем не менее, Формат sy соответствует обычному хранилищу в виде двумерного массива формы, соответствующей форме матрицы. В зависимости от значения аргумента UPLO используется либо верхняя треугольная часть, либо нижняя часть.

Чтобы использовать упакованный формат, взгляните на функцию dspev(),, где буквы sp соответствуют упакованному хранилищу симметричных матриц . Функция LAPACKE_dspev() от LAPACKE предоставляет удобный интерфейс для 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,...).

1 голос
/ 01 октября 2019

Согласно этому примеру , кажется, что 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 };

тогда я получил правильные собственные значения.

Нужно признать, что документация здесь вводит в заблуждение. Не понимаю, почему они не указывают четко схему хранения.

...