/*=================================================================
* full2sparse.c
* This example demonstrates how to populate a sparse
* matrix.  For the purpose of this example, you must pass in a
* non-sparse 2-dimensional argument of type double.

* Comment: You might want to modify this MEX-file so that you can use
* it to read large sparse data sets into MATLAB.
*
* This is a MEX-file for MATLAB.  
* Copyright (c) 1984-98 by The MathWorks, Inc.
* All rights reserved.
*=================================================================*/

/* $Revision: 1.1 $ */

#include <math.h> /* Needed for the ceil() prototype */
#include "mex.h"

/* If you are using a compiler that equates NaN to be zero, you must
 * compile this example using the flag  -DNAN_EQUALS_ZERO. For example:
 *
 *     mex -DNAN_EQUALS_ZERO full2sparse.c
 *
 * This will correctly define the IsNonZero macro for your C compiler. */

#if NAN_EQUALS_ZERO
#define IsNonZero(d) ((d)!=0.0 || mxIsNaN(d))
#else
#define IsNonZero(d) ((d)!=0.0)
#endif

#define PERCENT_SPARSE .20

void mexFunction(
		 int nlhs,       mxArray *plhs[],
		 int nrhs, const mxArray *prhs[]
		 )
{
    /* Declare variable */
    int j,k,m,n,nzmax,*irs,*jcs,cmplx,isfull;
    double *pr,*pi,*si,*sr;
    
    /* Check for proper number of input and output arguments */    
    if (nrhs != 1) {
	mexErrMsgTxt("One input argument required.");
    } 
    if(nlhs > 1){
	mexErrMsgTxt("Too many output arguments.");
    }
    
    /* Check data type of input argument  */
    if (!(mxIsDouble(prhs[0]))){
	mexErrMsgTxt("Input argument must be of type double.");
    }	
    
    if (mxGetNumberOfDimensions(prhs[0]) != 2){
	mexErrMsgTxt("Input argument must be two dimensional\n");
    }
	       
    /* Get the size and pointers to input data */
    m  = mxGetM(prhs[0]);
    n  = mxGetN(prhs[0]);
    pr = mxGetPr(prhs[0]);
    pi = mxGetPi(prhs[0]);
    cmplx = (pi==NULL ? 0 : 1);
    
    /* Allocate space for sparse matrix */
    /* NOTE:  Assume at most 20% of the data is sparse.  Use ceil
     * to cause it to round up. */
    nzmax=(int)ceil((double)m*(double)n*PERCENT_SPARSE);
    /* NOTE: The maximum number of non-zero elements cannot be less
       than the number of columns in the matrix. */
    if (n>nzmax){
	nzmax=n;
    }
    plhs[0] = mxCreateSparse(m,n,nzmax,cmplx);
    sr  = mxGetPr(plhs[0]);
    si  = mxGetPi(plhs[0]);
    irs = mxGetIr(plhs[0]);
    jcs = mxGetJc(plhs[0]);
    
    /* Copy nonzeros */
    k = 0; 
    isfull=0;
    for (j=0; (j<n ); j++) {
	int i;
	jcs[j] = k;
	for (i=0; (i<m ); i++) {
	    if (IsNonZero(pr[i]) || (cmplx && IsNonZero(pi[i]))) {
		/* Check to see if non-zero element will fit in 
		 * allocated output array.  If not, set flag and 
		 * print warning. */
		if (k>=nzmax){
		    isfull=1;
		    mexWarnMsgTxt("Truncating output, more than 20%% of input \
is non-zero data."); 
		    break;
		}
		sr[k] = pr[i];
		if (cmplx){
		    si[k]=pi[i];
		}
		irs[k] = i;
		k++;
	    }
	}
	if (isfull)
	  break;	
	pr += m;
	pi += m;
    }
    jcs[n] = k;
}