Chapter 5 M-Files
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
% 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
% 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
% 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
% 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
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
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
% 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)
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
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
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);
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;
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
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
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
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