/* \file PetscMatrixToMatlabMatrix.cpp
 * \brief: convert a sparse or dense Petsc matrix into a matlab matrix
 */


#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"

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

#include "../../../shared/shared.h"
#include <string>


int PetscMatrixToMatlabMatrix(mxArray** pdataref,Mat matrix){

	int i,j,k;

	/*output: */
	mxArray* dataref=NULL;
	mxArray* tdataref=NULL;
	int    rows,cols;

	int ncols;
	int nnz=0;
	int nzmax=0;
	const int* columns=NULL;
	const double* column_values=NULL;
	double* matrix_ptr=NULL;

	/*compress row format: */
	double* val=NULL;
	mwIndex*    col_ind=NULL;
	mwIndex*    row_ptr=NULL;

	/*petsc type: */
	MatType type;
	int*    idxm=NULL;
	int*    idxn=NULL;

	MatGetType(matrix,&type);
	MatGetSize(matrix,&rows,&cols);

	if (strcmp(type,MATSEQAIJ)==0){
		
		/*Dealing with a sparse sequential matrix: build ir, jc and val, as though it was a row comvalessed 
		 *format, instead of a column comvalessed format. We'll transpose later.*/

		/*First get nnz: */
		nnz=0;
		for(i=0;i<rows;i++){
			MatGetRow(matrix,i,&ncols,&columns,&column_values);
			nnz+=ncols;
			MatRestoreRow(matrix,i,&ncols,&columns,&column_values);
		}

		if(nnz){
			nzmax=nnz;
		}
		else{
			nzmax=1; //so a null matrix can be returned.
		}

		val=(double*)xcalloc(nzmax,sizeof(double));
		col_ind=(mwIndex*)xcalloc(nzmax,sizeof(mwIndex));
		row_ptr=(mwIndex*)xcalloc((rows+1),sizeof(mwIndex));

		j=0;

		for(i=0;i<rows;i++){

			/*Get row from petsc matrix: */
			MatGetRow(matrix,i,&ncols,&columns,&column_values);

			#ifdef _DEBUG_
			for(j=0;j<ncols;j++){
				printf("%i %i: %g\n",i,columns[j],column_values[j]);
			}
			#endif

			/*copy values: */
			if(ncols)memcpy( val+j, column_values,ncols*sizeof(double));
			  
			for(k=0;k<ncols;k++) col_ind[j+k]=(mwIndex)columns[k];
			row_ptr[i]=(mwIndex)j;
			
			j+=ncols;
			
			/*restore petsc row, otherwise we are leaking memory: */
			MatRestoreRow(matrix,i,&ncols,&columns,&column_values);
		}
		row_ptr[rows]=(mwIndex)nnz;

		/*Ok, allocate arrays: */
		dataref = mxCreateSparse((mwSize)0,(mwSize)0,(mwSize)0,mxREAL);
		
		if(nnz){
			/* free first to avoid mem leaks...: */
			mxFree(mxGetData(dataref));
			mxFree(mxGetIr(dataref));
			mxFree(mxGetJc(dataref));

			/* ...then set data: */
			mxSetM(dataref,(mwSize)cols);
			mxSetN(dataref,(mwSize)rows);
			mxSetNzmax(dataref,(mwSize)nzmax);
			mxSetData( dataref, val);
			mxSetIr(dataref,col_ind);
			mxSetJc(dataref,row_ptr);
		}
	

	}
	else{
		/*Dealing with a dense sequential matrix: recover pointer to the data in the Petsc matrix*/

		idxm=(int*)xmalloc(rows*sizeof(int));
		idxn=(int*)xmalloc(cols*sizeof(int));

		for(i=0;i<rows;i++)idxm[i]=i;
		for(i=0;i<cols;i++)idxn[i]=i;

		matrix_ptr=(double*)xmalloc(rows*cols*sizeof(double));
		MatGetValues(matrix,rows,idxm,cols,idxn,matrix_ptr);

		dataref = mxCreateDoubleMatrix(0,0,mxREAL);
		mxSetM(dataref,rows);
		mxSetN(dataref,cols);
		mxSetPr(dataref,(double*)matrix_ptr);	
	}


	/*Transpose matrcol_indes: */
	mexCallMATLAB(1,&tdataref,1,&dataref, "'");


	/*Assign output pointers: */
	*pdataref=tdataref;

	return 1;
}

#endif  //#ifdef _SERIAL_
