1

我有一个关于在 Octave 上运行 SOR 算法的学校项目，但我的效率非常低。所以我有这段代码：

for ii=1:n
    r = 1/A(ii,ii);
    for jj=1:n
        if (ii!=jj)
            A(ii,jj) = A(ii,jj)*r;
        end;
    end;
    b(ii,1) = b(ii,1)*r;
    x(ii,1) = b(ii,1);
end;

我怎样才能矢量化这个？我的第一次尝试是这样的：

for ii=1:n
    r = 1/A(ii,ii);
    A(find(eye(length(A))!=1)) = A(find(eye(length(A))!=1))*r;
    b(ii,1) = b(ii,1)*r;
    x(ii,1) = b(ii,1);
end;

但我不确定它有多大帮助。有没有更好和/或更有效的方法来做到这一点？

谢谢！

4

3 回答 3

4

我相信你可以完全避免循环。您必须看到您将 的对角线元素视为 的倒数A，那么您为什么要使用循环。直接做。这是第一步。删除Inf，现在您想乘以r相应行的相应非对角元素，对吗？

因此，请使用repmat构造这样的矩阵，该矩阵将沿列复制元素，因为您将相同的乘以r(1,2)、(1,3)、...、(1,n)。但R具有非零对角元素。因此使它们为零。现在你会得到你的 A 除了对角线元素将为零。因此，您只需将它们从 original 添加回来A。这可以通过A=A.*R+A.*eye(size(A,1)).

矢量化来自经验，最重要的是分析您的代码。在每个步骤中考虑是否要使用循环，如果不使用等效命令替换该步骤，则将遵循其他代码（例如，我构造了一个矩阵R，而您正在构造单个元素r。所以我只考虑转换r->R然后其余的代码将就位）。

代码如下：

R=1./(A.*eye(size(A,1))); %assuming matrix A is square and it does not contain 0 on the main diagonal
R=R(~isinf(R));
R=R(:);
R1=R;
R=repmat(R,[1,size(A,2)]);
R=R.*(true(size(A,1))-eye(size(A,1)));
A=A.*R+A.*eye(size(A,1)); %code verified till here since A comes out to be the same

b = b.*R1;
x=b;
于 2013-04-16T00:35:44.083 回答
2

我想有矩阵：

A (NxN)
b (Nx1)

编码：

d = diag(A);
A = diag(1 ./ d) * A + diag(d - 1);
b = b ./ d;
x = b;
于 2013-04-16T17:58:27.440 回答
2

随机碰到这个问题，乍一看看起来很有趣，因为这个问题被标记为vectorization问题。

我能够提出一个bsxfun基于矢量化的解决方案，该解决方案也使用diagonal indexing. 这个解决方案似乎让我3-4x加快了具有适当大小输入的循环代码。

假设你仍然有兴趣看到这个问题的加速改进，我很想知道你会得到什么样的加速。这是代码 -

diag_ind = 1:size(A,1)+1:numel(A);
diag_A = A(diag_ind(:));
A = bsxfun(@rdivide,A,diag_A);
A(diag_ind) = diag_A;
b(:,1) = b(:,1)./diag_A;
x(:,1) = b(:,1);

让我知道！

于 2014-09-30T19:47:04.753 回答