Медленное выполнение программы C/C++ c библ. LAPACK, Eigen

[ValfennayaVaflaya]

Здравствуйте, есть задача, решаемая методом конечных элементов, которая сводится к обобщенной задаче собственных значений $Ax=\lambda Bx$ с трехдиагональными, симметричными матрицами. Пишу на $C/C++$ в Visual Studio 2022, для решения используется LAPACK, а именно ф-я: LAPACKE_dggev, а также библиотека Eigen (GeneralizedEigenSolver).
Мой вопрос связан с тем, что вычисления для матриц 2000 $\times$ 2000 в первом случае (CLAPACK) заняли 40 (!) минут, во втором случае (Eigen) больше часа. У меня это вызывает беспокойство, учитывая, что на Python решение данной задачи занимает несколько минут. Полученные собственные числа совпадают в обоих случаях и похожи на правду. Я с таким сталкиваюсь впервые и поэтому даже не знаю в какую сторону двигаться, чтобы искать решение, поэтому буду рад любым версиям, объясняющим такое замедление.

[mihaild]

Для начала - проверьте, что библиотеки собраны с оптимизацией. И что долго выполняется именно эта функция (простейший вариант - позвать её дважды).
И 40 минут - это одно вычисление, или всё вместе?

[ValfennayaVaflaya]

1) Для LAPACK я использую готовое решение для VS отсюда https://icl.utk.edu/lapack-for-windows/lapack/index.html#lapacke, там все уже собрано, так что не могу ответить.
2) Да, проблема явно в ней, в коде больше ничего нет, только заполнение этих матриц числами, просто арифметические операции, printf перед и после нее это вроде подтверждают.
3) 40 минут - это на всю программу, но, как писал выше, больше код ничего производительного не делает, так что такое время уходит только на вычисления.

[mihaild]

Т.е. один вызов занимает 40 минут?
Проверьте всё-таки, что будет если сделать несколько вызовов (можно на матрице меньшего размера, чтобы не ждать часами), это дешевый способ проверить что действительно ничего не перепутано.
(ну и желательно покажите релевантную часть кода на плюсах и питоне)

[ValfennayaVaflaya]

Код небольшой, поэтому прикреплю весь, там решается ур-е Шредингера для одномерной ямы методом конечных элементов. Есть две матрицы, они собираются стандартным образом, затем вызывается ф-я библиотечная, выводятся собственные числа отрицательные, а в файл пишутся собственные ф-ии. На питоне, к сожалению, кода нет, так как он не мой, но сам видел, что задача не должна решатся 40 минут. Да, 40 минут на один вызов выходит.

Используется синтаксис C

#include <stdlib.h>
#include <stdio.h>
#include <cmath>
#include "lapacke.h"

/* Auxiliary routines prototypes */
extern void print_matrix(char* desc, lapack_int m, lapack_int n, double* a, lapack_int lda);
extern void print_matrix_int(char* desc, lapack_int m, lapack_int n, lapack_int* a, lapack_int lda);
extern void print_vector(char* desc, lapack_int n, double* a);
extern void print_vector_int(char* desc, lapack_int n, lapack_int* a);

extern double a_f(lapack_int k, lapack_int i, lapack_int j);
extern double b_f(lapack_int k, lapack_int i, lapack_int j);

extern void eig_values();

/* Parameters */
#define NN 500          // number of elements
#define M (NN+1)         // number of nodes
#define A 10.            // radius
#define h (2*A/NN)       // lattice step
#define V0 12            // extern potential (default: 12)

double MainMatrix[M * M];
double RightSideMatrix[M * M];
double X[M];
double V[M - 1];                      // Potential
lapack_int L[2 * (M - 1)];            // Matrix of the Nodes 

lapack_int INFO = 0;

double ALPHAR[M], ALPHAI[M], BETA[M]; 
double VL[M * M], VR[M * M];
double EIG[M];

/* Main program */
int main() {

    //----------------------------------------------------------------------------------//

    for (int i = 0; i < M; ++i) { X[i] = h * i - A; }

    for (int i = 0; i < NN; ++i) {
        if (abs((X[i] + X[i + 1]) / 2) < 1)
            V[i] = -V0;
        else V[i] = 0;
    }

    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < M; ++j) {
            MainMatrix[M * i + j] = 0;
            RightSideMatrix[M * i + j] = 0;
        }
    }
    
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < M; ++j)
            L[i * NN + j] = j + i;

    //print_matrix_int("L\n", 2, NN, L, NN);

    for (int k = 0; k < M - 1; ++k)
        for (int i = 0; i < 2; ++i)
            for (int j = 0; j < 2; ++j) {
                MainMatrix[L[i * NN + k] * M + L[j * NN + k]] += a_f(k, i, j);
                RightSideMatrix[L[i * NN + k] * M + L[j * NN + k]] += b_f(k, i, j);
            }

    // Correcting the equations (first and last) 

    for (int j = 0; j < M; ++j) {
        MainMatrix[M * 0 + j] = 0;
        MainMatrix[M * (M - 1) + j] = 0;
        RightSideMatrix[M * 0 + j] = 0;
        RightSideMatrix[M * (M - 1) + j] = 0;
    }

    MainMatrix[0] = 1;
    MainMatrix[M * (M - 1) + M - 1] = 1;

    //print_matrix("Main Matrix\n", M, M, MainMatrix, M);

    //print_matrix("Right Matrix\n", M, M, RightSideMatrix, M);

    printf("Before LAPACKE function\n");

    INFO = LAPACKE_dggev(LAPACK_ROW_MAJOR, 'N', 'V', M, MainMatrix, M, RightSideMatrix, M, ALPHAR, ALPHAI, BETA, VL, M, VR, M);
    printf("After LAPACKE function\n");
    
    eig_values();

    if (INFO == 0) printf("Ok\n");
    else printf("INFO = %d\n", INFO);

    FILE* file = fopen("FEM.txt", "w+");

    printf("Negative eigenvalues: ");
    for (int i = 0; i < M; ++i) {
        fprintf(file, " %6.2f", X[i]);
        for (int j = 0; j < M; ++j) {
            if (EIG[j] < 0) {                                 // write only eigenvectors with NEGATIVE eigenvalues...
                fprintf(file, " %6.2f", VR[i * M + j]);
            }
        }
        if (EIG[i] < 0) {
            printf(" %8.4f", EIG[i]);
        }
        fprintf(file, "\n");
    }
    printf("\n");

    fclose(file);

    //----------------------------------------------------------------------------------//

    exit(0);
}

void print_matrix(char* desc, lapack_int m, lapack_int n, double* a, lapack_int lda) {
    lapack_int i, j;
    printf("\n %s\n", desc);
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) printf(" %6.2f", a[i * lda + j]);
        printf("\n");
    }
}

void print_matrix_int(char* desc, lapack_int m, lapack_int n, lapack_int* a, lapack_int lda) {
    lapack_int i, j;
    printf("\n %s\n", desc);
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) printf(" %6d", a[i * lda + j]);
        printf("\n");
    }
}

void print_vector(char* desc, lapack_int n, double* a) {
    lapack_int j;
    printf("\n %s\n", desc);
    for (j = 0; j < n; j++)
        printf(" %6.2f", a[j]);
    printf("\n");
}

void print_vector_int(char* desc, lapack_int n, lapack_int* a) {
    lapack_int j;
    printf("\n %s\n", desc);
    for (j = 0; j < n; j++)
        printf(" %6i", a[j]);
    printf("\n");
}

double a_f(lapack_int k, lapack_int i, lapack_int j) {
    if (i == j) {
        return 1. / h + (2. / 3) * V[k] * h;
    }
    else
        return -1. / h + (1. / 3) * V[k] * h;
}

double b_f(lapack_int k, lapack_int i, lapack_int j) {
    if (i == j) {
        return (+2. / 3) * h;
    }
    else
        return (+1. / 3) * h;
}

void eig_values() {
    for (int i = 0; i < M; ++i) {
        EIG[i] = ALPHAR[i] / BETA[i];
    }
}
 

[mihaild]

На моей машине этот код отрабатывает за 0.6с и выдает -11.1520 -8.6534 -4.6963 -0.2158 - такие числа и должны быть?
Если да, то проблемы почти наверняка в сборке lapack.

ValfennayaVaflaya в сообщении #1589974 писал(а):

for (int i = 0; i < 2; ++i)
for (int j = 0; j < M; ++j)
L[i * NN + j] = j + i;

А вот тут у Вас UB. Вряд ли это влияет в данном случае, но всё равно нехорошо.

[zykov]

mihaild в сообщении #1589978 писал(а):

А вот тут у Вас UB. Вряд ли это влияет в данном случае, но всё равно нехорошо.

Там за границы массива вылезает при $i=1$ и $j=M-1$.
Сам массив имеет размер $2(M-1)$, последний индекс $2(M-1)-1$.

[ValfennayaVaflaya]

Да, числа такие же, но у меня вычисления занимают 10 секунд. Ошибка с выходом за пределы массива случайно получилась при копировании кода сюда, просто менял названия переменных для удобства и стер единицу. Хорошо, тогда буду разбираться как собирать самому lapack, спасибо.

[zgemm]

пользуйте DSBGV и будет вам в много порядков меньше арифметическая сложность