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

5
我想找到矩阵的逆。
我知道这需要先进行LU分解,然后进行求逆步骤,但是我在搜索苹果10.7文档时找不到所需的函数!
这篇文章似乎很有用:使用CBLAS/LAPACK在C中进行对称矩阵求逆,指出应该使用sgetrf_sgetri_函数。但是在搜索这些术语时,我在Xcode文档中找不到任何内容。
有没有人有这个矩阵操作的样板代码?
1个回答

15

苹果公司没有对LAPACK代码进行文档记录,我猜测这是因为他们只是实现了来自netlib.org的标准接口。很遗憾你不能在内置的Xcode文档中搜索这些函数名称,但解决方案非常简单:只需在URL中指定函数名称,例如对于dgetrf_(),转到http://www.netlib.org/clapack/what/double/dgetrf.c

要反转一个矩阵需要使用两个LAPACK函数:dgetrf_(),它执行LU分解,以及dgetri_(),它接收前一个函数的输出并执行实际的反转。

我使用Xcode创建了一个标准应用程序项目,添加了加速框架,创建了两个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;
}

1
LAPACK 的规范参考是 LAPACK 用户指南。(http://www.netlib.org/lapack/lug/) - Stephen Canon
我在扫描LAPACK库时遇到了一些困难,因为它非常晦涩难懂。如何将这段代码调整为单精度浮点数? - Nicolas Miari
哦,我找到了:sgetrf_和sgetri_(S代表“单精度”?) - Nicolas Miari
1
谢谢你的代码,它帮了我很大的忙。在你的matrix_invert函数中有一个bug,如果出现错误,你会跳过释放pivot和workspace的步骤并提前返回。 - X.Y.
我认为根据文档中的“维数(最小值为(M,N))”,枢轴应该是大小为N*sizeof(int) - combinatorial
真的,我分配了太多的枢轴空间!已修复。 - Daniel Farrell

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接