如何使用加速框架执行矩阵逆运算？

Posted 2023-02-25

【中文标题】如何使用加速框架执行矩阵逆运算？

【英文标题】：How to perform matrix inverse operation using the accelerate framework?

【发布时间】：2012-07-01 14:25:30

【问题描述】：

我想求矩阵的逆。

我知道这涉及到第一个 LU 因式分解，然后是求逆步骤，但我无法通过搜索苹果的 10.7 文档找到所需的函数！

这似乎是一篇有用的帖子Symmetric Matrix Inversion in C using CBLAS/LAPACK，指出应该使用sgetrf_ 和sgetri_ 函数。但是搜索这些术语我在 Xcode 文档中一无所获。

有人有这个矩阵运算的样板代码吗？

【参考方案1】：

Apple 根本没有记录 LAPACK 代码，我猜是因为他们只是实现了来自netlib.org 的标准接口。遗憾的是，您无法从内置 Xcode 文档中搜索这些函数名称，但是解决方案相当简单：只需在 URL 中指定函数名称，例如dgetrf_() 转至 http://www.netlib.org/clapack/what/double/dgetrf.c。

要对矩阵求逆，需要两个 LAPACK 函数：dgetrf_()，它执行 LU 分解，dgetri_()，它获取前一个函数的输出并进行实际的求逆。

我使用 Xcode 创建了一个标准的应用程序项目，添加了 Accelerate 框架，创建了两个 C 文件：matinv.h、matinv.c 并编辑了 main.m 文件以删除 Cocoa 内容：

// main.m

#import "matinv.h"

int main(int argc, char *argv[])

    int N = 3;
    double A[N*N];
    A[0] = 1; A[1] = 1; A[2] = 7;
    A[3] = 1; A[4] = 2; A[5] = 1;
    A[6] = 1; A[7] = 1; A[8] = 3;
    matrix_invert(N, A);
    //        [ -1.25  -1.0  3.25 ]
    // A^-1 = [  0.5       1.0  -1.5  ]
    //        [  0.25   0.0 -0.25 ] 
    return 0;

现在是头文件，

//  matinv.h

int matrix_invert(int N, double *matrix);

然后是源文件，

int matrix_invert(int N, double *matrix) 

    int error=0;
    int *pivot = malloc(N*sizeof(int)); // LAPACK requires MIN(M,N), here M==N, so N will do fine.
    double *workspace = malloc(N*sizeof(double));

    /*  LU factorisation */
    dgetrf_(&N, &N, matrix, &N, pivot, &error);

    if (error != 0) 
        NSLog(@"Error 1");
        free(pivot);
        free(workspace);
        return error;
    

    /*  matrix inversion */
    dgetri_(&N, matrix, &N, pivot, workspace, &N, &error);

    if (error != 0) 
        NSLog(@"Error 2");
        free(pivot);
        free(workspace);
        return error;
    

    free(pivot);
    free(workspace);
    return error;

【讨论】：

规范的 LAPACK 参考是 LAPACK 用户指南。 (netlib.org/lapack/lug)

我在扫描这个神秘的（至少可以说）LAPACK 库时遇到问题。如何将此代码调整为单精度浮点数？

哦，我找到了：sgetrf_ 和 sgetri_（S 代表“单精度”？）

感谢您的代码，它帮助了我。你的matrix_invert有bug，如果有错误，跳过pivot和workspace的释放，提前返回。

我认为枢轴应该是来自文档“维度 (min(M,N))”的尺寸 N*sizeof(int)