/* \file MatlabMatrixToPetscMatrix.cpp
 * \brief: convert a sparse or dense matlab matrix to a serial Petsc matrix:
 */

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
#include "../../c/shared/shared.h"

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

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

PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix){

	int dummy;
	PetscMat* matrix=new PetscMat();

	MatlabMatrixToPetscMat(&matrix->matrix, &dummy, &dummy, mxmatrix);

	return matrix;
}
int MatlabMatrixToPetscMat(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){

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

	/*output: */
	Mat matrix = NULL;

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

	/*petsc indices: */
	int *idxm = NULL;
	int *idxn = NULL;

	/*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);
		if(rows){
			nz=(int)((double)nnz/(double)rows);
		}
		else{
			nz=0;
		}

		ierr=MatCreateSeqAIJ(PETSC_COMM_SELF,rows,cols,nz,PETSC_NULL,&matrix);CHKERRQ(ierr);

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

		/*Now, start inserting data into sparse matrix: */
		count=0;
		for(i=0;i<cols;i++){
			for(j=0;j<(jc[i+1]-jc[i]);j++){
				MatSetValue(matrix,ir[count],i,pr[count],INSERT_VALUES);
				count++;
			}
		}
	}
	else{
		/*Dealing with dense matrix: recover pointer and size: */
		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
		rows=mxGetM(mxmatrix);
		cols=mxGetN(mxmatrix);

		/*transpose, as Petsc now does not allows MAT_COLUMN_ORIENTED matrices in MatSetValues: */
		tmatrix=xNew<double>(rows*cols);
		for(i=0;i<cols;i++){
			for(j=0;j<rows;j++){
				*(tmatrix+rows*i+j)=*(mxmatrix_ptr+cols*j+i);
			}
		}

		/*Create serial matrix: */
		ierr=MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&matrix);CHKERRQ(ierr);

		/*Insert mxmatrix_ptr values into petsc matrix: */
		idxm=xNew<int>(rows);
		idxn=xNew<int>(cols);

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

		ierr=MatSetValues(matrix,rows,idxm,cols,idxn,tmatrix,INSERT_VALUES); CHKERRQ(ierr);

		xDelete<double>(tmatrix);
	}

	/*Assemble matrix: */
	MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY); 
	MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY);

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

	return 1;
}