Trying to multiply symmetric band matrix times vector, i.e. simple multiplication

A*x = b

where A is symmetric band matrix, stored in format as here. So the lower part of it is not saved into the storage.

cols is a vector containing the indices of nonzero bands, e.g. cols(1,1)=1 is main diagonal band, cols(1,2)=2 is upper band located at 1 shift from main diagonal,

the second row of cols is length of the band, this value is constant, I initially stored it to get the matrix size.

The first part of the code

1. Reads the file, converts table into 2d array

2. Generates random vector x to test the result

3. Converts symmetric band form into regular matrix for the same purpose as above

So the idea is to take the times vector x and do the element wise product. Say, multiply each element of main diagonal of a matrix with corresponding vector element, then increment the rhs vector b.

For shifted bands I do slicing of the x vector and increment corresponding b

I managed to multiply the upper triangle with the x vector, the difference in max(abs(triu(asps)*x-b)) is zero

but when I try to do the multiplication of lower matrix part max(abs(asps*x-b)) gives some error.

What could be an issue? Maybe indexing is wrong?

Error should be somewhere at the line containing b(i width) after upper multiplication or when I slice the x array into xTimesLower

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

CodePudding user response：

it looks like the code is working just tested on integer matrices and difference is zero,

due to the non associativity of floating point addition, the results for type double differ by 10^-14 to 10^-15

below the integer code case

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