PLU decomposition Matlab version

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:  [e-mail]% Time:   2014-11-27 22:00A = 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 operationfor i = 1: mmval = 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 columnif mval == 0continue;    end% exchange process, in P, L, Uif row ~= itmp = 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;endfor j = i+1 : mratio = A(j, i) / mval;A(j, :) = A(j, :) - ratio * A(i, :);L1(j, i) = ratio;endendfor i = 1: min(m, n)L1(i, i) = 1.0;endfor i = 1: m    for j = 1: min(m, n)        L(i, j) = L1(i, j);    endendU1 = A;for i = 1: min(m, n)    for j = 1 : n        U(i, j) = U1(i, j);    endend