c - 基本线性代数算法(MATLAB 和 C)

Question

背景: 我正在将特定应用程序的特定 MATLAB 函数转换为 C，以减少总体内存开销并通过 MEX 提高特定算法的性能。包括的各种函数中的一些是，例如，分解/线性系统求解器算法及其在类似于 MATLAB 的左除运算符的层次结构中的实现。

我意识到 MATLAB 使用 JIT 编译器，许多人报告说很难超越 MATLAB 的性能等等，但是这里的目标是实际测试这个理论并真正找出是否可以提高性能或效率实现了。

我目前正在使用 GNU 的 GSL 库。

关于问题...

我没有得到基本对角矩阵求逆问题的正确结果。

在 MATLAB 中，代码:

x = A\B

其中 A 是一个稀疏的对角矩阵，其中对角线上的每个元素都非零，而 B 只是一个规则的填充矩阵，满足与 A 的矩阵乘法的大小要求。

要解决这个问题，难道不应该只是一个基本的 for 循环，将主对角线中的每个元素替换为该元素的倒数，或者 1/element？

然后当然是乘以矩阵 B。这似乎是一个简单的问题，但它不会产生与 MATLAB 相同的结果，即使输入矩阵相同。

C 代码:

    gsl_matrix* matrix1;  
    gsl_matrix* matrix1_copy;
    int index; 

    ...  //  A gets initialized and filled as a sparse diagonal matrix

    gsl_matrix_memcpy(matrix1_copy, matrix1);
    gsl_vector diag = gsl_matrix_diagonal(matrix1_copy).vector;
    for (index = 0; index < diag.size; index++) {
        gsl_matrix_set(matrix1_copy, index, index, 1/gsl_vector_get(&diag, index));
    }

这被证实用它们的逆元素替换了所有对角线元素。返回此矩阵，然后用作此函数中的第一个矩阵:

gsl_matrix* multiply(gsl_matrix* matrix1, gsl_matrix* matrix2) {
gsl_matrix* final_matrix = gsl_matrix_calloc(matrix1->size1, matrix2->size2);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, matrix1, matrix2, 0.0, final_matrix);
return final_matrix;
}

返回的矩阵被写入文件，因此我可以轻松读取它。整个过程既没有错误也没有警告，每一步似乎都完成了，但最终的矩阵与 MATLAB 中的对应矩阵不匹配；一开始很难说，因为这是一个非常大的稀疏矩阵，但是两个被比较的矩阵的开头是正确的，但是在 MATLAB 中输入矩阵的右下角后，我观察到 C 中的数字只是 0代码。我的精度设置为正确输出。

有什么想法吗？

Answer 1

您应该使用回代并将对角矩阵 A 视为特例。这样，您将使用相同的算法(和相同的单元测试!)涵盖更广泛的问题。

第一个额外的评论:由于您是从头开始开发算法，因此您必须使用简单的测试用例。用复杂的案例进行测试将被证明是困难的。我的意思是，您的测试用例必须非常小，以至于您可以手动验证结果。

第二个附加评论:Matlab 使用由 Intel aka 优化的 BLAS/LAPACK 库 Intel Math Kernel Library .如果您的目标是击败 Matlab，您应该使用 BLAS 或类似 BLAS 的库:goto BLAS , OpenBLAS , 或 FLAME .

最后的补充评论:你不是我认识的第一个想打败 Matlab 的人。从别人的错误中学习。选择一个目标平台并进行优化。深入了解您的处理器。使用多线程。使用 MPI。以此为圣经:Agner's optimization techniques