/*!\file MatToSerial.cpp
 * \brief gather a Petsc Mat matrix onto cpu 0
 */

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

#undef __FUNCT__ 
#define __FUNCT__ "MatToSerial"
void MatToSerial(double** poutmatrix,Mat matrix){

	int i;
	extern int my_rank;
	extern int num_procs;
	
	/*Petsc variables*/
	PetscInt lower_row,upper_row; 
	int range;
	int M,N; //size of matrix
	MPI_Status status;
	int* idxm=NULL;
	int* idxn=NULL; 
	double* local_matrix=NULL; /*matrix local to each node used for temporary holding matrix values*/
	int buffer[3];
	
	/*Output*/
	double* outmatrix=NULL;
	
	
	/*get matrix size: */
	MatGetSize(matrix,&M,&N);

	/*partition: */
	MatGetOwnershipRange(matrix,&lower_row,&upper_row);    
	upper_row--; 
	range=upper_row-lower_row+1;
	
	/*Local and global allocation*/
	if (my_rank==0)outmatrix=(double*)xmalloc(M*N*sizeof(double));
	
	if (range){
		local_matrix=(double*)xmalloc(N*range*sizeof(double));
		idxm=(int*)xmalloc(range*sizeof(int));  
		idxn=(int*)xmalloc(N*sizeof(int));  
	 
		for (i=0;i<N;i++){
			*(idxn+i)=i;
		}
		for (i=0;i<range;i++){
			*(idxm+i)=lower_row+i;
		}
	
		MatGetValues(matrix,range,idxm,N,idxn,local_matrix);     
	}

	/*Now each node holds its local_matrix containing range rows. 
	 * We send these rows to the matrix on node 0*/
	
	for (i=1;i<num_procs;i++){
		if (my_rank==i){ 
			buffer[0]=my_rank;
			buffer[1]=lower_row;
			buffer[2]=range;
			MPI_Send(buffer,3,MPI_INT,0,1,MPI_COMM_WORLD);   
			if (range)MPI_Send(local_matrix,N*range,MPI_DOUBLE,0,1,MPI_COMM_WORLD); 
		}
		if (my_rank==0){
			MPI_Recv(buffer,3,MPI_INT,i,1,MPI_COMM_WORLD,&status); 
			if (buffer[2])MPI_Recv(outmatrix+(buffer[1]*N),N*buffer[2],MPI_DOUBLE,i,1,MPI_COMM_WORLD,&status);
		}
	} 
	if (my_rank==0){ 
		//Still have the local_matrix on node 0 to take care of.
		memcpy(outmatrix,local_matrix,N*range*sizeof(double));

	} 
	
	/*Assign output pointer: */
	*poutmatrix=outmatrix;
	
	xfree((void**)&idxm);
	xfree((void**)&idxn);
	xfree((void**)&local_matrix);

}
