Решение линейного уравнения с помощью LAPACKE

Question

Я пытаюсь решить линейное уравнение (симметричное, трехугольное и положительное). Я должен использовать LAPACKE. Мой код выглядит следующим образом:

#include <lapacke.h>
#include <stdio.h>



void print_mtrx(double * mtrx, int n, int m)
{
    int i, j;

    for(i = 0; i < n; i++)
    {
        for(j = 0; j < m; j++)
        {
            printf("%f ", mtrx[i*m+j]);
        }
        printf("\n");
    }
    printf("\n");
}

int main()
{
    double matrix[5*5] = {
        2,  0,  0,  0,  0,
        0,  2,  0,  0,  0,
        0,  0,  2,  0,  0,
        0,  0,  0,  2,  0,
        0,  0,  0,  0,  2
    };

    double rozw[5] = {1,2,3,4,5};

    double matrix2[5*5] = {
        7,  0,  0,  0,  0,
        0,  7,  0,  0,  0,
        0,  0,  7,  0,  0,
        0,  0,  0,  7,  0,
        0,  0,  0,  0,  7
    };



    LAPACKE_dptsv(LAPACK_COL_MAJOR, 5, 5, matrix, matrix2, rozw, 5);

    print_mtrx(matrix, 5, 5);
    print_mtrx(matrix2, 5, 5);
    print_mtrx(rozw, 5, 1);

}

Функция LAPACKE, похоже, ничего не делает, без каких-либо ошибок. Основная проблема в том, что я понятия не имею, что означают параметры функции. Я долго искал, но нет реальной документации. Вот что мне удалось найти или угадать:

• int matrix_order - LAPACK_COL_MAJOR или LAPACK_ROW_MAJOR, как матрица представлена ​​в памяти
• lapack_int n - размер матрицы (т.е. количество столбцов)
• lapack_int nrhs - не уверен, возможно, размер вектора b
• double * d - матрица уравнений
• двойной * е - понятия не имею.
• double * b - вектор решений уравнений из d
• lapack_int ldb - ведущее направление b (почему? оно не идентично nrhs, которое само по себе идентично n?)

Как я могу найти реальное значение этих аргументов? Как я могу заставить мой код работать?

Answer 1

Когда дело доходит до документации для BLAS и / или LAPACK, Intel, вероятно, является наиболее полной. Вы можете найти документацию для ? Ptsv, которая объясняет, для чего предназначен каждый параметр.

(Подсказка: при поиске BLAS или LAPACK в Google обязательно s/d/c/z префикс.)

Вот соответствующий фрагмент:

Рутина решает для X действительная или сложная система линейных уравнений A*X = B, где A является n-от-n симметричная / эрмитова положительно-определенная трехдиагональная матрица, столбцы матрицы B отдельные правые части, а столбцы X соответствующие решения.

A учитывается как A = L*D*LT (настоящие ароматы) или A = L*D*LH (сложные ароматы), и факторизованная форма A Затем используется для решения системы уравнений A*X = B,

Входные параметры

n: Порядок матрицы A; n ≥ 0,

nrhs: Количество правых частей, количество столбцов в B; nrhs ≥ 0,

d: Массив, размер как минимум max(1, n), Содержит диагональные элементы трехдиагональной матрицы A,

e, b: Массивы: e(n - 1), b(ldb,*), Массив e содержит (n - 1) субдиагональные элементы A, Массив b содержит матрицу B чьи столбцы являются правыми частями для систем уравнений. Второе измерение b должен быть не менее max(1,nrhs),

ldb: Ведущее измерение b; ldb ≥ max(1, n),

Выходные параметры

d: Перезаписано n диагональные элементы диагональной матрицы D от L*D*LT (реальный) / L*D*LH (сложная) факторизация A,

e: Перезаписано (n - 1) субдиагональные элементы единичного бидиагонального фактора L от факторизации A,

b: Перезаписывается матрицей решения X,

info: Если info = 0, выполнение успешно. Если info = -i, iпараметр имеет недопустимое значение. Если info = i, ведущий несовершеннолетний порядка i (и, следовательно, матрица A само по себе) не является положительно определенным, и решение не было вычислено. Факторизация не была завершена, если i = n,

Answer 2

Поэтому необходимо изучить документацию по чистому LAPACK ( http://www.netlib.org/lapack/explore-html/d0/dea/dptsv_8f.html), а также игнорировать неправильное (в случае LAPACKE) замечание, что LDB будет больше или равно max(1,N).

Правильная программа выглядит следующим образом:

#include <lapacke.h>
#include <stdio.h>



void print_mtrx(double * mtrx, int n, int m)
{
    int i, j;

    for(i = 0; i < n; i++)
    {
        for(j = 0; j < m; j++)
        {
            printf("%f ", mtrx[i*m+j]);
        }
        printf("\n");
    }
    printf("\n");
}

int main()
{
    double diagonal[5] = {5,1,5,1,5};
    double subdiagonal[4] = {0,0,0,0};

    double solution[5] = {1,2,3,4,5};


    LAPACKE_dptsv(LAPACK_ROW_MAJOR, 5 /*size of matrix*/, 1 /*number of columns in solution*/,
                  diagonal, subdiagonal, solution, 1 /*leading dimension of solution vector*/);

    print_mtrx(solution, 5, 1);
}