function test_piro_band_special()
% 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') ;

    fprintf('\n--------------Unsymmetric case  ---------------- \n\n') ;
    % diagonal case not allowed in this interface.

    % bidiagonal case
    A = banded(10, 10, 1, 0) ;
    piro_band_verify(A, 0, 0, 0, 0) ;
    A = banded(10, 20, 1, 0) ;
    piro_band_verify(A, 0, 0, 0, 0) ;
    A = banded(20, 10, 1, 0) ;
    piro_band_verify(A, 0, 0, 0, 0) ;

    % tridiagonal case
    A = banded(10, 10, 1, 1) ;
    piro_band_verify(A, 0, 0, 1, 1) ;
    A = banded(10, 20, 1, 1) ;
    piro_band_verify(A, 0, 0, 1, 1) ;
    A = banded(20, 10, 1, 1) ;
    piro_band_verify(A, 0, 0, 1, 1) ;

    % upper banded case
    A = banded(10, 10, 5, 0) ;
    piro_band_verify(A, 2, 2, 0, 0) ;
    A = banded(10, 20, 6, 0) ;
    piro_band_verify(A, 3, 3, 0, 0) ;
    A = banded(20, 10, 5, 0) ;
    piro_band_verify(A, 4, 1, 0, 0) ;

    % lower banded case
    A = banded(10, 10, 1, 5) ;
    piro_band_verify(A, 0, 0, 1, 5) ;
    A = banded(10, 20, 1, 7) ;
    piro_band_verify(A, 0, 0, 7, 1) ;
    A = banded(20, 10, 1, 6) ;
    piro_band_verify(A, 0, 0, 3, 3) ;

    fprintf('\n--------------Symmetric case  ---------------- \n\n') ;
    % Only the upper half is stored.

    % tridiagonal case
    A = banded(10, 10, 1, 0) ;
    piro_band_sym_verify(A, 0, 0, 0, 0) ;

    A = banded(10, 10, 2, 0) ;
    piro_band_sym_verify(A, 1, 1, 0, 0) ;
    A = banded(10, 10, 5, 0) ;
    piro_band_sym_verify(A, 2, 2, 0, 0) ;
    A = banded(10, 10, 5, 0) ;
    piro_band_sym_verify(A, 2, 3, 0, 0) ;

    fprintf('\n-------------- Complex interface ---------------- ') ;
    fprintf('\n--------------Unsymmetric case  ---------------- \n\n') ;
    % diagonal case not allowed in this interface.

    % bidiagonal case
    A = banded(10, 10, 1, 0) ;
    A1 = banded(10, 10, 1, 0) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 0, 0) ;
    A = banded(10, 20, 1, 0) ;
    A1 = banded(10, 20, 1, 0) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 0, 0) ;
    A = banded(20, 10, 1, 0) ;
    A1 = banded(20, 10, 1, 0) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 0, 0) ;

    % tridiagonal case
    A = banded(10, 10, 1, 1) ;
    A1 = banded(10, 10, 1, 1) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 1, 1) ;
    A = banded(10, 20, 1, 1) ;
    A1 = banded(10, 20, 1, 1) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 1, 1) ;
    A = banded(20, 10, 1, 1) ;
    A1 = banded(20, 10, 1, 1) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 1, 1) ;

    % upper banded case
    A = banded(10, 10, 5, 0) ;
    A1 = banded(10, 10, 5, 0) ;
    A = A + i * A1 ;
    piro_band_verify(A, 2, 2, 0, 0) ;
    A = banded(10, 20, 6, 0) ;
    A1 = banded(10, 20, 6, 0) ;
    A = A + i * A1 ;
    piro_band_verify(A, 3, 3, 0, 0) ;
    A = banded(20, 10, 5, 0) ;
    A1 = banded(20, 10, 5, 0) ;
    A = A + i * A1 ;
    piro_band_verify(A, 4, 1, 0, 0) ;

    % lower banded case
    A = banded(10, 10, 1, 5) ;
    A1 = banded(10, 10, 1, 5) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 1, 5) ;
    A = banded(10, 20, 1, 7) ;
    A1 = banded(10, 20, 1, 7) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 7, 1) ;
    A = banded(20, 10, 1, 6) ;
    A1 = banded(20, 10, 1, 6) ;
    A = A + i * A1 ;
    piro_band_verify(A, 0, 0, 3, 3) ;

    fprintf('\n--------------Symmetric case  ---------------- \n\n') ;
    % Only the upper half is stored.

    % tridiagonal case
    A = banded(10, 10, 1, 0) ;
    A1 = banded(10, 10, 1, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_sym_verify(A, 0, 0, 0, 0) ;

    A = banded(10, 10, 2, 0) ;
    A1 = banded(10, 10, 2, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_sym_verify(A, 1, 1, 0, 0) ;
    A = banded(10, 10, 5, 0) ;
    A1 = banded(10, 10, 5, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_sym_verify(A, 2, 2, 0, 0) ;
    A = banded(10, 10, 5, 0) ;
    A1 = banded(10, 10, 5, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_sym_verify(A, 2, 3, 0, 0) ;


    fprintf('\n-------------- LAPACK style interface ---------------- \n\n') ;
    A = banded(10, 10, 0, 0) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 0, 0) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 0, 0) ;
    piro_band_lapack_verify(A) ;

    % bidiagonal case
    A = banded(10, 10, 1, 0) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 1, 0) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 1, 0) ;
    piro_band_lapack_verify(A) ;

    % tridiagonal case
    A = banded(10, 10, 1, 1) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 1, 1) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 1, 1) ;
    piro_band_lapack_verify(A) ;

    % upper banded case
    A = banded(10, 10, 5, 0) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 6, 0) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 5, 0) ;
    piro_band_lapack_verify(A) ;

    % lower banded case
    A = banded(10, 10, 0, 5) ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 0, 7) ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 0, 6) ;
    piro_band_lapack_verify(A) ;

    fprintf('\n--------------Symmetric case  ---------------- \n\n') ;
    % Only the upper half is stored.
    A = banded(10, 10, 0, 0) ;
    piro_band_lapack_sym_verify(A) ;

    % tridiagonal case
    A = banded(10, 10, 1, 0) ;
    piro_band_lapack_sym_verify(A) ;

    A = banded(10, 10, 2, 0) ;
    piro_band_lapack_sym_verify(A) ;
    A = banded(10, 10, 5, 0) ;
    piro_band_lapack_sym_verify(A) ;
    A = banded(10, 10, 6, 0) ;
    piro_band_lapack_sym_verify(A) ;

    % Only the lower half is stored.
    % tridiagonal case
    A = banded(10, 10, 0, 1) ;
    piro_band_lapack_sym_verify(A, 1) ;

    A = banded(10, 10, 0, 2) ;
    piro_band_lapack_sym_verify(A, 1) ;
    A = banded(10, 10, 0, 5) ;
    piro_band_lapack_sym_verify(A, 1) ;
    A = banded(10, 10, 0, 6) ;
    piro_band_lapack_sym_verify(A, 1) ;

    fprintf('\n-------------- Complex interface ---------------- \n\n') ;
    % diagonal case not allowed in this interface.
    A = banded(10, 10, 0, 0) ;
    A1 = banded(10, 10, 0, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 0, 0) ;
    A1 = banded(10, 20, 0, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 0, 0) ;
    A1 = banded(20, 10, 0, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

    % bidiagonal case
    A = banded(10, 10, 1, 0) ;
    A1 = banded(10, 10, 1, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 1, 0) ;
    A1 = banded(10, 20, 1, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 1, 0) ;
    A1 = banded(20, 10, 1, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

    % tridiagonal case
    A = banded(10, 10, 1, 1) ;
    A1 = banded(10, 10, 1, 1) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 1, 1) ;
    A1 = banded(10, 20, 1, 1) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 1, 1) ;
    A1 = banded(20, 10, 1, 1) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

    % upper banded case
    A = banded(10, 10, 5, 0) ;
    A1 = banded(10, 10, 5, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 6, 0) ;
    A1 = banded(10, 20, 6, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 5, 0) ;
    A1 = banded(20, 10, 5, 0) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

    % lower banded case
    A = banded(10, 10, 0, 5) ;
    A1 = banded(10, 10, 0, 5) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(10, 20, 0, 7) ;
    A1 = banded(10, 20, 0, 7) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;
    A = banded(20, 10, 0, 6) ;
    A1 = banded(20, 10, 0, 6) ;
    A = A + i * A1 ;
    piro_band_lapack_verify(A) ;

    fprintf('\n--------------Symmetric case  ---------------- \n\n') ;
    % Only the upper half is stored.
    A = banded(10, 10, 0, 0) ;
    A1 = banded(10, 10, 0, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A) ;

    % tridiagonal case
    A = banded(10, 10, 1, 0) ;
    A1 = banded(10, 10, 1, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A) ;

    A = banded(10, 10, 2, 0) ;
    A1 = banded(10, 10, 2, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A) ;
    A = banded(10, 10, 5, 0) ;
    A1 = banded(10, 10, 5, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A) ;
    A = banded(10, 10, 5, 0) ;
    A1 = banded(10, 10, 5, 0) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A) ;

    % Only the lower half is stored.
    % tridiagonal case
    A = banded(10, 10, 0, 1) ;
    A1 = banded(10, 10, 0, 1) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A, 1) ;

    A = banded(10, 10, 0, 2) ;
    A1 = banded(10, 10, 0, 2) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A, 1) ;
    A = banded(10, 10, 0, 5) ;
    A1 = banded(10, 10, 0, 5) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A, 1) ;
    A = banded(10, 10, 0, 6) ;
    A1 = banded(10, 10, 0, 6) ;
    d = diag(A1) ;
    A1 = A1 - diag(d) ;
    A = A + i * A1 ;
    piro_band_lapack_sym_verify(A, 1) ;


% 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 for symmetric case and verify the result.
function piro_band_sym_verify(A, nc, nr, ncl, nrl) 
    [m, n] = size(A) ;
    d = diag(A) ;
    A2 = A + A' - diag(d) ;
    sgood = svds(A2, n) ;
    blk = [ nc, nr, ncl, nrl] ;
    [b1, b2, U] = piro_band(A, 1, blk) ;

    Res = sparse(zeros(m, n)) ;
    Res = spdiags([[b2(2:n) ; 0] b1, b2] , -1:1, Res) ;

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

    A1 = U * Res * U' ;
    n2 = norm(A1 - A2) ;
    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) ;


% Find the symmetric SVD using LAPACK style call and verify the result.
function piro_band_lapack_sym_verify(A, uplo) 
    if (nargin < 2)
        uplo = 0 ;
    end

    [m, n] = size(A) ;
    d = diag(A) ;
    A2 = A + A' - diag(d) ;
    sgood = svds(A2, n) ;
    if (uplo == 0)
        [b1, b2, U] = piro_band_lapack(A, 1) ;
    else
        [b1, b2, U] = piro_band_lapack(A, 1, 3.14) ;
    end

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

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

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