对称带状矩阵与向量乘法实现问题排查

我正在实现对称带状矩阵与向量的乘法（A*x = b），其中对称带状矩阵A按维基百科的带状矩阵格式存储，仅保存上三角部分。cols向量存储非零带的索引及长度：第一行是带的偏移（比如cols(1,1)=1对应主对角线），第二行是对应带的长度。

代码流程为：读取文件转为数组，生成随机测试向量x，将带状格式转为常规矩阵用于验证。我的实现思路是逐元素相乘：主对角线元素与x对应元素相乘后累加到b；偏移带则通过切片x来更新b。目前上三角部分的乘法验证正确（max(abs(triu(asps)*x-b))为0），但实现下三角部分乘法时，max(abs(asps*x-b))出现误差，怀疑是索引问题，可能出在b(i+width)或xTimesLower的切片处。

附实现代码：

clc; clear; close all;
adnsT = readtable('C:\Users\user\Documents\denseCTAC.txt','ReadVariableNames', false);
cols = [1 2 3 9 10 11 19;...
    81 80 79 73 72 71 63];
adns = table2array(adnsT);
n = 81;
for i = 1:7
   len = cols(2,i);
   colShift = cols(1,i)-1;
   for j = 1:len
       asps(j,j+colShift)=adns(j,i);
   end
end
asps = asps' + asps - eye(n).*diag(asps);
% spy(asps);

x = rand(n);
x = x(:,1);
b = x*0;

for j = 1:7
    width = cols(1,j)-1;
    for i = 1:n-width
        if width == 0
            b(i) = b(i)+adns(i,1)*x(i);
        else
            xTimesUpper = x(1+width:n);
            xTimesLower = x(1:n-width);
            
            %%%% multiplying upper part
            b(i) = b(i) + adns(i,j)*xTimesUpper(i); % OK if compared to triuA*x-b
            
            %%%% multiplying lower part
            % b(i+width) = b(i+width) + adns(i,j)*xTimesLower(i); % shift issue?
        end
    end
end
max((asps)*x-b)
fprintf("Error=%d\n",max(abs(asps*x-b)));

你的下三角部分索引逻辑存在问题，核心是未正确对应对称矩阵的元素映射关系，以下是修复方案：

1. 对称矩阵的下三角元素A[i+width][i]与上三角元素A[i][i+width]值相同，均存储在adns(i,j)中。
2. 修复循环逻辑：
   • 直接使用cols中给出的带长度len = cols(2,j)作为循环上限，避免计算n-width带来的潜在越界
   • 移除冗余的xTimesUpper和xTimesLower切片，直接通过索引访问x元素，减少混淆

修复后的乘法循环代码：

for j = 1:7
    width = cols(1,j)-1;
    len = cols(2,j);
    for i = 1:len
        if width == 0
            b(i) = b(i) + adns(i,j)*x(i);
        else
            % 上三角部分：A[i][i+width] * x[i+width] → 累加到b[i]
            b(i) = b(i) + adns(i,j)*x(i + width);
            % 下三角部分：A[i+width][i] * x[i] → 累加到b(i+width)
            b(i + width) = b(i + width) + adns(i,j)*x(i);
        end
    end
end

3. 验证：修复后运行代码，max(abs(asps*x-b))会趋近于0（浮点精度范围内），说明乘法逻辑正确。

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