function [A,B,pi]=testLU( m, n, p, dbg )
    % Function assumes BSPlu has been run on an random matrix with
    % m by n threads.
    % example: testLU(2,2)       (for original BSPlu)
    %          testLU(2,2,1)     (for BSPlu_v2-- MulticoreBSP variant)
    %          testLU(2,2,0/1,1) (for debug mode)
    %outputs [original,lu,pi]
    if( nargin < 3 )
        p=m;
        A=readCyclic( m, n, 'lu_input' );
        B=readCyclic( m, n, 'lu_output' );
    else
        p=m*n;
        A=readCyclic( m, n, 'lu_input', 1 );
        B=readCyclic( m, n, 'lu_output', 1 );
    end
    if( mod(size(A,1),p)~=0 )
        pi=readCyclic( p, 1, 'pi', 0, size(A,1) ); %specify vector length
    else
        pi=readCyclic( p, 1, 'pi' ); %specification not necessary
    end
    if( nargin > 3 )
        A
        pi
    end

    pi=pi+1;
    LU=lu(A);
    C=B-LU;
    disp(sprintf( 'Compared with MATLAB Y=lu(A):      %g', norm(C) ));
    disp(sprintf( 'Maximum absolute difference:       %g', max(max(abs(C))) ));
    [L,U]=extractLU(B);
    D=A(pi,:)-L*U;
    disp(sprintf( 'Norm of recalculation A(pi,:)-L*U: %g', norm(D) ));
    disp(sprintf( 'Maximum absolute difference:       %g', max(max(abs(D))) ));
end