function test_piro_band_svd()
% Usage : test_piro_band_svd()
% Finds the SVD using the band reduction algorithm for random input matrices and
% compares the result with MATLAB's svds.
% random input matrices.
%

%   Copyright 2009, Sivasankaran Rajamanickam, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

    fprintf('All norms printed should be small\n') ;
    fprintf('n1 is the norm(svds(A) - piro_band_svd(A)) \n') ;
    fprintf('n2 is the norm(A - (U * B * VT))  \n\n') ;

    A = banded(10, 10, 5, 5) ;
    piro_band_verify_svd(A, 0) ;

    A = banded(10, 20, 5, 5) ;
    piro_band_verify_svd(A, 1) ;
    piro_band_verify_svd(A, 0) ;

    A = banded(20, 10, 5, 5) ;
    piro_band_verify_svd(A, 1) ;
    piro_band_verify_svd(A, 0) ;

    A = banded(20, 10, 5, 6) ;
    piro_band_verify_svd(A, 1) ;
    piro_band_verify_svd(A, 0) ;

    fprintf('\n------------------------ Complex SVD ------------ \n\n') 
    A = banded(10, 10, 5, 5) ;
    A1 = banded(10, 10, 5, 5) ;
    A = A + i * A1 ;
    piro_band_verify_svd(A, 1) ;
    piro_band_verify_svd(A, 0) ;

    A = banded(10, 20, 6, 9) ;
    A1 = banded(10, 20, 6, 9) ;
    A = A + i * A1 ;
    piro_band_verify_svd(A, 1) ;
    piro_band_verify_svd(A, 0) ;

    fprintf('\n------------------------ Econ SVD ------------ \n\n') 
    A = banded(50, 60, 6, 9) ;
    piro_band_verify_svd_options(A) ;

    A = banded(60, 50, 6, 9) ;
    piro_band_verify_svd_options(A) ;

    A = banded(50, 50, 6, 9) ;
    piro_band_verify_svd_options(A) ;

    fprintf('\n------------------------ Complex Econ SVD ------------ \n\n') 
    A = banded(50, 60, 6, 9) ;
    A1 = banded(50, 60, 6, 9) ;
    A = A + i * A1 ;
    piro_band_verify_svd_options(A) ;

    A = banded(60, 50, 6, 9) ;
    A1 = banded(60, 50, 6, 9) ;
    A = A + i * A1 ;
    piro_band_verify_svd_options(A) ;

    A = banded(50, 50, 6, 9) ;
    A1 = banded(50, 50, 6, 9) ;
    A = A + i * A1 ;
    piro_band_verify_svd_options(A) ;

function piro_band_verify_svd(A, econ) 
%
    if (econ == 1)
        [U, S, VT] = piro_band_svd(A, 0);
    else
        [U, S, VT] = piro_band_svd(A) ;
    end
    A1 = U * S * VT ;
    %A - A1

    [m, n] = size(A) ;
    smat = svds(A, min(m,n)) ;
    snow = svd(S) ;
    n1 = norm(smat - snow) ;
    n2 = norm(A1 - A) ;

    fprintf('A -- %d x %d, --------- n1 = %g,    n2 = %g\n', m, n, n1, n2) ;

function piro_band_verify_svd_options(A)
    [m, n] = size(A) ;
    smat = svds(A, min(m,n)) ;

    % verify econ options 
    snow = piro_band_svd(A, 0);
    n1 = norm(smat - snow) ;
    fprintf('A -- %d x %d, --------- n1 = %g \n', m, n, n1) ;

    [U, S, VT] = piro_band_svd(A, 0);
    A1 = U * S * VT ;
    snow = svd(S) ;
    n1 = norm(smat - snow) ;
    n2 = norm(A1 - A) ;
    fprintf('A -- %d x %d, --------- n1 = %g,    n2 = %g\n', m, n, n1, n2) ;

    [U1, S ] = piro_band_svd(A, 0);
    snow = svd(S) ;
    n1 = norm(smat - snow) ;
    n3 = norm(U - U1) ;
    fprintf('A -- %d x %d, --------- n1 = %g,    n3 = %g\n', m, n, n1, n3) ;

    % verify full options.
    [U, S, VT] = piro_band_svd(A);
    A1 = U * S * VT ;
    snow = svd(S) ;
    n1 = norm(smat - snow) ;
    n2 = norm(A1 - A) ;
    fprintf('A -- %d x %d, --------- n1 = %g,    n2 = %g\n', m, n, n1, n2) ;

    snow = piro_band_svd(A);
    n1 = norm(smat - snow) ;
    fprintf('A -- %d x %d, --------- n1 = %g \n', m, n, n1) ;

    [U1, S] = piro_band_svd(A);
    snow = svd(S) ;
    n1 = norm(smat - snow) ;
    n3 = norm(U - U1) ;
    fprintf('A -- %d x %d, --------- n1 = %g,    n3 = %g\n', m, n, n1, n3) ;