c++ - 在 GMRES 中使用 C 或 C++ 中大型矩阵的矩阵 vector 乘积的更快方法

Question

我有一个大而密集的矩阵 A，我的目标是使用迭代方法找到线性系统 Ax=b 的解(在 MATLAB 中是使用其内置 GMRES 的计划)。对于超过 10,000 行，这对我的计算机来说太多了，无法在内存中存储，但我知道 A 中的条目是由两个已知的长度为 N 的 vector x 和 y 构造的，并且条目满足:A(i,j) = .5*(x[i]-x[j])^2+([y[i]-y[j])^2 * log(x[i]-x[j] )^2+([y[i]-y[j]^2).

MATLAB 的 GMRES 命令接受一个可以计算矩阵 vector 乘积 A*x 的函数调用作为输入，这使我能够处理比我可以存储在内存中的更大的矩阵。为了编写 matrix-vecotr 乘积函数，我首先在 matlab 中逐行尝试并使用一些矢量化，但我避免生成整个数组 A(因为它太大)。不幸的是，这在我申请 GMRES 时相当缓慢。我的计划是为 MATLAB 编写一个 mex 文件，它是用 C 编写的，理想情况下应该比 matlab 代码快得多。我是 C 的新手，所以这很糟糕，我天真的尝试用 C 编写代码比我在 Matlab 中部分矢量化的尝试慢。

#include <math.h>
#include "mex.h"
void Aproduct(double *x, double *ctrs_x, double *ctrs_y, double *b, mwSize n)
{
    mwSize i;
    mwSize j;
    double val;
    for (i=0; i<n; i++) {
        for (j=0; j<i; j++) {
            val = pow(ctrs_x[i]-ctrs_x[j],2)+pow(ctrs_y[i]-ctrs_y[j],2);

            b[i] = b[i] + .5* val * log(val) * x[j];
        }
        for (j=i+1; j<n; j++) {
            val = pow(ctrs_x[i]-ctrs_x[j],2)+pow(ctrs_y[i]-ctrs_y[j],2);

            b[i] = b[i] + .5* val * log(val) * x[j];
        }
    }
}

以上是 matlab mex 文件代码的计算部分(如果我理解正确的话，它是稍微修改过的 C)。请注意，我跳过了 i=j 的情况，因为在这种情况下变量 val 将是 0*log(0)，对我来说它应该被解释为 0，所以我直接跳过它。

有没有更有效或更快速的写法？当我在 matlab 中通过 mex 文件调用这个 C 函数时，它非常慢，甚至比我使用的 matlab 方法还慢。这让我感到惊讶，因为我怀疑 C 代码应该比 matlab 快得多。

我正在与之比较的部分矢量化的替代 matlab 方法是

function Ax = Aprod(x,ctrs)
n = length(x);
Ax = zeros(n,1);
for j=1:(n-3)
    v = .5*((ctrs(j,1)-ctrs(:,1)).^2+(ctrs(j,2)-ctrs(:,2)).^2).*log((ctrs(j,1)-ctrs(:,1)).^2+(ctrs(j,2)-ctrs(:,2)).^2);

    v(j)=0;
    Ax(j) = dot(v,x(1:n-3);
end

(n-3是因为实际上多了3个组件，但是是分开处理的，所以去掉了那段代码)。这是部分向量化的，只需要一个 for 循环，所以它更快是有道理的。但是，我希望我可以使用 C+mex 文件走得更快。

如有任何建议或帮助，我们将不胜感激!谢谢!

编辑:我应该更清楚。我愿意接受任何可以帮助我使用 GMRES 反转我感兴趣的矩阵的更快方法，这需要一种更快的方法来执行矩阵 vector 乘积，而无需将数组显式加载到内存中。谢谢!

Answer 1

如果你有Parallel Computing Toolbox和 MATLAB Distributed Computing Server ，您可以直接使用反斜杠求解大型密集线性系统。 (如果您没有可用的集群，您可能想使用 Amazon EC2 machines )。像这样:http://www.mathworks.co.uk/help/distcomp/examples/benchmarking-a-b.html