Ошибка сегментации, когда свободный указатель MKL после вызова подпрограммы lapacke hesv

Я разрабатываю код для инвертирования матрицы с использованием интерфейса lapacke для языка C. Я использовал библиотеку intel mkl для выполнения такой работы. Тем не менее, когда я пробую простой тест инвертировать матрицу переменного размера N, я получаю неопределенное поведение, когда освобождаю указатель матрицы. Что странно, так это то, что он работает, если N = 3, но перестает работать с N >= 4, например, если я не освобождаю указатель, код работает отлично для любого значения N. Другие указатели, требуемые функцией, я могу освободить без каких-либо беда. Любые идеи?

Дополнительная информация: я использую функцию zhesv, которая может решать несколько линейных систем. API можно найти по адресу: https://software.intel.com/en-us/node/520994

Код, который я пытаюсь запустить:

/*
 * compile with intel mkl installation:
 *
 * gcc -o test exe/lapack_inversion_test.c  -L${MKLROOT}/lib/intel64 -Wl,--no-as-    needed -lmkl_intel_ilp64 -lmkl_gnu_thread -lmkl_core -lgomp -lm -ldl
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include <math.h>

#define PI 3.141592653589793 // to define matrix entries

#define LAPACK_ROW_MAJOR 101
#define LAPACK_COL_MAJOR 102

int main(int argc, char * argv[])
{
    int j, // counter
        i, // counter
        N, // The size of the Matrix
        k;

    double arg;

    sscanf(argv[1], "%d", &N);

    int * ipiv = (int *) malloc(N * sizeof(int));

    MKL_Complex16 x; x.real = 0; x.imag = 0;
    MKL_Complex16 * A = malloc(N * N * sizeof(MKL_Complex16));
    MKL_Complex16 * Acopy = malloc(N * N * sizeof(MKL_Complex16));
    MKL_Complex16 * Id = malloc(N * N * sizeof(MKL_Complex16));

    for (i = 0; i < N; i++)
    {   // Row major storage
        A[i * N + i].real = 1.5 + sin( 3 * PI * i / N );
        A[i * N + i].imag = 0;
        Acopy[i * N + i].real = 1.5 + sin( 3 * PI * i / N );
        Acopy[i * N + i].imag = 0;
        Id[i * N + i].real = 1;
        Id[i * N + i].imag = 0;
        for (j = 0; j < i; j++)
        {
            arg = 10 * PI * ((double) i * j) / (N * N);
            A[i * N + j].real = 2 * sin(arg) + 2;
            A[i * N + j].imag = 3 * cos(arg);
            A[j * N + i].real = 0; // Does not matter the upper tirangular
            A[j * N + i].imag = 0; // when call lapack routine with 'L'

            Acopy[i * N + j].real = 2 * sin(arg) + 2;
            Acopy[i * N + j].imag = 3 * cos(arg);
            Acopy[j * N + i].real = Acopy[i * N + j].real;
            Acopy[j * N + i].imag = - Acopy[i * N + j].imag;

            Id[i * N + j].real = 0; // Identity
            Id[i * N + j].imag = 0;
            Id[j * N + i].real = 0;
            Id[j * N + i].imag = 0;
        }
    }

    i = LAPACKE_zhesv(LAPACK_ROW_MAJOR, 'L', N, N, A, N, ipiv, Id, N);
    printf("\n\nLapacke returned : %d\n", i);

    printf("\n\nIf there was any problem print identity: \n");

    for (i = 0; i < N; i++)
    {
        printf("\n\t|");
        for (j = 0; j < N; j++)
        {
            x.real = 0;
            x.imag = 0;
            for (k = 0; k < N; k++)
            {
                x.real += Id[i*N + k].real * Acopy[k*N + j].real;
                x.real -= Id[i*N + k].imag * Acopy[k*N + j].imag;
                x.imag += Id[i*N + k].real * Acopy[k*N + j].imag;
                x.imag += Id[i*N + k].imag * Acopy[k*N + j].real;
            }
            printf(" (%6.3lf,%6.3lf) |", x.real, x.imag);
        }
    }

    // free(A); if one try to free does not work with N >= 4 !!!!
    free(Id);
    free(ipiv);

    printf("\n\n");
    return 0;
}

Если мы позволим коду освободить (A) (раскомментировать), он сделает:

$ ./тест 3

Лапак вернулся : 0

Если возникли проблемы с печатью удостоверения:

| ( 1.000,-0.000) | ( 0.000,-0.000) | (-0.000,-0.000) |
| ( 0.000, 0.000) | ( 1.000, 0.000) | ( 0.000, 0.000) |
| (-0.000,-0.000) | (-0.000, 0.000) | ( 1.000, 0.000) |

$ ./тест 4

Лапак вернулся : 0

Если возникли проблемы с печатью удостоверения:

| ( 1.000,-0.000) | ( 0.000,-0.000) | ( 0.000, 0.000) | (-0.000, 0.000) |
| (-0.000, 0.000) | ( 1.000, 0.000) | ( 0.000, 0.000) | (-0.000, 0.000) |
| (-0.000, 0.000) | (-0.000, 0.000) | ( 1.000,-0.000) | (-0.000,-0.000) |

Ошибка сегментации (дамп ядра)


person Alex Andriati    schedule 30.08.2018    source источник
comment
Значение i (возвращаемое из LAPACKE_zhesv) всегда равно 0?   -  person P.W    schedule 30.08.2018
comment
Не воспроизводимо. У меня работает от N=4 до N=8.   -  person P.W    schedule 30.08.2018
comment
Да, возвращаемое значение равно нулю, если можно выполнить инверсию.   -  person Alex Andriati    schedule 30.08.2018


Ответы (1)


Давайте посмотрим на объявление LAPACKE_zhesv:

lapack_int LAPACKE_zhesv(int matrix_layout, char uplo, lapack_int n, lapack_int nrhs, lapack_complex_double* a, lapack_int lda, lapack_int* ipiv, lapack_complex_double* b , lapack_int ldb);

В частности, он хочет lapack_int* ipiv, но вы объявляете его как int* ipiv. Затем вы связываете свою программу с mkl_intel_ilp64 и не определять MKL_ILP64. Ваша программа и MKL в конечном итоге используют несовместимые целочисленные типы: sizeof(int) = 4, тогда как sizeof(lapack_int) = 8.

Быстрое решение - связать с mkl_intel_lp64. Но вы должны переписать свой код, чтобы использовать lapack_int вместо int, тогда он будет корректно работать как с LP64, так и с ILP64.

person Evg    schedule 30.08.2018
comment
Большое спасибо, это сработало, просто изменив параметры компиляции. Я сейчас использую ссылку -lmkl_intel_lp64 вместо _ilp64. Таким образом, он работает как с lapack_int, так и с int *. - person Alex Andriati; 30.08.2018
comment
@AlexAndriati, если MKL_ILP64 не определено, то это фактически #define lapack_int int в заголовках MKL. - person Evg; 30.08.2018