/*!\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"


Mat NewMat(int M,int N,double* psparsity,int* pconnectivity,int* pnumberofdofspernode){

	/*output:*/
	Mat outmatrix=NULL;

	/*Recover parameters: */
	double sparsity=.001; //default

	int    connectivity=10; //default
	int    numberofdofspernode=2; //default

	/*local sizes: */
	int    m,n;

	/*row and col distribution: */
	int    d_nz,o_nz;
	int    nnz;

	/*Determine local sizes: */
	m=PetscDetermineLocalSize(M);
	n=PetscDetermineLocalSize(N);
	
	/*Branch on what was provided: */

	if(!psparsity && !pconnectivity && !pnumberofdofspernode){

		/*Ok, we have our default sparsity, create matrix using this info: */

		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

		MatCreateMPIAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 

	}
	else if(psparsity){
		
		/*Same as above, except we bypass the default sparsity: */
		if(pconnectivity || pnumberofdofspernode)
			throw ErrorException(__FUNCT__,exprintf("%s%p%s%p%s%p","cannot determine logic with following set of poniters-> sparsity: ",
						psparsity," connectivity: ",pconnectivity," numberofdofspernode: ",pnumberofdofspernode));

		sparsity=*psparsity;

		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

		MatCreateMPIAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 

	}
	else if(pconnectivity &&  pnumberofdofspernode){

		/*entirely different logic. We use the connectivity and numberofdofspernode to figure out how many non-zeros are 
		 *present per row: */

		if(psparsity)throw ErrorException(__FUNCT__,exprintf("%s%p%s%p%s%p","cannot determine logic with following set of poniters-> sparsity: ",
						psparsity," connectivity: ",pconnectivity," numberofdofspernode: ",pnumberofdofspernode));

		connectivity=*pconnectivity;
		numberofdofspernode=*pnumberofdofspernode;
	
		/*Figure out number of non zeros per row: */
		d_nz=(int)connectivity*numberofdofspernode/2;
		o_nz=(int)connectivity*numberofdofspernode/2;
		
		MatCreateMPIAIJ(MPI_COMM_WORLD,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix); 

	}
	else{
		throw ErrorException(__FUNCT__,exprintf("%s%p%s%p%s%p","cannot determine logic with following set of poniters-> sparsity: ",
						psparsity," connectivity: ",pconnectivity," numberofdofspernode: ",pnumberofdofspernode));
	}


	return outmatrix;
}
