23

我希望能够NxN使用 lapack 在 C/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;
}

您将如何完成它以3x3使用 dgetri_ 获得矩阵 M 的逆?

4

4 回答 4

26

这是在 C/C++ 中使用 lapack 计算矩阵逆的工作代码:

#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];
    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;
}
于 2010-08-19T18:48:12.273 回答
18

首先,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 应该是一个等于 N^2 的 int,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

    }
于 2010-08-19T08:49:53.760 回答
2

这是上面使用 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 
于 2017-01-04T23:16:20.540 回答
1

这是上面 Spencer Nelson 示例的工作版本。关于它的一个谜团是输入矩阵按行优先顺序排列,即使它似乎调用了底层的 fortran 例程dgetri。我被引导相信所有底层的 fortran 例程都需要列优先顺序,但我不是 LAPACK 专家,事实上,我正在使用这个示例来帮助我学习它。但是,除了一个谜:

示例中的输入矩阵是奇异的。LAPACK 试图通过3errorHandler. 我9将该矩阵中的 更改为 a 19,获得成功errorHandler0信号,并将结果与​​来自 的结果进行比较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 [原文如此] 包装器,大概可以让我们不必为了 fortran 的方便而到处获取地址,但我可能不会去尝试它,因为我有前进的道路。

编辑

这是一个绕过 LAPACKE [原文如此] 的工作版本,直接使用 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

编辑

这个谜似乎不再是一个谜。我认为计算是按列优先顺序完成的,因为它们必须,但我输入和打印矩阵,就好像它们是按行优先顺序一样。我有两个相互抵消的错误,所以即使它们是列式的,事情看起来也像是行式的。

于 2016-05-20T15:56:55.103 回答