使用C中的lapack计算矩阵的逆矩阵

Question

我希望能够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 的逆？

Answer 1

这是使用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+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;
}

Answer 2

首先,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,如果第i个参数具有不正确的输入值,则INFO = -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

    }