function [ P,L,U ] = PALUdecompV2( A )
%PALUdecompV2 decompose A = P^{-1}LU where P is stored as a vector
    [m, n]=size(A);
    if m ~= n, error('Input must be a square matrix.'), end
    L=zeros(n); U=zeros(n);
    P = 1:n;  % P starts as a vector with elements 1,2,... n

    % remember A^(0) is A
    AofK = A;  
    for k = 1:n
        
        % this is the new part, do the pivot...
        rk = findLargestinCol(AofK,k);
        if rk ~= k
            AofK = swapRows( AofK, rk, k );
            L = swapRows( L, rk, k );
            % swap two elements of P (instead of swapping rows)
            storeMe = P(rk); P(rk) = P(k); P(k) = storeMe;
        end
        
        % this is the OLD algorithm
        for j = k:n
            U(k,j) = AofK(k,j);
        end
        % check that we don't divide by zero(!)
        if U(k,k) == 0
            error('** A^(k-1)_{k,k}==0 in PALU decomp')
        end
        for i = k:n
           L(i,k) = AofK(i,k)/U(k,k);
        end
        for i = k:n
            for j = k:n
                AofK(i,j) = AofK(i,j) - L(i,k)*U(k,j);
            end
        end
    end  % of the loop over k

end % of the function

