Assignment 1 Scripts

P1

% P1: Matrix Products

global OpCount           % For flop counting

% Experiment Parameters (May want to vary these as you develop your MatProd implementation.)

t = 20;                  % The number of matrices in the product.
SizeLimit = 100;         % Each matrix dimension is randomly picked from the interval [1,SizeLimit].
NProbs = 500;            % The number of examples that are run.

% So we all do the same 'random' problems...
rand('seed',0)
randn('seed',0)

error = []; quotient = [];
for eg = 1:NProbs
    
   % Generate a random set of matrix dimensions
   z = ceil(SizeLimit*rand(t+1,1)); m = z(1:t); n = z(2:t+1);
   
   % Set up a cell array of random matrices
   A = cell(t,1);
   for k=1:t
      A{k} = randn(m(k),n(k));
   end
   
   % The left-to-right method
   B = A{1};
   OpCount0 = 0;
   for k=2:t
      B = B*A{k};
      OpCount0 = OpCount0 + 2*m(1)*m(k)*n(k);
   end
   
   % The Simple Divide and Conquor Method
   OpCount = 0;
   C = MatProd(A);
   
   % Make a note of the error and the flop ratio
   error = [error norm(B-C,1)/norm(B,1)];
   quotient = [quotient OpCount/OpCount0];
   
end

% Graphicaly display the results of the simulation ...

close all
plot(quotient)
title('Ratio of D&C flops to Left-to-Right Flops')
xlabel('Example Index')
ylabel('Ratio')

figure
edges = [0:.1:1];
N = histc(quotient,edges);
favorable = round(100*(sum(quotient<=1)/NProbs));
bar(edges,N,'histc')
axis([0 1 0 1.2*max(N)])
title(sprintf('%2d percent of the ratios are favorable. Their distribution...',favorable)) 
ylabel('Frequency')
xlabel('Flops of D&C Method  /  Flops of Left-to-Right Method')
m1= sprintf('t = %1d',t);
m2 = sprintf('SizeLimit = %1d',SizeLimit);
m3 = sprintf('NProbs = %1d',NProbs);
text(.7,max(N),{m1, m2, m3 } )

figure
semilogy(error)
title('Relative error of D&C method assuming Left-to-Right method is exact.')
ylabel('Relative Error')
xlabel('Example Index')

P3

% P3: Block Matrix Manipulation

% Tests the Hamiltonian squared function HamSq

clc
randn('seed',0)
close all

% A small example to check for correctness..

n = 4;
A = randn(n,n);
F = randn(n,n); F = F+F';
G = randn(n,n); G = G+G';
H = [A F ; G  -A'];
KNaive = H*H
K = HamSq(H)

disp(' ')
% Fill the following comment so it specifies the right quotient for your implementation
disp('(Naive Method Flops) / (HamSq Method Flops) = XXX')

% Benchmark against the naive method...

ratio = [];

% nProbs = number of problems per timing. Should be large enough so that
% the benchmark on your computer is not affected by clock granularity.
nProbs = 5;

% nVals = problem dimensions
nVals = 200:50:500;

% You may want to reduce nProbs and nVals while you are
% still developing your HamSq implementation.

for n = nVals
   % Generate a problem
   A = randn(n,n);
   F = randn(n,n); F = F+F';
   G = randn(n,n); G = G+G';
   H = [A F ; G  -A'];
   
   % Time HamSq 
   tic
   for k=1:nProbs
       K = HamSq(H);
   end
   tHamSq = toc;
   
   % Time the Naive Method
   tic
   for k=1:nProbs
       K = H*H;
   end
   tNaive = toc;
   ratio = [ratio tNaive/tHamSq];
end
plot(nVals,ratio)
xlabel('n')
ylabel('tNaive / tHamSq')
title('Comparing the time for K = H^2 with time for K = HamSq(H)')
   

P4

% P4: Perturbing Stochastic Matrices

% Generate a small example, display, and report its condition number.

rand('seed',0)
clc
n = 5;
A = rand(n,n);
for k=1:n
    A(:,k) = A(:,k)/norm(A(:,k),1);
end
A = A
kappa = cond(A)
disp(' ')

%
q = 2;
for p=1:n
   theta = MakeSing(A,p,q);
   mu = (1 - theta*A(p,q))/(1 - A(p,q));
   Apq = A;
   Apq(:,q) = [mu*A(1:p-1,q) ; theta*A(p,q) ; mu*A(p+1:n,q)];
   % Check column sums and condition
   error = norm(ones(1,n)*Apq - ones(1,n));
   kappa0 = cond(Apq);
   disp(sprintf('p = %1d, q = %1d, error = %5.1e, cond = %5.1e',p,q,error,kappa0))
   Apq = Apq
   disp(' ')
end