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


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

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

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

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

int MatlabVectorToPetscVector(Vec* pvector,int* pvector_rows,const mxArray* mxvector){

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

	/*output: */
	Vec vector=NULL;

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

	/*petsc indices: */
	int* idxm=NULL;
	
	/*Ok, first check if we are dealing with a sparse or full vector: */
	if (mxIsSparse(mxvector)){

		/*Dealing with sparse vector: recover size first: */
		mxvector_ptr=(double*)mxGetPr(mxvector);
		rows=mxGetM(mxvector);
		cols=mxGetN(mxvector);
		nnz=mxGetNzmax(mxvector);
		nz=(int)((double)nnz/(double)rows);

		ierr=VecCreateSeq(PETSC_COMM_SELF,rows,&vector);CHKERRQ(ierr);

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

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

	}
	else{

		/*Dealing with dense vector: recover pointer and size: */
		mxvector_ptr=(double*)mxGetPr(mxvector);
		rows=mxGetM(mxvector);
		cols=mxGetN(mxvector);

		/*Create serial vector: */
		ierr=VecCreateSeq(PETSC_COMM_SELF,rows,&vector);CHKERRQ(ierr);

		/*Insert mxvector_ptr values into petsc vector: */
		idxm=xNew<int>(rows);

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

		ierr=VecSetValues(vector,rows,idxm,mxvector_ptr,INSERT_VALUES);CHKERRQ(ierr);

	}

	/*Assemble vector: */
	VecAssemblyBegin(vector);
	VecAssemblyEnd(vector);

	/*Assign output pointer: */
	*pvector=vector;
	*pvector_rows=rows;

	return 1;
}
