Chapter 7 M-File

M-File Home

Script Files

 ShowLSFit Fits a line to the square root function. ShowLSq Sensitivity of LS solutions. ShowSparseLS Examines "\"  LS solving with sparse arrays. ShowQR Illustrates QRrot. ShowBVP Solves a two-point boundary value problem. CholBench Benchmarks various Cholesky implementations. CholErr Sensitivity of symmetric positive definite systems. ShowCholBlock Explores block size selection in CholBlock.

Function Files

 Rotate Computes a rotation to zero bottom of 2-vector. QRrot QR factorization via rotations. LSq Least squares solution via QR factorization. CholTrid Cholesky factorization (tridiagonal version). CholTridSol Solves a factored tridiagonal system. CholScalar Cholesky factorization (scalar version). CholDot Cholesky factorization (dot product version). CholSax Cholesky factorization (saxpy version). CholGax Cholesky factorization (gaxpy version). CholRecur Cholesky factorization (recursive version). CholBlock Cholesky factorization (block version). MakeScalar Makes a block matrix a scalar matrix.

ShowLSFit

```% Script File: ShowLSFit
% Displays two LS fits to the function f(x) = sqrt(x) on [.25,1].

close all
z = linspace(.25,1);
fz = sqrt(z);
for m = [2 100 ]
x = linspace(.25,1,m)';
A = [ones(m,1) x];
b = sqrt(x);
xLS = A\b;
alpha = xLS(1);
beta  = xLS(2);
figure
plot(z,fz,z,alpha+beta*z,'--')
title(sprintf('m = %2.0f,  alpha = %10.6f,  beta = %10.6f',m,alpha,beta))
end
```

ShowLSq

```% Script File: Show LSq
% Illustrates LSq on problems with user-specified
% dimension and condition.
m = input('Enter m: ');
n = input('Enter n: ');
logcond = input('Enter the log of the desired condition number: ');
condA = 10^logcond;
[Q1,R1] = qr(randn(m,m));
[Q2,R2] = qr(randn(n,n));
A = Q1(:,1:n)*diag(linspace(1,1/condA,n))*Q2;
clc
disp(sprintf('m = %2.0f, n = %2.0f, cond(A) = %5.3e',m,n,condA))
b = A*ones(n,1);
[x,res] = LSq(A,b);
disp(' ')
disp('Nonzero residual problem.')
disp(' ')
disp('   Exact Solution          Computed Solution')
disp('-----------------------------------------')
for i=1:n
disp(sprintf(' %20.16f   %20.16f',1,x(i)))
end
b = b+Q1(:,m);
[x,res] = LSq(A,b);
disp(' ')
disp(' ')
disp(' ')
disp('Zero residual problem.')
disp(' ')
disp('   Exact Solution          Computed Solution')
disp('-----------------------------------------')
for i=1:n
disp(sprintf(' %20.16f   %20.16f',1,x(i)))
end
```

```% Script File: ShowSparseLS
% Illustrate sparse backslash solving for LS problems

clc
n = 50;
disp(' m      n   full A flops   sparse A flops')
disp('----------------------------------------')
for m = 100:100:1000
A = tril(triu(rand(m,m),-5),5);
p = m/n;
A = A(:,1:p:m);
A_sparse = sparse(A);
b = rand(m,1);
% Solve an m-by-n LS problem where the A matrix has about
% 10 nonzeros per column. In column j these nonzero entries
% are more or less A(j*m/n+k,j), k=-5:5.
flops(0)
x = A\b;
f1 = flops;
flops(0)
x = A_sparse\b;
f2 = flops;
disp(sprintf('%4d  %4d  %10d  %10d  ',m,n,f1,f2))
end
```

ShowQR

```% Script File: ShowQR
% Illustrates how the QR factorization is found via Givens Rotations.

m = 5;
n = 3;
disp(' ')
disp('Show upper triangularization of ')
A = rand(m,n);
for j=1:n
for i=m:-1:j+1
input('Strike Any Key to Continue');
clc
%Zero A(i,j)
Before = A
[c,s] = Rotate(A(i-1,j),A(i,j));
A(i-1:i,j:n) = [c s; -s c]  *A(i-1:i,j:n);
disp(sprintf('Zero A(%g,%g)',i,j))
After = A
end
end
disp('Upper Triangularization Complete.')
```

ShowBVP

```% Script File: ShowBVP
% Illustrates the numerical solution to
%
%        -D^2 u + u(x) = 2xsin(x) - 2cos(x)    u(0) = u(pi) = 0.
%
% Exact solution = u(x) = x*sin(x).

close all
n =  100;
x =  linspace(0,pi,n)';
hx = pi/(n-1);
d =  2*ones(n-2,1) + hx^2;
e = -ones(n-2,1);
[g,h] = CholTrid(d,e);
b = hx^2*( 2*x(2:n-1).*sin(x(2:n-1)) - 2*cos(x(2:n-1)));
umid  = CholTridSol(g,h,b);
u = [0;umid;0];
plot(x,u,x,x.*sin(x))
err = norm(u - x.*sin(x),'inf');
title('Solution to   -D^2 u + u = 2xsin(x) - 2cos(x), u(0)=u(pi)=0')
xlabel(sprintf(' n = %3.0f    norm(u - xsin(x),''inf'') = %10.6f',n,err))
```

CholBench

```% Script File: CholBench
% Benchmarks five different Cholesky implemntations.

clc
n = input('Enter n: ');
X = rand(2*n,n);
A = X'*X;
disp(' ');
disp(sprintf('n = %2.0f',n))
disp(' ');
disp('Algorithm   Relative Time     Flops')
disp('-----------------------------------')

flops(0); tic; G = CholScalar(A); t1 = toc; f = flops;
disp(sprintf('CholScalar      %6.3f        %5.0f',t1/t1,f));

flops(0); tic; G = CholDot(A);    t2 = toc; f = flops;
disp(sprintf('CholDot         %6.3f        %5.0f',t2/t1,f));

flops(0); tic; G = CholSax(A);    t3 = toc; f = flops;
disp(sprintf('CholSax         %6.3f        %5.0f',t3/t1,f));

flops(0); tic; G = CholGax(A);    t4 = toc; f = flops;
disp(sprintf('CholGax         %6.3f        %5.0f',t4/t1,f));

flops(0); tic; G = CholRecur(A);  t5 = toc; f = flops;
disp(sprintf('CholRecur       %6.3f        %5.0f',t5/t1,f));

disp(' ')
disp('Relative time = Time/CholScalar Time')
```

CholErr

```% Script File: CholErr
% Error when solving a Hilbert Matrix System.

clc
disp('  n     cond(A)       relerr ')
disp('----------------------------------')
for n = 2:12
A = hilb(n);
b = randn(n,1);
G = CholScalar(A);
y = LTriSol(G,b);
x = UTriSol(G',y);
condA = cond(A);
xTrue = invhilb(n)*b;
relerr = norm(x-xTrue)/norm(xTrue);
disp(sprintf(' %2.0f   %10.3e    %10.3e  ',n,condA,relerr))
end
```

ShowCholBlock

```% Script File: ShowCholBlock
% Effect of block size on performance of block Cholesky.

n = 192;
A = randn(n,n);
A = A'*A;
G0 = chol(A)';
clc
disp(sprintf('n =  %2.0f',n))
disp(' ')
disp(' Block Size   Time / Unblocked Time    ')
disp('---------------------------------------')
k =0;
for p=[ 8 12 16 24 32 48 96]
tic;
G = CholBlock(A,p);
k=k+1;
t(k) = toc;
disp(sprintf('     %2.0f       %6.3f  ',p,t(k)/t(1)))
end
```

Rotate

```  function [c,s] = Rotate(x1,x2)
% [c,s] = Rotate(x1,x2)
% x1 and x2 are real scalars and c and s is a cosine-sine pair so
% -s*x1 + c*x2 = 0.

if x2==0
c = 1;
s = 0;
else
if abs(x2)>=abs(x1)
cotangent = x1/x2;
s = 1/sqrt(1+cotangent^2);
c = s*cotangent;
else
tangent = x2/x1;
c = 1/sqrt(1+tangent^2);
s = c*tangent;
end
end
```

QRrot

```  function [Q,R] = QRRot(A)
% [Q,R] = QRRot(A)
% The QR factorization of an m-by-n matrix A. (m>=n).
% Q is m-by-m orthogonal and R is m-by-n upper triangular.

[m,n] = size(A);
Q = eye(m,m);
for j=1:n
for i=m:-1:j+1
%Zero A(i,j)
[c,s] = Rotate(A(i-1,j),A(i,j));
A(i-1:i,j:n) = [c s; -s c]*A(i-1:i,j:n);
Q(:,i-1:i) = Q(:,i-1:i)*[c s; -s c]';
end
end
R = triu(A);
```

LSq

```  function [xLS,res] = LSq(A,b)
% [xLS,res] = LSq(A,b)
% Solution to the LS problem min norm(Ax-b) where A is a full
% rank m-by-n matrix with m>=n and b is a column m-vector.
% xLS is the n-by-1 vector that minimizes the norm(Ax-b) and
% res = norm(A*xLS-b).

[m,n] = size(A);
for j=1:n
for i=m:-1:j+1
%Zero A(i,j)
[c,s] = Rotate(A(i-1,j),A(i,j));
A(i-1:i,j:n) = [c s; -s c]*A(i-1:i,j:n);
b(i-1:i) = [c s; -s c]*b(i-1:i);
end
end
xLS = UTriSol(A(1:n,1:n),b(1:n));
if m==n
res = 0;
else
res = norm(b(n+1:m));
end
```

```  function [g,h] = CholTrid(d,e)
% G = CholTrid(d,e)
% Cholesky factorization of a symmetric, tridiagonal positive definite matrix A.
% d and e are column n-vectors with the property that
% A = diag(d) + diag(e(2:n),-1) + diag(e(2:n),1)
%
% g and h are column n-vectors with the property that the lower bidiagonal
% G = diag(g) + diag(h(2:n),-1) satisfies A = GG^T.

n = length(d);
g = zeros(n,1);
h = zeros(n,1);
g(1) = sqrt(d(1));
for i=2:n
h(i) = e(i)/g(i-1);
g(i) = sqrt(d(i) - h(i)^2);
end
```

```  function x = CholTridSol(g,h,b)
% x = CholTridSol(g,h,b)
% Solves the linear system G*G'x = b where b is a column n-vector and
% G is a nonsingular lower bidiagonal matrix. g and h are column n-vectors
% with G = diag(g) + diag(h(2:n),-1).

n = length(g);
y = zeros(n,1);

% Solve Gy = b
y(1) = b(1)/g(1);
for k=2:n
y(k) = (b(k) - h(k)*y(k-1))/g(k);
end

% Solve G'x = y
x = zeros(n,1);
x(n) = y(n)/g(n);
for k=n-1:-1:1
x(k) = (y(k) - h(k+1)*x(k+1))/g(k);
end
```

CholScalar

```  function G = CholScalar(A)
% G = CholScalar(A)
% Cholesky factorization of a symmetric and positive definite matrix A.
% G is lower triangular so A = G*G'.

[n,n] = size(A);
G = zeros(n,n);
for i=1:n
% Compute G(i,1:i)
for j=1:i
s = A(j,i);
for k=1:j-1
s = s - G(j,k)*G(i,k);
end
if j < i
G(i,j) = s/G(j,j);
else
G(i,i) = sqrt(s);
end
end
end
```

CholDot

```  function G = CholDot(A)
% G = CholDot(A)
% Cholesky factorization of a symmetric and positive definite matrix A.
% G is lower triangular so A = G*G'.

[n,n] = size(A);
G = zeros(n,n);
for i=1:n
% Compute G(i,1:i)
for j=1:i
if j==1
s = A(j,i);
else
s = A(j,i) - G(j,1:j-1)*G(i,1:j-1)';
end
if j < i
G(i,j) = s/G(j,j);
else
G(i,i) = sqrt(s);
end
end
end
```

CholSax

```  function G = CholSax(A)
% G = CholSax(A)
% Cholesky factorization of a symmetric and positive definite matrix A.
% G is lower triangular so A = G*G'.

[n,n] = size(A);
G = zeros(n,n);
s = zeros(n,1);
for j=1:n
s(j:n) = A(j:n,j);
for k=1:j-1
s(j:n) = s(j:n) - G(j:n,k)*G(j,k);
end
G(j:n,j) = s(j:n)/sqrt(s(j));
end
```

CholGax

```  function G = CholGax(A)
% G = CholGax(A)
% Cholesky factorization of a symmetric and positive definite matrix A.
% G is lower triangular so A = G*G'.

[n,n] = size(A);
G = zeros(n,n);
s = zeros(n,1);
for j=1:n
if j==1
s(j:n) = A(j:n,j);
else
s(j:n) = A(j:n,j) - G(j:n,1:j-1)*G(j,1:j-1)';
end
G(j:n,j)  = s(j:n)/sqrt(s(j));
end
```

CholRecur

```  function G = CholRecur(A)
% G = CholRecur(A)
% Cholesky factorization of a symmetric and positive definite matrix A.
% G is lower triangular so A = G*G'.

[n,n] = size(A);
G = zeros(n,n);
if n==1
G = sqrt(A);
else
G(1:n-1,1:n-1) = CholRecur(A(1:n-1,1:n-1));
G(n,1:n-1)     = LTriSol(G(1:n-1,1:n-1),A(1:n-1,n))';
G(n,n)         = sqrt(A(n,n) - G(n,1:n-1)*G(n,1:n-1)');
end
```

```  function G = CholBlock(A,p)
% G = CholBlock(A,p)
% Cholesky factorization of a symmetric and positive definite n-by-n matrix A.
% G is lower triangular so A = G*G'. p is the block size and must divide n.

% Represent A and G as m-by-m block matrices where m = n/p.
[n,n] = size(A);
m = n/p;
A = MakeBlock(A,p);
G = cell(m,m);
for i=1:m
for j=1:i
S = A{i,j};
for k=1:j-1
S = S - G{i,k}*G{j,k}';
end
if j < i
G{i,j} = (G{j,j}\S')';
else
G{i,i} = CholScalar(S);
end
end
end
% Convert G to a scalar matrix.
G = MakeScalar(G);
```

```  function A = MakeScalar(A_block)
% A = MakeScalar(A_block)
% Represents the m-by-m block matrix A_block as an n-by-n matrix of scalrs with
% where each block is p-by-p and n=mp.

[m,m] = size(A_block);
[p,p] = size(A_block{1,1});
for i=1:m
for j=1:m
if ~isempty(A_block{i,j})
A(1+(i-1)*p:i*p,1+(j-1)*p:j*p) = A_block{i,j};
end
end
end
```