% function C = TridTril(A,B)
% A: tridiagonal matrix
% B: Lower triangular matrix
% n = input('enter size of problem:');
disp(sprintf(' n f1 f2 f3 f4 t1 t2 t3 t4'))
disp(sprintf('------------------------------------------------------------------------------'))
for n = 8:8:80
v1 = rand(n-1,1);
v2 = rand(n,1);
v3 = rand(n-1,1);
A = diag(v1,-1)+diag(v2)+diag(v3,1);
% define nx3 array to hold the nontrivial diagonals of A as columns
A1 = zeros(n,3);
for i = 2:n-1
A1(i,1:3)=A(i,i-1:i+1);
end
A1(1,1)=0; A1(n,3)=0; A1(1,2:3)=A(1,1:2); A1(n,1:2)=A(n,n-1:n);
B=tril(rand(n,n));
flops(0);tic; C1 = A*B ; f1 = flops;
t1=toc;
flops(0); tic;C2 = sparse(A)*sparse(B) ; f2 = flops;
t2=toc;
C3=zeros(n,n);
flops(0);tic;
C3(1,:)= A1(1,2)*B(1,:)+A1(1,3)*B( 2,:);
for j = 2:n-1
C3(j,:)=A1(j,1)*B(j-1,:)+A1(j,2)*B(j,:)+A1(j,3)*B(j+1,:);
end
C3(n,:)=A1(n,1)*B(n-1,:)+A1(n,2)*B(n,:);
f3 = flops;t3=toc;
E3 = C3-C1;
error3 = norm(E3,inf);
C4=zeros(n,n);
flops(0);tic;
C4(1,1: 2)= A1(1,2)*B(1,1: 2)+A1(1,3)*B( 2,1: 2);
for j = 2:n-1
C4(j,1:j+1)=A1(j,1)*B(j-1,1:j+1)+A1(j,2)*B(j,1:j+1)+A1(j,3)*B(j+1,1:j+1);
end
C4(n,1:n )=A1(n,1)*B(n-1,1:n )+A1(n,2)*B(n,1:n );
f4 = flops;t4=toc;
E4 = C4-C1;
error4 = norm(E4,inf);
disp(sprintf(' %4.0f %7.0f %5.0f %7.0f %5.0f %8.3f %8.3f %8.3f %8.3f',n,f1,f2,f3,f4,t1,t2,t3,t4))
end