Related to this question Thanks to the idea of @veryreverie, I was eventually able to multiply the block diagonal matrix by a vector with some modifications. (Matrix structure is at the end of post). However something is wrong. When I run the algorithm in MATLAB, there no issues.
The difference in direct multiplication and my result of C'*x - bn is within 1e-15 range due to machine addition non-associativity. However, when trying the same thing in Fortran it differs significantly, up to NaN value:
f90 code:
module RHS
implicit none
contains
subroutine RHS_eval(m,n,C,MM,q,LAP,BCterm,ADV,bn)
real*8, intent (in) :: x(:),C(:,:)
real*8, intent (out) :: bn(:)
integer*8 :: row, col
integer*8, intent (in) :: m,n
integer*8 :: offset, blkSz, blkShft, i
integer*8 :: xstart, ystart
bn = 0
bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) + x(1:(m-1)*(n-1));
bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) - x(m:(m-1)*n);
blkSz = m;
offset = (m-1)*(n);
blkShft = m-1;
do i = 1,(n-1)
xstart = offset + blkSz*(i-1)+1;
ystart = (i-1)*(m-1)+1;
bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) - x(xstart:xstart+blkShft-1);
bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) + x(xstart+1:xstart+blkShft);
enddo
end subroutine RHS_eval
end module RHS
MATLAB code:
n = 15; m = 10;
x = (10*rand((m-1)*n+(n-1)*m,1));
bn = zeros((n-1)*(m-1),1);
bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) + x(1:(m-1)*(n-1));
bn(1:(m-1)*(n-1)) = bn(1:(m-1)*(n-1)) - x(m:(m-1)*n);
blkSz = m;
offset = (m-1)*(n);
blkShft = m-1;
for i = 1:(n-1)
xstart = offset + blkSz*(i-1)+1;
ystart = (i-1)*(m-1)+1;
bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) - x(xstart:xstart+blkShft-1);
bn(ystart:ystart+blkShft-1) = bn(ystart:ystart+blkShft-1) + x(xstart+1:xstart+blkShft);
end
We take the transpose of matrix C and multiply from RHS by a vector. There is no need to store the matrix since its structure is known.
Fixed a tiny bit, at early iterations it kinda works as expected, the number shown is the difference between matmul and my multiplication algo
maxval(abs(bn - matmul(transpose(C),(x))))
But starting at around 60th step starts to blow up
This issue only happens for the case when m and n are not equal