matlab 对称带状矩阵向量乘法

nwlls2ji  于 2022-11-15  发布在  Matlab
关注(0)|答案(1)|浏览(208)

尝试将对称带状矩阵乘以向量,即简单乘法
Ax=b
其中A是对称频带矩阵,以here的格式存储。因此,它的较低部分不会保存到存储中。
cols是包含非零带折射率的矢量,例如cols(1,1)=1是主对角线带,cols(1,2)=2是位于主对角线+1位移处的上带,
cols的第二行是带的长度,这个值是常量,我最初存储它是为了获得矩阵大小。
代码的第一部分
1.读取文件,将表转换为二维数组
1.生成随机向量x进行结果检验
1.将对称频带形式转换为正则矩阵,目的与上述相同
所以我们的想法是取时间向量x,然后做元素乘积。比方说,将矩阵的主对角线的每个元素与相应的向量元素相乘,然后递增RHS向量b
对于移动的波段,我对x向量进行切片,并递增相应的b
我设法将上三角形与x向量相乘,max(abs(triu(asps)*x-b))的差值为零
但当我尝试做下矩阵部分max(abs(asps
x-b))的乘法时,给出了一些错误。
会有什么问题呢?也许索引是错误的?
在高位乘法后或将x数组切片为xTimesLower时,包含b(i+width)的行应出现错误
Band Matrix file

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)));
bmvo0sr5

bmvo0sr51#

看起来代码正在工作,只是在整数矩阵上测试,差异为零,
由于浮点加法的非结合性,类型Double的结果相差10^-14到10^-15
在整型代码的情况下

clc; clear; close all;
adnsT = readtable('denseCTAC.txt','ReadVariableNames', false);
cols = [1 2 3 9 10 11 19;...
    81 80 79 73 72 71 63];
adns = floor(table2array(adnsT)*100);
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 = randi(10,n);
x = x(:,1);
b = x*0;
% x = x';

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
            b(i+width) = b(i+width) + adns(i,j)*xTimesLower(i);
            %%%% multiplying lower part
%             b(i+width) = b(i+width) + adns(i,j)*xTimesLower(i); % shift issue?
        end
    end
end
% max(triu(asps)*x-b)
max(abs(asps*x-b))

% fprintf("Error=%d\n",max(abs(asps*x-b)));

相关问题