files = {
    'tumorAntiAngiogenesis_2.mtx', ...
    'olm500.mtx', ...
    'adder_dcop_05.mtx', ...
    'bayer10.mtx', ...
    'rajat01.mtx', ...
    'rajat19.mtx', ...
    'west0497.mtx', ...
    'gent113.mtx', ...
    'hangGlider_2.mtx', ...
    'Tina_AskCal.mtx', ...
    'GD98_a.mtx', ...
    'a0.mtx', ...
    'a1.mtx', ...
    'Tina_AskCal_perm.mtx', ...
    'ash219.mtx', ...
    'lp_e226.mtx', ...
    'r2.mtx', ...
    'LFAT5.mtx', ...
    'west0067.mtx', ...
    'arrow.mtx', ...
    'a2.mtx', ...
    'az88.mtx', ...
    'young1c.mtx', ...
    's32.mtx', ...
    'lp_share1b.mtx', ...
    'cage3.mtx', ...
    'b1_ss.mtx', ...
    'lfat5b.mtx', ...
    'c32.mtx', ...
    'bfwa62.mtx', ...
    'Ragusa16.mtx', ...
    'lp_e226_transposed.mtx', ...
    'row0.mtx', ...
    'Groebner_id2003_aug.mtx', ...
    'c2.mtx', ...
    'a4.mtx', ...
    'problem.mtx', ...
    'permuted_b1_ss.mtx', ...
    'pwr01b.mtx', ...
    'n3c4-b4.mtx', ...
    'Franz6_id1959_aug.mtx', ...
    'Ragusa16_pattern.mtx', ...
    'temp.mtx', ...
    'cage5.mtx', ...
    'lpi_itest6.mtx', ...
    'a04.mtx', ...
    'lpi_galenet.mtx', ...
    'nnc1374.mtx', ...
    '494_bus.mtx', ...
    'dwt_878.mtx', ...
    'xenon1.mtx', ...
    'c-62.mtx'} ;

warning ('off', 'MATLAB:singularMatrix') ;
warning ('off', 'MATLAB:nearlySingularMatrix') ;

for k = 1:length (files)
    file = files {k} ;
    fprintf ('\nmatrix: %s\n', file) ;
    % mread is in the CHOLMOD/MATLAB interface:
    A = mread (['../Matrix/' file]) ;
    [m n] = size (A) ;
    if (m ~= n)
        fprintf ('    rectangular\n') ;
        continue ;
    end
    b = (1:n)' ;
    anorm = norm (A, 1) ;
    if (anorm == 0)
        anorm = 1 ;
    end

    % backslash: uses iterative refinement if using LU factorization
    x = A\b ;
    xnorm = norm (x, 1) ;
    if (xnorm == 0)
        xnorm = 1 ;
    end
    rel_resid = norm (A*x-b,1) / (anorm * xnorm) ;
    fprintf ('    anorm %8.2e  rel_resid %8.2e ', anorm, rel_resid) ;

    % just use umfpack, but with row scaling.  No iterative refinement
    [L,U,P,Q,D] = lu (A) ;
    % P*inv(D)*A*Q = L*U
    x = Q * (U \ (L \ (P * (D\b)))) ;
    xnorm = norm (x, 1) ;
    if (xnorm == 0)
        xnorm = 1 ;
    end
    rel_resid = norm (A*x-b,1) / (anorm * xnorm);
    fprintf (' %8.2e', rel_resid) ;
    clear L U P Q D

    % umfpack, no row scaling.  No iterative refinement
    [L,U,P,Q] = lu (A) ;
    % P*A*Q = L*U
    x = Q * (U \ (L \ (P * (b)))) ;
    xnorm = norm (x, 1) ;
    if (xnorm == 0)
        xnorm = 1 ;
    end
    rel_resid = norm (A*x-b,1) / (anorm * xnorm) ;
    fprintf (' %8.2e\n', rel_resid) ;

    clear A Z x b L U P Q D Z

end