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.

 


CircBench

% 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

ProdBound

% 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

Show2DQuad

% 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

FFTflops

% 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

Circulant1


  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

Circulant2

  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

MatVecRO

  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

MatVecCO

  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

MatMatDot

  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

MatMatSax

  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

MatMatVec


  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);

wCompNC

  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