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

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


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

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=xNew<double>(M*N);
	
	if (range){
		local_matrix=xNew<double>(N*range);
		idxm=xNew<int>(range);  
		idxn=xNew<int>(N);  
	 
		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;
	xDelete<int>(idxm);
	xDelete<int>(idxn);
	xDelete<double>(local_matrix);
}
