/*!\file:  NewMat.cpp
 * \brief create matrix using the Petsc library
 */ 

#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"

#include "./petscpatches.h"

#include "../../../shared/shared.h"
#include "../../../include/include.h"
#include "../../mpi/patches/mpipatches.h"

/*NewMat(int M,int N){{{*/
Mat NewMat(int M,int N){

	/*output:*/
	Mat outmatrix=NULL;

	/*parameters: */
	double sparsity=.001; //default
	int    m,n;
	int    d_nz,o_nz,nnz;

	/*Determine local sizes: */
	m=DetermineLocalSize(M);
	n=DetermineLocalSize(N);
	
	nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
	d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
	o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2

	#ifdef _HAVE_PETSCDEV_
	MatCreateAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 
	#else
	MatCreateMPIAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 
	#endif

	return outmatrix;
}
/*}}}*/
/*NewMat(int M,int N,double sparsity){{{*/
Mat NewMat(int M,int N,double sparsity){

	/*output:*/
	Mat outmatrix=NULL;

	/*parameters: */
	int    m,n;
	int    d_nz,o_nz;
	int    nnz;

	/*Determine local sizes: */
	m=DetermineLocalSize(M);
	n=DetermineLocalSize(N);
	
	nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
	d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
	o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2

	#ifdef _HAVE_PETSCDEV_
	if(sparsity==1){
		MatCreateDense(MPI_COMM_WORLD,m,n,M,N,NULL,&outmatrix); 
	}
	else{
		MatCreateAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 
	}
	#else
	MatCreateMPIAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 
	#endif

	return outmatrix;
}
/*}}}*/
/*NewMat(int M,int N,int connectivity,int numberofdofspernode){{{*/
Mat NewMat(int M,int N,int connectivity,int numberofdofspernode){

	/*output:*/
	Mat outmatrix=NULL;

	/*parameters: */
	int    m,n;
	int    d_nz,o_nz;
	int    nnz;

	#if _PETSC_MAJOR_ >= 3 
	const MatType type;
	#else
	MatType type;
	#endif

	/*Determine local sizes: */
	m=DetermineLocalSize(M);
	n=DetermineLocalSize(N);

	/*Figure out number of non zeros per row: */
	d_nz=(int)connectivity*numberofdofspernode/2;
	o_nz=(int)connectivity*numberofdofspernode/2;

	MatCreate(MPI_COMM_WORLD,&outmatrix);
	MatSetSizes(outmatrix,m,n,M,N);
	MatSetFromOptions(outmatrix);

	/*preallocation  according to type: */
	MatGetType(outmatrix,&type);
	
	#if _PETSC_MAJOR_ == 2 
	if((strcmp(type,"mpiaij")==0) || (strcmp(type,"aijmumps")==0)){
		MatMPIAIJSetPreallocation(outmatrix,d_nz,NULL,o_nz,NULL);
	}
	#else
	if((strcmp(type,"mpiaij")==0) || (strcmp(type,"mpidense")==0)){
		MatMPIAIJSetPreallocation(outmatrix,d_nz,NULL,o_nz,NULL);
	}
	#endif

	return outmatrix;
}
/*}}}*/
