用于 GMRES 的大矩阵的 C 或 C++ 矩阵矢量积的更快方法

我有一个大而密集的矩阵A，我的目标是使用迭代方法找到线性系统Ax=b的解（在MATLAB中，计划使用其内置的GMRES）。对于超过 10,000 行，这对于我的计算机来说太多了，无法存储在内存中，但我知道 A 中的条目是由长度为 N 的两个已知向量 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 命令接受可以计算矩阵向量积 A*x 的函数调用作为输入，这允许我处理比我可以存储在内存中的更大的矩阵。为了编写矩阵-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 反转我感兴趣的这个矩阵，这需要一种更快的方法来执行矩阵向量积，而无需显式将数组加载到内存中。谢谢！