基于压缩稀疏行(CSR)格式的Jacobi迭代法实现问题求助

我尝试基于压缩稀疏行（CSR）格式实现Jacobi迭代法，但无法得到正确输出。使用的是4×4三对角矩阵，期望输出为 [28.1987 47.3978 48.5979 25.7986]，以下是我的MATLAB代码：

clear all;
close all;
clc;

H=4;
a=2;
b=-1;
c=-1;
A = diag(a*ones(1,H)) + diag(b*ones(1,H-1),1) + diag(c*ones(1,H-1),-1);%Matrix A
n = size(A,1); % no of rows
m = size(A,2); % no of columns

V = [];
C = [];
R = [];

counter=1;
R= [counter];
for i=1:n
  for j=1:m
    if (A(i,j) ~= 0)
      V = [V A(i,j)];
      C = [C j];
      counter=counter+1;
    end
    R(i+1)=counter;
  end
end

b = [9,18,24,3];
x_new = [1 ; 1 ; 1 ; 1];

eps = 1e-5;  % 1 x 10^(-10).
error = 1000; % use any large value greater than eps to make sure that the loop can work

counter2=1;
while (error > eps)

  x_old = x_new;

  for i=1:length(R)-1  %modified
    t = 0;
    for j=R(i):R(i+1)-1 %modified
      if (C(j)~=i)   %not equal
   
        t = t + x_old(C(j))*A(i,C(j)); %modified
      end
    end
    x_new(i,1) = (b(i) - t)/A(i,C(j)); % is a row vector
  end
  error = norm(x_new-x_old);
  counter2=counter2+1;
end

x_new % print x

问题点分析 #

• CSR数组构造错误：原代码中R(i+1)=counter放在内层j循环的每次迭代中，导致每遍历一列就更新一次行指针，不符合CSR格式定义。CSR的R数组应记录每行非零元素的起始索引，正确做法是处理完第i行所有列后，再将counter赋值给R(i+1)。
• 未利用CSR存储的非零值：计算t时直接引用原矩阵A(i,C(j))，既浪费CSR的优势，也容易出错，应该使用CSR的V(j)来获取非零元素值。
• 分母引用错误：计算x_new(i)时，分母A(i,C(j))中的j是内层循环结束后的最后一个索引，并非第i行的对角元位置。正确的分母应为第i行的对角元，可以直接取A(i,i)，或者在CSR中定位对角元对应的V值。
• 向量维度不统一：b定义为行向量，与x_new的列向量维度不一致，虽然MATLAB会自动兼容，但统一维度更规范。

修正后的代码 #

clear all;
close all;
clc;

H=4;
a=2;
b=-1;
c=-1;
A = diag(a*ones(1,H)) + diag(b*ones(1,H-1),1) + diag(c*ones(1,H-1),-1);%Matrix A
n = size(A,1); % no of rows

% 正确构造CSR格式
V = [];
C = [];
R = [1]; % 行指针起始为1
counter = 1;
for i=1:n
    for j=1:n
        if A(i,j) ~= 0
            V = [V, A(i,j)];
            C = [C, j];
            counter = counter + 1;
        end
    end
    R = [R, counter]; % 处理完一行后更新行指针
end

% 统一b为列向量
b = [9; 18; 24; 3];
x_new = ones(n, 1); % 初始解向量

eps = 1e-5;
error = 1000;
counter2 = 1;

while error > eps
    x_old = x_new;
    
    for i=1:n
        t = 0;
        diag_val = A(i,i); % 直接取对角元，或者从CSR中查找
        % 遍历第i行的所有非零元素
        for j=R(i):R(i+1)-1
            if C(j) ~= i
                t = t + x_old(C(j)) * V(j); % 使用CSR的V数组
            else
                diag_val = V(j); % 也可以从CSR中获取对角元
            end
        end
        x_new(i) = (b(i) - t) / diag_val;
    end
    
    error = norm(x_new - x_old);
    counter2 = counter2 + 1;
end

disp('迭代结果：');
disp(x_new);

修正说明 #

1. 调整CSR的R数组构造逻辑，确保每行处理完毕后再更新行指针，符合CSR格式规范。
2. 计算t时使用CSR存储的V(j)替代原矩阵A，充分利用稀疏存储的优势。
3. 从CSR中获取对角元值作为分母，避免引用错误的索引。
4. 将b改为列向量，统一向量维度。

运行修正后的代码，即可得到期望的输出结果。

内容的提问来源于stack exchange，提问作者rishrish