Подтвердить что ты не робот

Вычисление обратной матрицы, использующей lapack в C

Я хотел бы иметь возможность вычислить инверсию общей матрицы NxN в C/С++, используя lapack.

Мое понимание заключается в том, что способ сделать инверсию в лапаке с помощью функции dgetri, однако, я не могу понять, что все его аргументы должны быть.

Вот код, который у меня есть:

void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);

int main(){

    double M [9] = {
        1,2,3,
        4,5,6,
        7,8,9
    };

    return 0;
}

Как бы вы завершили его, чтобы получить обратную матрицу M 3x3, используя dgetri _?

4b9b3361

ОТВЕТЫ

Ответ 1

Во-первых, M должен быть двумерным массивом, например double M[3][3]. Ваш массив, математически говоря, вектор 1x9, который не является обратимым.

  • N является указателем на int для порядок матрицы - в этом случае, N = 3.

  • A является указателем на LU факторизации матрицы, которая вы можете получить, запустив LAPACK подпрограмма dgetrf.

  • LDA является целым числом для "ведущего элемент "матрицы, что позволяет вы выбираете подмножество большего матрицы, если вы хотите просто инвертировать немного шт. Если вы хотите инвертировать вся матрица, LDA должна быть просто равным N.

  • IPIV - это сводные индексы матрица, другими словами, это список инструкций о том, какие строки для обмена чтобы инвертировать матрицу. IPIV должен быть создан LAPACK подпрограмма dgetrf.

  • LWORK и WORK - это "рабочие пространства", используется LAPACK. Если вы инвертируете вся матрица, LWORK должна быть int равно N ^ 2, а WORK - двойной массив с элементами LWORK.

  • INFO - это просто переменная статуса для скажите, работает ли операция успешно завершена. Поскольку не все матрицы обратимы, я бы рекомендуем отправить это кому-то вид системы проверки ошибок. INFO = 0 для успешной работы, INFO = -i, если i-й аргумент имел неправильное входное значение, а INFO > 0, если матрица не обратима.

Итак, для вашего кода я бы сделал что-то вроде этого:

int main(){

    double M[3][3] = { {1 , 2 , 3},
                       {4 , 5 , 6},
                       {7 , 8 , 9}}
    double pivotArray[3]; //since our matrix has three rows
    int errorHandler;
    double lapackWorkspace[9];

    // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix 
    // called A, sending the pivot indices to IPIV, and spitting error 
    // information to INFO.
    // also don't forget (like I did) that when you pass a two-dimensional array
    // to a function you need to specify the number of "rows"
    dgetrf_(3,3,M[3][],3,pivotArray[3],&errorHandler);
    //some sort of error check

    dgetri_(3,M[3][],3,pivotArray[3],9,lapackWorkspace,&errorHandler);
    //another error check

    }

Ответ 2

Вот рабочий код для вычисления обратного к матрице с использованием lapack в C/С++:

#include <cstdio>

extern "C" {
    // LU decomoposition of a general matrix
    void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);

    // generate inverse of a matrix given its LU decomposition
    void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);
}

void inverse(double* A, int N)
{
    int *IPIV = new int[N+1];
    int LWORK = N*N;
    double *WORK = new double[LWORK];
    int INFO;

    dgetrf_(&N,&N,A,&N,IPIV,&INFO);
    dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);

    delete IPIV;
    delete WORK;
}

int main(){

    double A [2*2] = {
        1,2,
        3,4
    };

    inverse(A, 2);

    printf("%f %f\n", A[0], A[1]);
    printf("%f %f\n", A[2], A[3]);

    return 0;
}

Ответ 3

Вот рабочая версия выше, использующая интерфейс OpenBlas для LAPACKE. Ссылка с библиотекой openblas (LAPACKE уже содержится)

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

// inplace inverse n x n matrix A.
// matrix A is Column Major (i.e. firts line, second line ... *not* C[][] order)
// returns:
//   ret = 0 on success
//   ret < 0 illegal argument value
//   ret > 0 singular matrix

lapack_int matInv(double *A, unsigned n)
{
    int ipiv[n+1];
    lapack_int ret;

    ret =  LAPACKE_dgetrf(LAPACK_COL_MAJOR,
                          n,
                          n,
                          A,
                          n,
                          ipiv);

    if (ret !=0)
        return ret;


    ret = LAPACKE_dgetri(LAPACK_COL_MAJOR,
                       n,
                       A,
                       n,
                       ipiv);
    return ret;
}

int main()
{
    double A[] = {
        0.378589,   0.971711,   0.016087,   0.037668,   0.312398,
        0.756377,   0.345708,   0.922947,   0.846671,   0.856103,
        0.732510,   0.108942,   0.476969,   0.398254,   0.507045,
        0.162608,   0.227770,   0.533074,   0.807075,   0.180335,
        0.517006,   0.315992,   0.914848,   0.460825,   0.731980
    };

    for (int i=0; i<25; i++) {
        if ((i%5) == 0) putchar('\n');
        printf("%+12.8f ",A[i]);
    }
    putchar('\n');

    matInv(A,5);

    for (int i=0; i<25; i++) {
        if ((i%5) == 0) putchar('\n');
        printf("%+12.8f ",A[i]);
    }
    putchar('\n');
}

Пример:

% g++ -I [OpenBlas path]/include/ example.cpp [OpenBlas path]/lib/libopenblas.a
% a.out

+0.37858900  +0.97171100  +0.01608700  +0.03766800  +0.31239800 
+0.75637700  +0.34570800  +0.92294700  +0.84667100  +0.85610300 
+0.73251000  +0.10894200  +0.47696900  +0.39825400  +0.50704500 
+0.16260800  +0.22777000  +0.53307400  +0.80707500  +0.18033500 
+0.51700600  +0.31599200  +0.91484800  +0.46082500  +0.73198000 

+0.24335255  -2.67946180  +3.57538817  +0.83711880  +0.34704217 
+1.02790497  -1.05086895  -0.07468137  +0.71041070  +0.66708313 
-0.21087237  -4.47765165  +1.73958308  +1.73999641  +3.69324020 
-0.14100897  +2.34977565  -0.93725915  +0.47383541  -2.15554470 
-0.26329660  +6.46315378  -4.07721533  -3.37094863  -2.42580445

Ответ 4

Вот рабочая версия примера Спенсера Нельсона выше. Одна загадка в том, что входная матрица находится в строчном порядке, хотя она, как представляется, вызывает базовую процедуру fortran dgetri. Я убежден, что все базовые подпрограммы fortran требуют порядка столбцов, но я не эксперт по LAPACK, на самом деле, я использую этот пример, чтобы помочь мне изучить его. Но эта одна тайна в стороне:

Входная матрица в примере является особой. LAPACK пытается сказать вам, что вернув 3 в errorHandler. Я изменил 9 в этой матрице на 19, получив сообщение errorHandler of 0, и сравнил результат с результатом Mathematica. Сравнение было также успешным и подтвердило, что матрица в примере должна быть в строчном порядке, как представлено.

Вот рабочий код:

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

int main() {
    int N = 3;
    int NN = 9;
    double M[3][3] = { {1 , 2 , 3},
                       {4 , 5 , 6},
                       {7 , 8 , 9} };
    int pivotArray[3]; //since our matrix has three rows
    int errorHandler;
    double lapackWorkspace[9];

    // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix
    // called A, sending the pivot indices to IPIV, and spitting error information
    // to INFO. also don't forget (like I did) that when you pass a two-dimensional
    // array to a function you need to specify the number of "rows"
    dgetrf_(&N, &N, M[0], &N, pivotArray, &errorHandler);
    printf ("dgetrf eh, %d, should be zero\n", errorHandler);

    dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler);
    printf ("dgetri eh, %d, should be zero\n", errorHandler);

    for (size_t row = 0; row < N; ++row)
    {   for (size_t col = 0; col < N; ++col)
        {   printf ("%g", M[row][col]);
            if (N-1 != col)
            {   printf (", ");   }   }
        if (N-1 != row)
        {   printf ("\n");   }   }

    return 0;   }

Я построил и запустил его следующим образом на Mac:

gcc main.c -llapacke -llapack
./a.out

Я сделал nm в библиотеке LAPACKE и нашел следующее:

liblapacke.a(lapacke_dgetri.o):
                 U _LAPACKE_dge_nancheck
0000000000000000 T _LAPACKE_dgetri
                 U _LAPACKE_dgetri_work
                 U _LAPACKE_xerbla
                 U _free
                 U _malloc

liblapacke.a(lapacke_dgetri_work.o):
                 U _LAPACKE_dge_trans
0000000000000000 T _LAPACKE_dgetri_work
                 U _LAPACKE_xerbla
                 U _dgetri_
                 U _free
                 U _malloc

и похоже, что существует оболочка LAPACKE [sic], которая, по-видимому, освобождает нас от необходимости принимать адреса везде для удобства fortran, но я, вероятно, не собираюсь его пытаться, потому что у меня есть путь вперед.

ИЗМЕНИТЬ

Вот рабочая версия, которая обходит LAPACKE [sic], напрямую используя LAPACK fortran. Я не понимаю, почему входной ряд строк дает правильные результаты, но я снова подтвердил его в Mathematica.

#include <stdio.h>
#include <stddef.h>

int main() {
    int N = 3;
    int NN = 9;
    double M[3][3] = { {1 , 2 ,  3},
                       {4 , 5 ,  6},
                       {7 , 8 , 19} };
    int pivotArray[3]; //since our matrix has three rows
    int errorHandler;
    double lapackWorkspace[9];
    /* from http://www.netlib.no/netlib/lapack/double/dgetrf.f
      SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
      *
      *  -- LAPACK routine (version 3.1) --
      *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
      *     November 2006
      *
      *     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
      *     ..
      *     .. Array Arguments ..
      INTEGER            IPIV( * )
      DOUBLE PRECISION   A( LDA, * )
    */

    extern void dgetrf_ (int * m, int * n, double * A, int * LDA, int * IPIV,
                         int * INFO);

    /* from http://www.netlib.no/netlib/lapack/double/dgetri.f
       SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
       *
       *  -- LAPACK routine (version 3.1) --
       *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
       *     November 2006
       *
       *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, LWORK, N
       *     ..
       *     .. Array Arguments ..
       INTEGER            IPIV( * )
       DOUBLE PRECISION   A( LDA, * ), WORK( * )
    */

    extern void dgetri_ (int * n, double * A, int * LDA, int * IPIV,
                         double * WORK, int * LWORK, int * INFO);

    // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix
    // called A, sending the pivot indices to IPIV, and spitting error information
    // to INFO. also don't forget (like I did) that when you pass a two-dimensional
    // array to a function you need to specify the number of "rows"
    dgetrf_(&N, &N, M[0], &N, pivotArray, &errorHandler);
    printf ("dgetrf eh, %d, should be zero\n", errorHandler);

    dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler);
    printf ("dgetri eh, %d, should be zero\n", errorHandler);

    for (size_t row = 0; row < N; ++row)
    {   for (size_t col = 0; col < N; ++col)
        {   printf ("%g", M[row][col]);
            if (N-1 != col)
            {   printf (", ");   }   }
        if (N-1 != row)
        {   printf ("\n");   }   }
    return 0;   }

построено и работает следующим образом:

$ gcc foo.c -llapack
$ ./a.out
dgetrf eh, 0, should be zero
dgetri eh, 0, should be zero
-1.56667, 0.466667, 0.1
1.13333, 0.0666667, -0.2
0.1, -0.2, 0.1

ИЗМЕНИТЬ

Тайна больше не кажется тайной. Я думаю, что вычисления выполняются в порядке столбцов, как и должно быть, но я ввожу и печатаю матрицы так, как если бы они были в порядке строк. У меня есть две ошибки, которые отменяют друг друга, поэтому вещи выглядят как строки, даже если они столбцы.