function [A,p] = ludec(A,pivotrule,show)
% ludec: performs LU decomposition with or without pivoting on a matrix
%
% usage: [A,p] = ludec(A,pivotrule,show)
%
% inputs:        A: the matrix (should be square)
%        pivotrule: (optional) the type of pivoting required
%             show: unless this is 0 then the program displays the matrix at
%                   each stage of the decomposition

if nargin<3,show=1;end

if nargin<2,
  disp('Choose your pivoting strategy');
  disp('1) No pivoting');
  disp('2) Non-zero pivoting');
  disp('3) Partial pivoting');
  disp('4) Scaled partial pivoting');
  pivotrule = input('Please type in your selection: ');
  if ~any(pivotrule==[1 2 3 4])
    disp(['Your choice ',int2str(pivotrule),' was not recognised.']);
    pivotrule = input('Please type in a different selection: ');
  end
end

if(diff(size(A))),error('Your matrix A is not square.'); end

n = max(size(A));

p = 1:n;

switch pivotrule
  case 1 % No pivoting
    for k = 1:n-1
      for i = k+1:n
        mult = A(i,k)/A(k,k);
        A(i,k) = mult;
        for j = k+1:n
          A(i,j) = A(i,j) - mult*A(k,j);
        end
      end
      if(show)
        disp(['At step ',int2str(k),' the matrix is']);
        disp(A);
      end
    end
  case 2 % Non-zero pivoting only
    for k = 1:n-1
      j = min(find(A(p(k:n),k))); % look for first non-zero element in the kth
                                % column of A beneath the diagonal
      j = j+k-1;
      p([k,j]) = p([j,k]);    % swap rows k and j in p
      for i = k+1:n
        mult = A(p(i),k)/A(p(k),k);
        A(p(i),k) = mult;
        for j = k+1:n
          A(p(i),j) = A(p(i),j) - mult*A(p(k),j);
        end
      end
      if(show)
        disp(['At step ',int2str(k),' the matrix is']);
        disp(A);
        disp(['and the permutation vector is']);
        disp(['[',int2str(p(:)'),']']);
      end
    end
  case 3 % Partial pivoting
    for k = 1:n-1
      [tmp,j] = max(abs(A(p(k:n),k))); 
      % Look for the largest element in the kth
      % column of A beneath the diagonal.
      % Matlab automatically picks out the first.
      j = j+k-1;
      p([k,j]) = p([j,k]);     % Swap rows k and j in p.
      for i = k+1:n
        mult = A(p(i),k)/A(p(k),k);
        A(p(i),k) = mult;
        for j = k+1:n
          A(p(i),j) = A(p(i),j) - mult*A(p(k),j);
        end
      end
      if(show)
        disp(['At step ',int2str(k),' the matrix is']);
        disp(A);
        disp(['and the permutation vector is']);
        disp(['[',int2str(p(:)'),']']);
      end
    end
  case 4 % Scaled partial pivoting
    s = max(abs(A')); % a vector of maximum row elements.
    for k = 1:n-1
      [tmp,j] = max(abs(A(p(k:n),k)./s(p(k:n))'));
      j = j+k-1;
      % Look for the largest scaled element in the kth column of A beneath
      % the diagonal.
      % Matlab automatically picks out the first of these.
      p([k,j]) = p([j,k]);  % Swap rows k and j in p.
      for i = k+1:n
        mult = A(p(i),k)/A(p(k),k);
        A(p(i),k) = mult;
        for j = k+1:n
          A(p(i),j) = A(p(i),j) - mult*A(p(k),j);
        end
      end
      if(show)
        disp(['At step ',int2str(k),' the matrix is']);
        disp(A);
        disp(['and the permutation vector is']);
        disp(['[',int2str(p(:)'),']']);
      end
    end
end