Chapter 5 M-Files

M-File Home

Script Files

 CircBench Benchmarks Circulant1 and Circulant2. MatBench Benchmarks various matrix-multiply methods. AveNorms Compares various norms on random vectors. ProdBound Examines the error in three-digit matrix multiplication. ShowContour Displays various contour plots of SampleF. Show2DQuad Illustrates CompQNC2D. FFTflops Compares FFT and DFT flops. StrassFlops Examines Strassen multiply.

Function Files

 Circulant1 Scalar-level circulant matrix set-up. Circulant2 Vector-level circulant matrix set-up. MatVecRO Row-oriented matrix-vector product. MatVecCO Column-Oriented matrix-vector product. MatMatDot Dot-product matrix-matrix product. MatMatSax Saxpy matrix-matrix product. MatMatVec Matrix-vector matrix-matrix product. MatMatOuter Outer product matrix-matrix product. MakeBlock Makes a cell array representation of a block matrix. ShowSparse Illustrates sparse. ShowNonZeros Displays the sparsity pattern of a matrix. Prod3Digit Three-digit matrix-matrix product. SampleF A Gaussian type function of two variables. CompQNC2D Two-dimensional Newton-Cotes rules. wCompNC Weight vector for composite Newton-Cotes rules. SampleF2 A function of two variables with strong local maxima. DFT Discrete Fourier transform. FFTRecur A recursive radix-2 FFT. Strass Recursive Strassen matrix multiply.

% Script File: CircBench
% Benchmarks Circulant1 and Circulant2.

clc
nRepeat = 20;
disp('t1 = time for Circulant1.')
disp('t2 = time for Circulant2.')
disp(' ')
disp('  n  t1/t2 ')
disp('----------------')
for n=[100 200 400]
a = 1:n;
tic;  for k=1:nRepeat, C = Circulant1(a); end, t1 = toc;
tic;  for k=1:nRepeat, C = Circulant2(a); end, t2 = toc;
disp(sprintf(' %4.0f   %5.3f',n,t1/t2))
end


MatBench

% Script File: MatBench
% Compare different matrix multiplication methods.

clc
nRepeat = 10;
disp(' n     Dot     Saxpy    MatVec   Outer    Direct')
disp('------------------------------------------------')
for n = [100 200 400]
A = rand(n,n);
B = rand(n,n);
tic; for k=1:nRepeat,C = MatMatDot(A,B);   end, t1 = toc/nRepeat;
tic; for k=1:nRepeat,C = MatMatSax(A,B);   end, t2 = toc/nRepeat;
tic; for k=1:nRepeat,C = MatMatVec(A,B);   end, t3 = toc/nRepeat;
tic; for k=1:nRepeat,C = MatMatOuter(A,B); end, t4 = toc/nRepeat;
tic; for k=1:nRepeat,C = A*B;              end, t5 = toc/nRepeat;
disp(sprintf('%2.0f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f ',n,t1,t2,t3,t4,t5))
end


AveNorms

% Script File: AveNorms
% Examines average value of norms.

clc
disp('r1 = norm(x,1)/norm(x,''inf'')')
disp('r2 = norm(x,2)/norm(x,''inf'')')
disp(' ')
disp('    n          r1           r2')
disp('----------------------------------')
r = 100;
for n = 10:10:100
s1 = 0;
s2 = 0;
for k=1:r
x = randn(n,1);
s1 = s1 + norm(x,1)/norm(x,'inf');
s2 = s2 + norm(x,2)/norm(x,'inf');
end
disp(sprintf(' %4d    %10.3f   %10.3f',n,s1/r,s2/r))
end


% Script File: ProdBound
% Examines the error in 3-digit matrix multiplication.

clc
eps3 = .005;       % 3-digit machine precision
nRepeat = 10;      % Number of trials per n-value
disp('   n   1-norm factor ')
disp('------------------------')
for n = 2:10
s = 0;
for r=1:nRepeat
A = randn(n,n);
B = randn(n,n);
C = Prod3Digit(A,B);
E = C - A*B;
s = s+ norm(E,1)/(eps3*norm(A,1)*norm(B,1));
end
disp(sprintf('%4.0f     %8.3f   ',n,s/nRepeat))
end


ShowContour

% Script File: ShowContour
% Illustrates various contour plots.

close all
% Set up array of function values.
x = linspace(-2,2,50)';
y = linspace(-1.5,1.5,50)';
F = SampleF(x,y);

% Number of contours set to default value:
figure
Contour(x,y,F)
axis equal

% Five contours:
figure
contour(x,y,F,5);
axis equal

% Five contours with specified values:
figure
contour(x,y,F,[1 .8  .6  .4  .2])
axis equal

% Four contours with manual labeling:
figure
c= contour(x,y,F,4);
clabel(c,'manual');
axis equal


% Script File: Show2DQuad
% Integral of SampleF2 over [0,2]x[0,2] for various 2D composite
% QNC(7) rules.

clc
m = 7;
disp(' Subintervals    Integral            Time')
disp('------------------------------------------------')
for n = [32 64 128 256]
tic, numI2D = CompQNC2D('SampleF2',0,2,0,2,m,n,m,n); time = toc;
disp(sprintf(' %7.0f  %17.15f  %11.4f',n,numI2D,time))
end


% Script File: FFTflops
% Compares the ordinary DFT with the FFT.

clc
x = randn(1024,1)+sqrt(-1)*randn(1024,1);
disp('   n        DFT       FFTRecur      FFT ')
disp('           Flops       Flops       Flops')
disp('------------------------------------------')
for n = [2 4 8 16 32 64 128 256 512 1024]
flops(0); y = dft(x(1:n));       F1 = flops;
flops(0); y = FFTRecur(x(1:n));  F2 = flops;
flops(0); y = fft(x(1:n));       F3 = flops;
disp(sprintf(' %4d %10d %10d  %10d', n,F1,F2,F3))
end


StrassFlops

% Script File: StrassFlops
% Compares Strassen method flops with ordinary A*B flops.

clc
nmin = 32;
nmax = 1024;
A0 = rand(nmax,nmax);
B0 = rand(nmax,nmax);
disp('  n   Strass Flops/Ordinary Flops')
disp('-------------------------------------')
n = 32;
while (n<=nmax)
A = A0(1:n,1:n);
B = B0(1:n,1:n);
flops(0);
C = Strass(A,B,nmin);
f = flops;
disp(sprintf('  %4d       %6.3f', n,f/(2*n^3)))
n = 2*n;
end


  function C = Circulant1(a)
% C = Circulant1(a) is a circulant matrix with first row equal to a.

n = length(a);
C = zeros(n,n);
for i=1:n
for j=1:n
C(i,j) = a(rem(n-i+j,n)+1);
end
end

  function C = Circulant2(a)
% C = Circulant2(a) is a circulant matrix with first row equal to a.

n = length(a);
C = zeros(n,n);
C(1,:) = a;
for i=2:n
C(i,:) = [ C(i-1,n) C(i-1,1:n-1) ];
end


  function y = MatVecR0(A,x)
% y = MatVecRO(A,x)
% Computes the matrix-vector product y = A*x (via saxpys) where
% A is an m-by-n matrix and x is a column n-vector.

[m,n] = size(A);
y = zeros(m,1);
for i=1:m
y(i) = A(i,:)*x;
end


  function y = MatVecC0(A,x)
% y = MatVecCO(A,x)
% This computes the matrix-vector product y = A*x (via saxpys) where
% A is an m-by-n matrix and x is a columnn-vector.

[m,n] = size(A);
y = zeros(m,1);
for j=1:n
y = y + A(:,j)*x(j);
end


  function C = MatMatDot(A,B)
% C = MatMatDot(A,B)
% This computes the matrix-matrix product C =A*B (via dot products) where
% A is an m-by-p matrix, B is a p-by-n matrix.

[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for j=1:n
% Compute j-th column of C.
for i=1:m
C(i,j) = A(i,:)*B(:,j);
end
end


  function C = MatMatSax(A,B)
% C = MatMatSax(A,B)
% This computes the matrix-matrix product C = A*B (via saxpys) where
% A is an m-by-p matrix, B is a p-by-n matrix.

[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for j=1:n
% Compute j-th column of C.
for k=1:p
C(:,j) = C(:,j) + A(:,k)*B(k,j);
end
end


  function C = MatMatVec(A,B)
% C = MatMatVec(A,B)
% This computes the matrix-matrix product C = A*B (via matrix-vector products)
% where A is an m-by-p matrix, B is a p-by-n matrix.

[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for j=1:n
% Compute j-th column of C.
C(:,j) = C(:,j) + A*B(:,j);
end


MatMatOuter

  function C = MatMatOuter(A,B)
% C = MatMatOuter(A,B)
% This computes the matrix-matrix product C = A*B (via outer products) where
% A is an m-by-p matrix, B is a p-by-n matrix.

[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for k=1:p
% Add in k-th outer product
C = C + A(:,k)*B(k,:);
end


MakeBlock

  function A = MakeBlock(A_scalar,p)
% A = MakeBlock(A_scalar,p)
% A_scalar is an n-by-n matrix and p divides n. A is an (n/p)-by-(n/p)
% cell array that represents A_scalar as a block matrix with p-by-p blocks.

[n,n] = size(A_scalar);
m = n/p;
A = cell(m,m);
for i=1:m
for j=1:m
A{i,j} = A_scalar(1+(i-1)*p:i*p,1+(j-1)*p:j*p);
end
end


ShowSparse

% Script File: ShowSparse
% Looks at the effect of the sparse operation on tridiagonal matrix products.

close all
nvals = 100:100:1000;
repeatVal = [200 200  100  50  20  10  10  10  10 5];
for k=1:length(nvals)
n = nvals(k); x = randn(n,1);
A = diag(2*(1:n)) - diag(1:n-1,-1) + diag(1:n-1,1);
S = sparse(A);
flops(0), y = A*x;  fullA_flops(k)   = flops;
flops(0), y = S*x;  sparseA_flops(k) = flops;
end
semilogy(nvals,fullA_flops,nvals,sparseA_flops)
xlabel('n')
title('(Tridiagonal A) \times (Vector)')
legend('Flops with Full A','Flops with Sparse A',2)


ShowNonZeros

  function ShowNonZeros(A)
% Displays the sparsity structure of an n-by-n matrix A.
% It is best that n<=30.

close all
[n,n] = size(A);
s = 1/n;
m = round(200/n);
figure
axis([0 1 0 1])
axis ij equal off
HA = 'HorizontalAlignment';
VA = 'VerticalAlignment';
for i=1:n
for j=1:n
if A(i,j)==0
text(i*s,(j-.3)*s,'.',HA,'center');
else
text(i*s,j*s,'X','FontWeight','bold','FontSize',m,HA,'center');
end
end
end


Prod3Digit

  function C = Prod3Digit(A,B)
% C = Prod3Digit(A,B)
% This computes the matrix product C = A*B in 3-digit floating point
% arithmetic. A is m-by-p and B is p-by-n.

[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for i=1:m
for j=1:n
s = represent(0);
for k=1:p
aik = Represent(A(i,k));
bkj = Represent(B(k,j));
s   = Float(s,Float(aik,bkj,'*'),'+');
end
C(i,j) = Convert(s);
end
end


SampleF

  function F = SampleF(x,y)
% x is a column n-vector, y is a column m-vector and
% F  is an m-by-n matrix with F(i,j) = exp(-(x(j)^2 + 3y(i)^2)).

n = length(x);
m = length(y);
A = -((2*y.^2)*ones(1,n) + ones(m,1)*(x.^2)')/4;
F = exp(A);


CompQNC2D

  function numI2D = CompQNC2D(fname,a,b,c,d,mx,nx,my,ny)
% numI2D = CompQNC2D(fname,a,b,c,d,mx,nx,my,ny)
%
% fname is a string that names a function of the form f(x,y).
% If x and y are vectors, then it should return a matrix F with the
% property that F(i,j) = f(x(i),y(j)), i=1:length(x), j=1:length(y).
%
% a,b,c,d are real scalars.
% mx and my are integers that satisfy 2<=mx<=11, 2<=my<=11.
% nx and ny are positive integers
%
% numI2D  approximation to the integral of f(x,y) over the rectangle [a,b]x[c,d].
% The compQNC(mx,nx) rule is used in the x-direction and the compQNC(my,ny)
% rule is used in the y-direction.

[omega,x] = wCompNC(a,b,mx,nx);
[mu,y]    = wCompNC(c,d,my,ny);
F = feval(fname,x,y);
numI2D = (b-a)*(d-c)*(omega'*F*mu);


  function [omega,x] = wCompNC(a,b,m,n)
% [omega,x] = wCompNC(a,b,m,n)
%
% m is an integer with 2<=m<=11 and n is a positive integer.
%
% omega is a column vector of weights for the composite m-point
% Newton-Cotes rule with n equal-length subintervals across [a,b]
% x is column vector of abscissas.

N = n*(m-1)+1;
x = linspace(a,b,N)';
w = NCWeights(m);
omega = zeros(N,1);
for k=1:(m-1):N-m+1
omega(k:k+m-1) = omega(k:k+m-1) + w;
end
omega = omega/n;


SampleF2

  function F = SampleF2(x,y)
% A function of 2 variables with strong local maxima at (.3,.4) and (.7,.3).

F = zeros(length(y),length(x));
for j=1:length(x)
F(:,j) =  1./( ((x(j)-.3).^2 + .1)*((y-.4).^2 + .1)) + 1./( ((x(j)-.7).^2 + .2)*((y-.3).^2 + .2)) - 6;
end


DFT

  function y = DFT(x)
% y = DFT(x)
% y is the discrete Fourier transform of a column n-vector x.

n = length(x);
y = x(1)*ones(n,1);
if n > 1
v = exp(-2*pi*sqrt(-1)/n).^(0:n-1)';
for k=2:n
z = rem((k-1)*(0:n-1)',n ) +1;
y = y + v(z)*x(k);
end
end


FFTRecur

  function y = FFTRecur(x)
% y = FFTRecur(x)
% y is the discrete Fourier transform of a column n-vector x where
% n is a power of two.

n = length(x);
if n ==1
y = x;
else
m = n/2;
yT = FFTRecur(x(1:2:n));
yB = FFTRecur(x(2:2:n));
d = exp(-2*pi*sqrt(-1)/n).^(0:m-1)';
z = d.*yB;
y = [ yT+z ; yT-z ];
end


Strass

  function C = Strass(A,B,nmin)
% C = Strass(A,B,nmin)
% This computes the matrix-matrix product C = A*B (via the Strassen Method) where
% A is an n-by-n matrix, B is a n-by-n matrix and n is a power of two. Conventional
% matrix multiplication is used if n < nmin where nmin is a positive integer.

[n,n] = size(A);
if n < nmin
C = A*B;
else
m = n/2; u = 1:m; v = m+1:n;
P1 = Strass(A(u,u)+A(v,v),B(u,u)+B(v,v),nmin);
P2 = Strass(A(v,u)+A(v,v),B(u,u),nmin);
P3 = Strass(A(u,u),B(u,v)-B(v,v),nmin);
P4 = Strass(A(v,v),B(v,u)-B(u,u),nmin);
P5 = Strass(A(u,u)+A(u,v),B(v,v),nmin);
P6 = Strass(A(v,u)-A(u,u),B(u,u) + B(u,v),nmin);
P7 = Strass(A(u,v)-A(v,v),B(v,u)+B(v,v),nmin);
C = [ P1+P4-P5+P7   P3+P5; P2+P4 P1+P3-P2+P6];
end