B9

Called Functions

MyDetLU
MyInv
bslashtxNP
bslashtxP
lutx

A33

[L,U,p,sig]=lutx(randn(5))

for i=7:9
   A=randn(i);
   norm(det(A)-MyDetLU(A))
end

L =

    1.0000         0         0         0         0
   -0.6893    1.0000         0         0         0
    0.5872   -0.7820    1.0000         0         0
   -0.4469   -0.3260    0.5504    1.0000         0
   -0.5923    0.0658    0.3511   -0.2129    1.0000


U =

   -1.3504   -0.2036    0.6520   -0.5767   -0.7348
         0   -1.5601    1.1501   -0.2964   -0.3156
         0         0    1.3929   -0.2182   -1.0427
         0         0         0    1.4418    0.8682
         0         0         0         0    0.9248


p =

     4
     3
     2
     1
     5


sig =

     1


ans =

   2.6645e-15


ans =

   8.8818e-15


ans =

   1.0658e-14

A34

for i = 6:9
   A=randn(i);
   norm(inv(A)-MyInv(A))
end

ans =

   1.3623e-15


ans =

   1.0728e-15


ans =

   2.9190e-13


ans =

   1.7541e-14

A35

[L,U] = lutx(randn(4))

for i = 6:9
   A=randn(i);
   [L,U] = lutx(A);
   [L2,U2] = lu(A);
   norm(L-L2)
end

L =

   -0.4216    1.0000         0         0
   -0.4804   -0.5037    0.6771    1.0000
    1.0000         0         0         0
    0.1985   -0.7455    1.0000         0


U =

   -1.9737    0.5152    0.3154    0.1098
         0   -0.4595    1.2739    0.4750
         0         0    1.8326   -1.3224
         0         0         0   -0.1372


ans =

   3.4953e-16


ans =

   4.0917e-16


ans =

   2.7859e-16


ans =

   3.0588e-16

A36

A=[10 -7 0;-3 2.099 6;5 -1 5];
b=[7 3.901 6]';

for i=6:-1:3
   digits(i);
   xP = bslashtxP(A,b)
   xNP= bslashtxNP(A,b)
end

function det = MyDetLU(A)
    [L,U,p,sig] = lutx(A);
    det=sig*prod(diag(U));
end

function x = MyInv(A)
% MyInv

n = length(A);
x = zeros(n);
E = eye(n);

% Triangular factorization
[L,U,p] = lutx(A);

for i=1:n
    b=E(:,i);
    % Permutation and forward elimination
    y = forward(L,b(p));

    % Back substitution
    x(:,i) = backsubs(U,y);
end

% ------------------------------

function x = forward(L,x)
% FORWARD. Forward elimination.
% For lower triangular L, x = forward(L,b) solves L*x = b.
[n,n] = size(L);
x(1) = x(1)/L(1,1);
for k = 2:n
   j = 1:k-1;
   x(k) = (x(k) - L(k,j)*x(j))/L(k,k);
end


% ------------------------------

function x = backsubs(U,x)
% BACKSUBS.  Back substitution.
% For upper triangular U, x = backsubs(U,b) solves U*x = b.
[n,n] = size(U);
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
   j = k+1:n;
   x(k) = (x(k) - U(k,j)*x(j))/U(k,k);
end

function x = bslashtxNP(A,b)
% BSLASHTX  Solve linear system (backslash)
% x = bslashtx(A,b) solves A*x = b



[n,n] = size(A);
if isequal(triu(A,1),zeros(n,n))
   % Lower triangular
   x = forward(A,b);
   return
elseif isequal(tril(A,-1),zeros(n,n))
   % Upper triangular
   x = backsubs(A,b);
   return
elseif isequal(A,A')
   [R,fail] = chol(A);
   if ~fail
      % Positive definite
      y = forward(R',b);
      x = backsubs(R,y);
      return
   end
end

% Triangular factorization
[L,U] = (lutxNP(A));

% Permutation and forward elimination
y = (forward(L,b));

% Back substitution
x = backsubs(U,y);


% ------------------------------

function x = forward(L,x)
% FORWARD. Forward elimination.
% For lower triangular L, x = forward(L,b) solves L*x = b.
[n,n] = size(L);
x(1) = x(1)/L(1,1);
for k = 2:n
   j = 1:k-1;
   x(k) = vpa(vpa(x(k) - vpa(L(k,j)*x(j)))/L(k,k));
end


% ------------------------------

function x = backsubs(U,x)
% BACKSUBS.  Back substitution.
% For upper triangular U, x = backsubs(U,b) solves U*x = b.
[n,n] = size(U);
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
   j = k+1:n;
   x(k) = vpa(vpa(x(k) - vpa(U(k,j)*x(j)))/U(k,k));
end

function x = bslashtxP(A,b)
% BSLASHTX  Solve linear system (backslash)
% x = bslashtx(A,b) solves A*x = b



[n,n] = size(A);
if isequal(triu(A,1),zeros(n,n))
   % Lower triangular
   x = forward(A,b);
   return
elseif isequal(tril(A,-1),zeros(n,n))
   % Upper triangular
   x = backsubs(A,b);
   return
elseif isequal(A,A')
   [R,fail] = chol(A);
   if ~fail
      % Positive definite
      y = forward(R',b);
      x = backsubs(R,y);
      return
   end
end

% Triangular factorization
[L,U,p] = lutxP(A);

% Permutation and forward elimination
y = forward(L,b(p));

% Back substitution
x = backsubs(U,y);


% ------------------------------

function x = forward(L,x)
% FORWARD. Forward elimination.
% For lower triangular L, x = forward(L,b) solves L*x = b.
[n,n] = size(L);
x(1) = x(1)/L(1,1);
for k = 2:n
   j = 1:k-1;
   x(k) = vpa(vpa(x(k) - vpa(L(k,j)*x(j)))/L(k,k));
end


% ------------------------------

function x = backsubs(U,x)
% BACKSUBS.  Back substitution.
% For upper triangular U, x = backsubs(U,b) solves U*x = b.
[n,n] = size(U);
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
   j = k+1:n;
   x(k) = vpa(vpa(x(k) - vpa(U(k,j)*x(j)))/U(k,k));
end

function [L,U,p,sig] = lutx(A)
%LUTX  Triangular factorization, textbook version
%   [L,U,p] = lutx(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U, and a permutation vector p,
%   so that L*U = A(p,:)



[n,n] = size(A);
p = (1:n)';
sig = 1;
for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   [r,m] = max(abs(A(k:n,k)));
   m = m+k-1;

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)

      % Swap pivot row
      if (m ~= k)
         A([k m],:) = A([m k],:);
         p([k m]) = p([m k]);
         sig=sig*-1;
      end

      % Compute multipliers
      i = k+1:n;
      A(i,k) = A(i,k)/A(k,k);

      % Update the remainder of the matrix
      j = k+1:n;
      A(i,j) = A(i,j) - A(i,k)*A(k,j);
   end
end

% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);

if nargout == 2
    for i=1:n
        q(i)=find(p==i);
    end
    L = L(q,:);
end

B9

Contents

Called Functions

A33

A34

A35

A36