function test_piro_band()
% Usage : test_piro_band()
% Tests the band reduction method with 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) - svds(B)) \n') ;
    fprintf('n2 is the norm(A - (U * B * VT))  \n\n') ;
    A = banded(10, 10, 5, 5) ;
    piro_band_verify(A, 2, 2, 2, 2) ;
    A = banded(10, 10, 6, 5) ;
    piro_band_verify(A, 3, 2, 2, 2) ;
    A = banded(100, 100, 30, 30) ;
    piro_band_verify(A, 10, 10, 10, 10) ;
    A = banded(100, 100, 30, 40) ;
    piro_band_verify(A, 15, 15, 17, 9) ;
    A = banded(10, 20, 6, 9) ;
    piro_band_verify(A, 3, 2, 3, 2) ;
    A = banded(20, 10, 7, 8) ;
    piro_band_verify(A, 3, 3, 3, 3) ;
    A = banded(20, 30, 11, 11) ;
    piro_band_verify(A, 3, 5, 3, 5) ;
    A = banded(30, 20, 11, 11) ;
    piro_band_verify(A, 3, 5, 3, 5) ;

    fprintf('\n-------------- LAPACK style interface ---------------- \n\n') ;
    A = banded(10, 10, 5, 5) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 10, 6, 5) ;
    piro_band_lapack_verify(A) ;
    A = banded(100, 100, 30, 30) ;
    piro_band_lapack_verify(A) ;
    A = banded(100, 100, 30, 40) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 6, 9) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 7, 8) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 30, 11, 11) ;
    piro_band_lapack_verify(A) ;
    A = banded(30, 20, 11, 11) ;
    piro_band_lapack_verify(A) ;

    fprintf('\n-------------- Complex interface ---------------- \n\n') ;
    A = banded(10, 10, 5, 5) ;
    A1 = banded(10, 10, 5, 5) ;
    A = A + i * A1 ;
    piro_band_verify(A, 2, 2, 2, 2) ;

    A = banded(10, 20, 6, 9) ;
    A1 = banded(10, 20, 6, 9) ;
    A = A + i * A1 ;
    piro_band_verify(A, 3, 2, 3, 2) ;

    fprintf('\n-------------- Complex interface ---------------- \n\n') ;
    A = banded(10, 10, 5, 5) ;
    A1 = banded(10, 10, 5, 5) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

    A = banded(10, 20, 6, 9) ;
    A1 = banded(10, 20, 6, 9) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

% Find the SVD and verify the result.
function piro_band_verify(A, nc, nr, ncl, nrl) 
    [m, n] = size(A) ;
    sgood = svds(A, n) ;
    blk = [ nc, nr, ncl, nrl] ;
    [b1, b2, U, V] = piro_band(A, 0, blk) ;

    Res = sparse(zeros(m, n)) ;
    if (m < n)
      Res = spdiags([b1, [b2(2:m) ; 0]] , 0:1, Res) ;
    else
      Res = spdiags([b1, b2] , 0:1, Res) ;
    end

    s = svds(Res, n) ;
    n1 = norm(sgood-s) ;

    A1 = U * Res * V' ;
    n2 = norm(A1 - A) ;
    fprintf('A -- %d x %d, --------- n1 = %g,      n2 = %g \n',  m, n, n1, n2) ;

% Find the SVD using LAPACK style call and verify the result.
function piro_band_lapack_verify(A) 
    [m, n] = size(A) ;
    sgood = svds(A, n) ;
    [b1, b2, U, VT] = piro_band_lapack(A) ;

    Res = sparse(zeros(m, n)) ;
    if (m < n)
      Res = spdiags([b1, b2] , 0:1, Res) ;
    else
      Res = spdiags([b1, [0 ; b2(1:n-1) ]] , 0:1, Res) ;
    end

    s = svds(Res, n) ;
    n1 = norm(sgood-s) ;

    A1 = U * Res * VT ;
    n2 = norm(A1 - A) ;
    fprintf('A -- %d x %d, --------- n1 = %g,    n2 = %g\n', m, n, n1, n2) ;