/* \file MatlabMatrixToDoubleMatrix.cpp
 * \brief: convert a sparse or dense matlab matrix to a double* pointer
 */


#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif


#ifdef _SERIAL_

/*Petsc includes: */
#include "petscmat.h"
#include "petscvec.h"
#include "petscksp.h"

/*Matlab includes: */
#include "mex.h"

#include "../../../shared/shared.h"

int MatlabMatrixToDoubleMatrix(double** pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){

	int rows, cols;
	double* mxmatrix_ptr=NULL;
	int ierr;
	int i,j;

	/*output: */
	double* matrix=NULL;

	/*matlab indices: */
	mwIndex*    ir=NULL;
	mwIndex*    jc=NULL;
	double* pr=NULL;
	int     count;
	int     nnz;
	int     nz;

	/*Ok, first check if we are dealing with a sparse or full matrix: */
	if (mxIsSparse(mxmatrix)){

		/*Dealing with sparse matrix: recover size first: */
		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
		rows=mxGetM(mxmatrix);
		cols=mxGetN(mxmatrix);
		nnz=mxGetNzmax(mxmatrix);
		nz=(int)((double)nnz/(double)rows);

		matrix=(double*)xcalloc(rows*cols,sizeof(double));

		/*Now, get ir,jc and pr: */
		pr=mxGetPr(mxmatrix);
		ir=mxGetIr(mxmatrix);
		jc=mxGetJc(mxmatrix);

		/*Now, start inserting data into double* matrix: */
		count=0;
		for(i=0;i<cols;i++){
			for(j=0;j<(jc[i+1]-jc[i]);j++){
				*(matrix+rows*ir[count]+i)=pr[count];
				count++;
			}
		}

	}
	else{


		/*Dealing with dense matrix: recover pointer and size: */
		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
		rows=mxGetM(mxmatrix);
		cols=mxGetN(mxmatrix);

		
		/*Create serial matrix: */
		matrix=(double*)xcalloc(rows*cols,sizeof(double));

		for(i=0;i<rows;i++){
			for(j=0;j<cols;j++){
				*(matrix+cols*i+j)=*(mxmatrix_ptr+rows*j+i);
			}
		}
		
		#ifdef _DEBUG_
		for(i=0;i<rows;i++){
			for(j=0;j<cols;j++){
				printf("%g ",*(matrix+cols*i+j));
			}
			printf("\n");
		}
		#endif

	}

	/*Assign output pointer: */
	*pmatrix=matrix;
	*pmatrix_rows=rows;
	*pmatrix_cols=cols;

	return 1;
}
#endif
