function Feb5()

n = 5;
A = randn(n,n);
C = randn(n,n);
b = randn(n,1);
d = randn(n,1);
% Method 1..
finalSol = [A C ; zeros(n,n) A']\[b ; d];
% Method 2...
[L,U] = GE(A);
% Solve A'y = d for y.   U'L'y = d
z = LTriSol(U',d);
y = UTriSOl(L',z);
% Solve Ax = b - C*y for x...
w = LTriSol(L,b-C*y);
x = UTriSol(U,w);
final = [x;y]
[finalSol final]



  function [L,U] = GE(A)
% [L,U] = GE(A)
%
% The LU factorization without pivoting. If A is n-by-n and has an
% LU factorization, then L is unit lower triangular and U is upper
% triangular so A = LU. 

[n,n] = size(A); 
for k=1:n-1
   A(k+1:n,k) = A(k+1:n,k)/A(k,k);
   A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end
L = eye(n,n) + tril(A,-1);
U = triu(A);

  function x = UTriSol(U,b)
% x = UTriSol(U,b)
%
% Solves the nonsingular upper triangular system  Ux = b.
% where U is n-by-n, b is n-by-1, and X is n-by-1.

n = length(b);
x = zeros(n,1);
for j=n:-1:2
   x(j) = b(j)/U(j,j);
   b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = b(1)/U(1,1);

  function x = LTriSol(L,b)
% x = LTriSol(L,b)
%
% Solves the nonsingular lower triangular system  Lx = b 
% where L is n-by-n, b is n-by-1, and x is n-by-1.

n = length(b);
x = zeros(n,1);
for j=1:n-1
   x(j) = b(j)/L(j,j);
   b(j+1:n) = b(j+1:n) - L(j+1:n,j)*x(j);
end
x(n) = b(n)/L(n,n);

