/*!\file PlapackToPetsc.cpp
 * \brief convert plapack matrix to petsc matrix
 */

#include "../../petsc/petscincludes.h"
#include "../plapackincludes.h"

int PlapackToPetsc(Mat* A,int local_mA,int local_nA,int mA,int nA,MatType type,PLA_Obj a,PLA_Template templ,int nprows,int npcols,int nb,COMM comm){

	int i;

	int lower_row,upper_row;
	int* row_nodes=NULL;
	int* col_nodes=NULL;
	int* idxm=NULL;
	int* idxn=NULL;
	int myrow,mycol;
	int count,idxm_count,idxn_count;
	int d_nz,o_nz;
	int i0,i1;
	double* local_buffer=NULL;

	/*Create matrix A (right now, we just have an allocated pointer A*/
	if (strcasecmp_eq(type,MATMPIAIJ)){
		/*See pretty large here, because we have no idea how many nnz per row*/
		d_nz=nA/2;
		o_nz=nA/2;
		MatCreateMPIAIJ(comm,local_mA,local_nA, mA,nA,d_nz,PETSC_NULL,o_nz,PETSC_NULL,A);
	}
	else if(strcasecmp_eq(type,MATMPIDENSE)){
		MatCreateMPIDense(comm,local_mA,local_nA, mA,nA,PETSC_NULL,A); 
	}

	MatGetOwnershipRange(*A,&lower_row,&upper_row);
	upper_row--;

	/*Build the Plapack row and column indices corresponding to the local_buffer stored in a. 
	  We need those indices to directly plug local_buffer into the Petsc matrix A. We do not 
	  use PLA_axpy_global_to_matrix to extract a local matrix, because this routine hangs when the 
	  problem size becomes big. We rely therefore on MatAssembly from Petsc to gather the plapack 
	  matrix into a Petsc Matrix.*/

	/*Vector physically based block cyclic distribution: */
	row_nodes=xNew<int>(mA);
	col_nodes=xNew<int>(nA);
	for (i=0;i<mA;i++){
		i0=i/nb;
		*(row_nodes+i)=i0%nprows;
	}
	for (i=0;i<nA;i++){
		i0=i/nb;
		i1=i0/nprows;
		*(col_nodes+i)=i1%npcols;
	}

	PLA_Temp_comm_row_rank(templ,&mycol);
	PLA_Temp_comm_col_rank(templ,&myrow);

	idxm=xNew<int>(mA);
	count=0;
	for (i=0;i<mA;i++){
		if(*(row_nodes+i)==myrow){
			*(idxm+count)=i;
			count++;
		}
	}
	idxm_count=count;

	idxn=xNew<int>(nA);
	count=0;
	for (i=0;i<nA;i++){
		if(*(col_nodes+i)==mycol){
			*(idxn+count)=i;
			count++;
		}
	}
	idxn_count=count;

	/*Get local buffer: */
	PLA_Obj_local_buffer(a,(void**)&local_buffer);

	/*Insert into invA matrix. Use col oriented insertion, for Plapack is column oriented*/
	MatSetOption(*A,MAT_COLUMN_ORIENTED);
	MatSetValues(*A,idxm_count,idxm,idxn_count,idxn,local_buffer,INSERT_VALUES);

	/*assemble: */
	MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);

	/*Free ressources:*/
	xDelete<int>(row_nodes);
	xDelete<int>(col_nodes);
	xDelete<int>(idxm);
	xDelete<int>(idxn);
}
