【发布时间】：2014-08-05 02:46:04
【问题描述】：
背景： 我正在将特定应用程序的特定 MATLAB 函数转换为 C，以通过 MEX 减少总体内存开销并提高特定算法的性能。包括的各种功能中的一些是，例如，分解/线性系统求解器算法及其在类似于 MATLAB 左除运算符的层次结构中的实现。
我意识到 MATLAB 使用 JIT 编译器，许多人报告说很难超越 MATLAB 的性能等等，但是这里的目标是实际测试这个理论并真正找出性能或效率的任何提高是否可以达到了。
我目前正在使用 GNU 的 GSL 库。
解决问题...
对于基本的对角矩阵求逆问题，我没有得到正确的结果。
在 MATLAB 中，代码：
x = A\B
其中 A 是一个稀疏的对角矩阵，其中对角线上的每个元素都非零，而 B 只是一个正则填充矩阵，满足矩阵与 A 相乘的大小要求。
为了解决这个问题，不应该只是一个基本的 for 循环，用该元素的倒数或 1/元素替换主对角线中的每个元素吗？
然后当然是乘以矩阵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 的数字代码。我的精度设置为正确输出。
有什么想法吗？
【问题讨论】：
-
要检查的几件事。 1）在Matlab中使用相同的算法，即做
x=AA*B，其中AA=1./A。我很确定x=A/B不会做inv(A)*B，而是在不计算逆矩阵的情况下解决系统A*x=B。由于舍入错误导致执行的操作不相同，因此您可能会得到一些不同的值。 2) 为了排除 GSL 库代码的不正确使用，自己使用循环进行乘法（即 x(i,k)=sum(1/A(i,j)*B(j,k)，其中 sum 超过 j） . 3) 写入文件时，这些值可能变为零。检查您如何写入/读取文件。