Assignment 1  Test Scripts

P1

% Script P1
% Examines PowerDiff

% Correctness
clc
A = randn(20,20); u = randn(20,1); v = randn(20,1);
[X,Y] = PowerDiff(A,u,v,3);
error = norm(A^3 + X*Y' - (A+u*v')^3);
disp(sprintf('Error in an n=20, k=3 example: %10.3e',error))
    
% Benchmark
n = 1000;
A = randn(n,n);u = randn(n,1);v = randn(n,1);
nRepeat = 20;
tic
for i=1:nRepeat
   z= A*u;
end
t1 = toc;
disp(' ')
disp('  k    Time(PowerDiff)/Time(Mat*vec) ')
disp('-----------------------------------------------')
for k=[2 4 8 16 32]
    tic
    for i=1:nRepeat
       [X,Y] = PowerDiff(A,u,v,k);
    end
    t2 = toc;
    ratio = t2/t1;
    disp(sprintf('  %2d       %6.2f   ',k,ratio))
end

P2

% P2
% Checks ThreeWayProd

clc

% Correctness

m = 10; A = randn(m,m);
n = 20; B = randn(n,n);
p = 5;  C = randn(p,p);
x = randn(m*n*p,1);
y0 = kron(A,kron(B,C))*x;
y = ThreeWayProd(A,B,C,x);

RelError = norm(y - y0)/norm(y0);
disp(sprintf('Relative Error = %10.3e',RelError))

% Compare ThreeWayProd with naive algorithm...

n = 15;
A = randn(n,n); B = randn(n,n); C = randn(n,n); x = randn(n^3,1);
nRepeat = 50; % For reliability given clock granularity
tic 
for i = 1:nRepeat
   y = kron(A,kron(B,C))*x;
end
t1 = toc;
tic 
for i=1:nRepeat
   y = ThreeWayProd(A,B,C,x);
end
t2 = toc;
disp(sprintf('\nTime(Naive) / Time(ThreeWayProd) = %5.1f  (n=%1d)',t1/t2,n))

% Increase n by factor of 10 and see how work changes...

n1 = 10*n;
A = randn(n1,n1); B = randn(n1,n1); C = randn(n1,n1); x = randn(n1^3,1);
tic
y = ThreeWayProd(A,B,C,x);
t3 = toc;
ratio = nRepeat*t3/t2;
disp(sprintf('\nIf n is multiplied by 10 then work goes up by a factor of %6.1f',ratio))