/* ============================ test_piro_band.c ============================ */

/* -----------------------------------------------------------------------------
 * PIRO_BAND.  Version 0.9.
 * Copyright (C) 2009, Sivasankaran Rajamanickam, Timothy A. Davis.
 * PIRO_BAND is licensed under Version 2.0 of the GNU
 * General Public License.  See gpl.txt for a text of the license.
 * PIRO_BAND is also available under other licenses; contact authors for 
 * details. http://www.cise.ufl.edu/research/sparse
 * -------------------------------------------------------------------------- */

/* Simple test for band reduction of a small unsymmetric matrix. 
 * */

#include "piro_band.h"

/* 
 * Test the band reduction routine using a matrix stored in a file.
 * Input:
 *  sizefile - Size file name. The size file will have the dimensions of the
 *             the band matrix, the lower and upper bandwidths.
 *  pfile    - Matrix file name. The file with the numerical entries 
 *             (real/complex) of the band matrix stored in packed band format.
 * */
static int PIRO_BAND(test)(
    char *sizefile, 
    char *pfile
)
{
    int m, n, bl, bu ;
    int ldab  ;
    double *A , *work, *U, *V, *C, *Atemp ;
    double *D, *E, *diag ;
    double *Aptr, *AllA ;
    double *mtemp ;
    int blk[4] ;
    int dsize ;
    int err ;

    /* Get the size of the matrix from the file */
    if (!piro_band_get_matrix_size_dri(sizefile, &m, &n, &bl, &bu)) 
    {
	printf("Error in reading %s file\n", sizefile) ;
	return 1 ;
    }
    printf("m = "ID", n= "ID", bl="ID", bu="ID"\n", m, n, bl, bu) ;

    ldab = bl+bu+1 ;

    /* Allocate memory two times the required memory for the matrix and 
     * keep a copy as band reduction will destroy the input matrix. */
    AllA = (double *) malloc( 2 * n * (ldab) * sizeof(double)) ;
    if (!AllA)
    {
	printf("Out of memory") ;
	return 1 ;
    }

    Aptr = AllA ;
    A = Aptr ;
    Aptr += (n * ldab) ;
    Atemp = Aptr ;

    /* Get the numerical values from the file */
    if (!piro_band_get_matrix_dri(pfile, ldab, n, A))
    {
	printf("Error in reading %s file\n", pfile) ;
	free(AllA) ;
	return 1 ;
    }

    memcpy(Atemp, A, n * ldab * sizeof(double)) ; 

    /* Allocate space for the singular vectors */
    if (!piro_band_init_input_dri(m, n, &mtemp, &U, &V, &C)) 
    {
        printf("Error in initializing input \n") ;
        free(AllA) ;
        return 1 ;
    }

    /* Allocate memory bidiagonal (D, E). */
    diag = (double *) malloc( 2 * MIN(m, n) * sizeof(double)) ;
    if (!diag)
    {
	printf("Out of memory") ;
        free(AllA) ;
        free(mtemp) ;
    }
    D = diag ;
    E = diag + MIN(m, n) ;
    E[0] = 0.0 ;

    /* Get the suggested block size for the given matrix and allocate 
     * memory for the workspace.
     * */
    piro_band_get_blocksize(m, n, bl, bu, 1, blk) ;

    dsize = blk[0]*blk[1] > blk[2]*blk[3] ? blk[0]*blk[1] : blk[2]*blk[3] ;
    work = (double *) malloc(2 * dsize * sizeof(double)) ;
    if (!work)
    {
        free(AllA) ;
        free(mtemp) ;
        free(diag) ;
        printf("Out of memory") ;
        return 1 ;
    }

    /* Reduce the input matrix to the bidiagonal form */
    err = piro_band_reduce_dri(blk, m, n, m, bl, bu, A, ldab, D, E+1, 
            U, m, V, n, C, m, work, 0) ;
    if (err != 0)
    {
        printf("test Band reduction failed "ID"\n", err) ;
    }

    /* Check that the A - (U * S * V') is small enough */
    if (!piro_band_chk_output_dri(m, n, bu, bl, ldab, 0, D, E, U, V, C, 
                                Atemp, 0))
    {
	free(AllA) ;
	free(mtemp) ;
        free(diag) ;
	if (!work) free(work) ;
	printf("Chk_output failed\n") ;
	return 1 ;
    }

    /* Free Allocated space */
    free(AllA) ;
    free(mtemp) ;
    free(diag) ;
    if (!work) free(work) ;
    return 0 ;
}


/* Run demo for PIRO_BAND */
int main(void)
{
    char *fname ;
    int blks[4] ;

    fname = (char *) malloc (32 * sizeof(char)) ;
    if (!fname)
    {
        printf("Out of memory\n") ;
        return ;
    }

    COPY_FNAME(fname, "A2") ;
    PIRO_BAND(test)("A2size", fname) ;

    free(fname) ;
    return 0 ;
}