function [A,P] = lupp(A);
%
% compute a LU factorization of A with partial pivoting
%
% calling sequence: [A,P] = lupp(A);
%
% input:
%   A = square matrix
% 
% output:
%   A contains L and U factors (original data overwritten)
%   P = permutation of [1,2,..n] indicating row permutation,
%       such that PA = LU
%