/* ============================ storebandmex.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
 * -------------------------------------------------------------------------- */

/* Mex interface to store a sparse/full matrix in packed band storage.
 * Usage:
 *      B = storeband(A) ;
 *      where A is a band matrix of size m x n with upper bandwidth bu and lower
 *      bandwidth bl. B is a packed band matrix of size (bl + bu + 1) x n.
 *      */

#include "piro_band_matlab.h"


void mexFunction
(
    int nlhs,
    mxArray *plhs [],
    int nrhs,
    const mxArray *prhs []
)
{
    Int ub, lb ;	
    Int m, n ;	
    Int *Ap, *Ai ;
    double *Ax, *Axi  ;
    double *Cx, *Cxi ;
    double *rCx ;
    Int crow ;
    Int i ;
    Int msize ;
    bool iscomplex ;
 
    if (nlhs != 1 || nrhs != 1)
    {
        mexErrMsgTxt("Incorrect no of arguments to the mex function\n") ;
    }    

    Ax = mxGetPr(prhs[0]) ;
    iscomplex = mxIsComplex(prhs[0]) ;
    Axi = NULL ;
    if (iscomplex)
    {
	Axi = mxGetPi(prhs[0]) ;
    }

    m = mxGetM(prhs[0]) ;
    n = mxGetN(prhs[0]) ;

    /* Estimate the bandwidth and store the matrix in packed band format */
    if (mxIsSparse(prhs[0]))
    {
	Ap = (Int *) mxGetJc(prhs[0]) ;
	Ai = (Int *) mxGetIr(prhs[0]) ;
        PIRO_BAND(find_bandwidth)(m, n, Ap, Ai, &lb, &ub) ;

        crow = ub+lb+1 ;
        msize = iscomplex ? 2 * crow * n : crow * n ;
        rCx = (double *) mxMalloc(msize * sizeof(double)) ;
	PIRO_BAND(storeband_withzeroes)( Ap, Ax, Axi, m, n, rCx, ub, lb) ;
    }
    else
    {
        lb = 0 ;
        ub = 0 ;
        PIRO_BAND(find_full_bandwidth)(m, n, Ax, &lb, &ub) ;
        if (iscomplex)
        {
            PIRO_BAND(find_full_bandwidth)(m, n, Axi, &lb, &ub) ;
        }

        crow = ub+lb+1 ;
        msize = iscomplex ? 2 * crow * n : crow * n ;
        rCx = (double *) mxMalloc(msize * sizeof(double)) ;
	PIRO_BAND(storeband_withzeroes_full)( Ax, Axi, m, n, rCx, ub, lb) ;
    }

    /* Copy the result to the MATLAB data structures */
    Cxi = NULL ;
    if (iscomplex)
    {
	plhs[0] = mxCreateDoubleMatrix (crow, n, mxCOMPLEX) ;
	Cxi = mxGetPi(plhs[0]) ;
    }
    else
    {
	plhs[0] = mxCreateDoubleMatrix (crow, n, mxREAL) ;
    }
    Cx = mxGetPr(plhs[0]) ;

    for (i = 0 ; i < crow * n ; i++)
    {
	if (iscomplex)
	{
	    Cx[i] = rCx[2*i] ;
	    Cxi[i] = rCx[2*i+1] ;
	}
	else
	{
	    Cx[i] = rCx[i] ;
	}
    }

    /* Free work space */
    mxFree(rCx) ;

}