Assignment 1 Scripts

P1

% Script P1
% Tests PolyVec

% Set seed so we all do the same random problems..

randn('state',0);

% Correctness...
p = 3; n = 15; m = 4;
A = triu(randn(n,n),-p); c = randn(m,1);
M = (c(1)*eye(n,n) + c(2)*A + c(3)*A^2 + c(4)*A^3);
y0 = M(:,1);
y = PolyVec(c,A,p);
error = norm(y - y0);
clc
disp(sprintf('Error = %10.4e',error))

% Efficiency...
n = 500; m = 20;
A0 = randn(n,n); 
c = randn(m,1);

% Repeat factor for more reliable benchmarks
% Might want a smaller value while refining your code.

nRepeat = 100; 

% Get a calibration time...
x = randn(n,1);
tic 
for j=1:nRepeat
   y = A0*x;
end; 
tbase = toc;
disp(' ')
disp('n = 500 Example')
disp(' ')
disp('  p    Normalized time')
disp('--------------------------')
for p=5:5:20
   A = triu(A0,-p);
   tic; 
   for j=1:nRepeat
      y = PolyVec(c,A,p);
   end
   t = toc;
   disp(sprintf('%3d   %10.3f',p,t/tbase))
end

P2

% Script P2
% Tests KrylovProd

clc
n = 100;
randn('seed',0)

% Get a random orthogonal matrix Q vis at the qr factorization..

[Q,R] = qr(randn(n,n));

v = randn(n,1);
x = randn(n,1);

% Explicit M method...
M = v;
for k=2:n
   M = [M Q*M(:,k-1)];
end
y0 = M'*(M*x);

% Compare...
y = KrylovProd(Q,v,x);
disp(sprintf('KrylovProd Error = %10.3e',norm(y-y0)))

P3

% Script P3
% Tests pNormImage

close all
A = [2 5;3 -3];
% Look at what A does to unit p-norm vectors for p=1,2,10, and 30..
for p=[1 2 10 30]
   figure
   pNormImage(A,p);
   hold on
   plot([-8 8],[0 0],'--')
   plot([0 0],[-8 8],'--')
   axis([-8 8 -8 8])
   title(sprintf('A = [%3.0f %3.0f ; %3.0f %3.0f]   p = %3.0f',A,p))
end

P4

% Script P4
% Tests FastProd

t = 9;

% Generate W_{2^t} explicitly
W = 1;
for k=1:9
   W = [ W W; W -W];
end

clc
randn('seed',0)
x = randn(2^t,1);

% Compare the explicit-W method with FastProd...
y0 = W*x;
y = FastProd(x);
disp(sprintf('FastProd Error = %10.3e',norm(y-y0)))