function [P, L, U] = plu(A)

% The implementation of PLU Factorization
% L is lower triangular and U is upper triangular
% P is permutation matrix

% Author: Zhenlin Du(Johnsondu)
% Email:  qlduzhlin@126.com
% Time:   2014-11-27 22:00

A = double(A);
[m, n] = size(A);
L1 = zeros(m, n);
L = zeros(m, min(m, n));
U1 = zeros(n, n);
U = zeros(min(m, n), n);
P = eye(m);

% row operation
for i = 1: m
	mval = 0.0;
	row = i;
	
	% find maximum number in current column
    for k = i : min(i, n)
        for j = i: m
            if abs(mval) < abs(A(j, k))
                mval = A(j, k);
                row = j;
            end
        end
    end
	
	% if current maximum number is zero
	% process the next column
	if mval == 0
		continue;
    end
	
	% exchange process, in P, L, U
	if row ~= i
		tmp = A(i, :);
		A(i, :) = A(row, :);
		A(row, :) = tmp;
		tmp = P(i, :);
		P(i, :) = P(row, :);
		P(row, :) = tmp;
		tmp = L1(i, :);
		L1(i, :) = L1(row, :);
		L1(row, :) = tmp;
	end
	
	for j = i+1 : m
		ratio = A(j, i) / mval;
		A(j, :) = A(j, :) - ratio * A(i, :);
		L1(j, i) = ratio;
	end	
end

for i = 1: min(m, n)
	L1(i, i) = 1.0;
end

for i = 1: m
    for j = 1: min(m, n)
        L(i, j) = L1(i, j);
    end
end

U1 = A;
for i = 1: min(m, n)
    for j = 1 : n
        U(i, j) = U1(i, j);
    end
end
	

PLU decomposition Matlab version