function s6a3(n)

% test input, set default
if nargin < 1 | n<=2 | mod(n,1)~=0 n=300; disp('using default'); end

% b) compare sparse and full computations
rand(0);	
x=rand(n,1);	% x has random entries in (0,1)
disp('my sparsemult');
tic;
ysm = sparsemult(x,n);
toc;

disp('full matrix');
tic;
Af = diag(1:n) + diag(1:n-1,1) + diag(1:n-1,-1);
yf = Af*x;
toc;

disp('MATLAB sparse');
tic;
As = sparse(diag(1:n) + diag(1:n-1,1) + diag(1:n-1,-1));
ys = As*x;
toc;

% c) Vector Iteration to compute largest eigenvalue (absolut value)
x = eye(n,1); 		% starting value for iteration
l = 0; l0 = 1;
tol = 1.e-10;		% stopping tolerance
c  = 0;			% count iterations
while (abs(l-l0)>tol)
  c  = c+1;
  l0 = l; x0 = x;	% save previous values
  x  = sparsemult(x,n); % vector iteration x_k --> x_{k+1} = A*x_k
  l = x0'*x/(x0'*x0);	% Rayleigh quotient
  x  = x/norm(x);
end
disp(['tol=',num2str(tol,4),', after ',int2str(c),...
      ' iterations: \lamda_{max}=',num2str(l,12)]);
if (n<10) max(abs(eig(Af))), else OPTIONS.disp=0; eigs(As,1,OPTIONS), end


%%%%%%%%%%%%%%% SPARSE MATRIX VECTOR MULTIPLICATION %%%%%%%%%%%

function y=sparsemult(x,n)

% computes the sparse Matrix-Vector-Product y=Ax for 
% A=diag(1:n) + diag(1:n-1,1) + diag(1:n-1:-1) 
%                / d_1*x_1 + d_1*x_2,                    i = 1
% y_i = (Ax)_i = | d_{i-1}*x_{i-1} + d_i*(x_i+x_{i+1}),  2 <= i <= n-1
%                \ d_{n-1}*x_{n-1} + d_n*x_n,            i = n
% in:	x	vector to be multiplied
% 	n	A is in n x n, n >= 2
% out:	y	on output y = Ax

d = [1:n]';
y = d.*x + [d(1:n-1).*x(2:n);0] + [0;d(1:n-1).*x(1:n-1)];


%%% OUTPUT %%%
%  >> s6a3(300);  
%  my sparsemult
%  elapsed_time =
%     0.00229100000000
%  full matrix
%  elapsed_time =
%     0.06270100000000
%  MATLAB sparse
%  elapsed_time =
%     0.04803500000000
%  tol=1e-10, after 933 iterations: \lamda_{max}=869.641106459
%  ans =
%       8.696411064606707e+02
