% glu.m: general LU factorization without pivoting. 

  clear;
  format short;

% Step 0: Assign the matrix A and initialize two vectors d and w. 
%         The components of the vector d are prescribed (nonzero)
%         values of the diagonal entries. The vector w is a logical
%         array: if w(k) = 1, then L(k,k) = d(k); if w(k) = -1, then
%         U(k,k) = d(k); if w(k) = 0, then L(k,k) = U(k,k). The last
%         case corresponds to the Cholesky factorization.  
 
  n = 4;
  A = zeros(n,n);
  for i = 1:n
    for j = 1:n
      A(i,j) = 1/(i+j-1);
    end 
  end 
  
  d = zeros(n,1);
  for i = 1:n
    d(i) = 1;
  end  

  w = zeros(n,1);
  for i = 1:n
    w(i) = 1;          % this gives the Doolittle factorization
%   w(i) = -1;         % this gives the Crout factorization
%   w(i) = 0;          % this gives the Cholesky factorization
  end 

% Step 1: Basic Gaussian Elimination.

  L = zeros(n,n);      % initialize the matrix L
  U = zeros(n,n);      % initialize the matrix U
  for k = 1:n
    if w(k) == 1
      L(k,k) = d(k);
      U(k,k) = A(k,k);
      for s = 1:(k-1)
        U(k,k) = U(k,k)-L(k,s)*U(s,k);
      end 
      U(k,k) = U(k,k)/L(k,k);
      for j = (k+1):n
        U(k,j) = A(k,j);
        for s = 1:(k-1)
          U(k,j) = U(k,j)-L(k,s)*U(s,j);
        end
        U(k,j) = U(k,j)/L(k,k);
      end
      for i = (k+1):n
        L(i,k) = A(i,k);
        for s = 1:(k-1)
          L(i,k) = L(i,k)-L(i,s)*U(s,k);
        end
        L(i,k) = L(i,k)/U(k,k);
      end
    elseif w(k) == -1
%     (fill out this part for the Crout factorization)
    else
%     (fill out this part for the Cholesky factorization) 
    end
  end
  L                    % output L
  U                    % output U