%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  This program solves Problem 1 in
%  Your Homework Assignment
%  Computer Memory and Arithmetic:
%  A Look under the Hood
%  Dianne P. O'Leary
%  problem1.m    02/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We define a random matrix A of size 1024x1024, and a random vector x.
% We initialize two vectors to 0, to allocate space for the product Ax,
% computed in two ways.
%
% We compute y1 using the "inner product" formulation of the algorithm,
% computing each component of y1 as the inner product of a row of A
% with the vector x.
%
% We compute y2 using the "saxpy" formulation of the algorithm,
% successively adding columns of A weighted by elements of x.
%
% We time the two operations, using tic/toc and using cputime,
% and then display the resulting times and the differences in the
% computed products.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 1024;

A = rand(n,n);
x = rand(n,1);
y1 = zeros(n,1);
y2 = zeros(n,1);

t0 = cputime;

for i=1:n,
   y1(i) = A(i,:)*x;
end

t1 = cputime;

for i=1:n,
   y2 = y2 + x(i)*A(:,i);
end

t2 = cputime - t1;
t1 = t1 - t0;

y1 = zeros(n,1);
y2 = zeros(n,1);

tic

for i=1:n,
   y1(i) = A(i,:)*x;
end

time1 = toc;

tic

for i=1:n,
   y2 = y2 + x(i)*A(:,i);
end

time2 = toc;

disp(sprintf('Elapsed time for inner product algorithm: %f sec',time1))
disp(sprintf('Elapsed time for    saxpy      algorithm: %f sec',time2))
disp(sprintf('  '))
disp(sprintf('CPU     time for inner product algorithm: %f sec',t1))
disp(sprintf('CPU     time for    saxpy      algorithm: %f sec',t2))
disp(sprintf('  '))
disp(sprintf('Difference in the computed answers:       %e',norm(y1-y2)))