function [r] = sbmvmul(M,B,v)
% [r] = sbmvmul(M,B,v)
%
% Complex-symmetric and block tridiagonal matrix-vector multiplication.
% Inner product version.
%
% Input
%   M   M(:,:,i) are the main diagonal blocks
%   B   B(:,:,i) are the subdiagonal blocks
%   v   column vector
% Outputs
%   r   product vector so that
%         r = J*v,
% where J is complex-symmetric and block tridiagonal whose main diagonal
% blocks are M(:,:,i) and subdiagonal blocks are B(:,:,i).
%
[b,b,k] = size(M);              % get block size and number of diagonal blocks
n = b*k;                        % matrix size
%
r = zeros(n,1);                 % initialize a column vector
%
r(1:b) = M(:,:,1)*v(1:b) + B(:,:,1).'*v(b+1:2*b);
%
for i = 1:k-2
    low = i*b +1;               % lower bound
    up = (i+1)*b;               % upper bound
    r(low:up) = B(:,:,i)*v((low-b):(low-1)) + M(:,:,i+1)*v(low:up) + ...
                B(:,:,i+1).'*v(up+1:up+b);
end
%
r(n-b+1:n) = B(:,:,k-1)*v((n-2*b+1):(n-b)) + M(:,:,k)*v(n-b+1:n);