% Script P4
% Tests Sprod

% Generate a random example...
  n = 1000;
  randn('seed',0)
  d = randn(n,1);
  u = randn(n,1);
  v = randn(n,1);
  x = randn(n,1);

% The conventional method...
  A = diag(d) + tril(u*v',-1) + triu(v*u',1);
  yc = A*x;
  
% The fast method...
  [y,nFlops] = Sprod(d,u,v,x);
  
% Display results
  clc
  relErr = norm(yc-y)/norm(yc);
  fprintf('relErr = %10.3e\n\nnFlops = %1d\n\n',relErr,nFlops)
  format long
  First5Components = y(1:5)