/*!\file:  Matrix.h
 * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
 * implements our underlying matrix format.
 */ 

#ifndef _MATRIX_H_
#define _MATRIX_H_

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
/*}}}*/

enum matrixtype { PetscMatType, SeqMatType };

template <class doubletype> class Vector;

template <class doubletype> 
class Matrix{

	public:

		#ifdef _HAVE_PETSC_
		PetscMat *pmatrix;
		#endif
		SeqMat<doubletype>   *smatrix;
		int       type;

		/*Matrix constructors, destructors*/
		/*FUNCTION Matrix(){{{*/
		Matrix(){

			#ifdef _HAVE_PETSC_
			pmatrix=NULL;
			#endif
			smatrix=NULL;

			type=PetscMatType; //default
			#ifndef _HAVE_PETSC_
			type=SeqMatType;
			#endif

		}
		/*}}}*/
		/*FUNCTION Matrix(int M,int N,int in_type){{{*/
		#ifdef _HAVE_PETSC_
		Matrix(int M,int N,int in_type=PetscMatType){
		#else
		Matrix(int M,int N,int in_type=SeqMatType){
		#endif

			#ifdef _HAVE_PETSC_
			pmatrix=NULL;
			#endif
			smatrix=NULL;
			type=in_type;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix=new PetscMat(M,N);
				#else
				_error_("Petsc matrix format not usable, as Petsc has not been compiled!");
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix=new SeqMat<doubletype>(M,N);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz,int type){{{*/
		#ifdef _HAVE_PETSC_
		Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz,int in_type=PetscMatType){
		#else
		Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz,int in_type=SeqMatType){
		#endif

			#ifdef _HAVE_PETSC_
			pmatrix=NULL;
			#endif
			smatrix=NULL;
			type=in_type;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);
				#else
				_error_("Petsc matrix format not usable, as Petsc has not been compiled!");
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix=new SeqMat<doubletype>(m,n,M,N,d_nnz,o_nnz);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION Matrix(int M,int N,IssmDouble sparsity,int in_type){{{*/
		#ifdef _HAVE_PETSC_
		Matrix(int M,int N,double sparsity,int in_type=PetscMatType){
		#else
		Matrix(int M,int N,double sparsity,int in_type=SeqMatType){
		#endif

			#ifdef _HAVE_PETSC_
			pmatrix=NULL;
			#endif

			smatrix=NULL;
			type=in_type;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix=new PetscMat(M,N,sparsity);
				#else
				_error_("Petsc matrix format not usable, as Petsc has not been compiled!");
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix=new SeqMat<doubletype>(M,N,sparsity);
			}
			else _error_("Matrix type: " << type << " not supported yet!");
		}
		/*}}}*/
		/*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int in_type){{{*/
		#ifdef _HAVE_PETSC_
		Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity,int in_type=PetscMatType){
		#else
		Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity,int in_type=SeqMatType){
		#endif

			#ifdef _HAVE_PETSC_
			pmatrix=NULL;
			#endif
			smatrix=NULL;
			type=in_type;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
				#else
				_error_("Petsc matrix format not usable, as Petsc has not been compiled!");
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix=new SeqMat<doubletype>(serial_mat,M,N,sparsity);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode,int in_type){{{*/
		#ifdef _HAVE_PETSC_
		Matrix(int M,int N,int connectivity,int numberofdofspernode,int in_type=PetscMatType){
		#else
		Matrix(int M,int N,int connectivity,int numberofdofspernode,int in_type=SeqMatType){
		#endif

			#ifdef _HAVE_PETSC_
			pmatrix=NULL;
			#endif
			smatrix=NULL;
			type=in_type;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
				#else
				_error_("Petsc matrix format not usable, as Petsc has not been compiled!");
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix=new SeqMat<doubletype>(M,N,connectivity,numberofdofspernode);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION ~Matrix(){{{*/
		~Matrix(){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				delete this->pmatrix;
				#endif
			}
			else if(type==SeqMatType){
				delete this->smatrix;
			}

		}
		/*}}}*/

		/*Matrix specific routines:*/
		/*FUNCTION Echo{{{*/
		void Echo(void){
			_assert_(this);

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->Echo();
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->Echo();
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION AllocationInfo{{{*/
		void AllocationInfo(void){
			_assert_(this);
			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->AllocationInfo();
				#endif
			}
			else if(type==SeqMatType){
				//this->smatrix->AllocationInfo();
			}
			else _error_("Matrix type: " << type << " not supported yet!");
		}/*}}}*/
		/*FUNCTION Assemble{{{*/
		void Assemble(void){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->Assemble();
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->Assemble();
			}
			else{
				_error_("Matrix type: " << type << " not supported yet!");
			}
		}
		/*}}}*/
		/*FUNCTION Norm{{{*/
		IssmDouble Norm(NormMode norm_type){

			IssmDouble norm=0;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				norm=this->pmatrix->Norm(norm_type);
				#endif
			}
			else if(type==SeqMatType){
				norm=this->smatrix->Norm(norm_type);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

			return norm;
		}
		/*}}}*/
		/*FUNCTION GetSize{{{*/
		void GetSize(int* pM,int* pN){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->GetSize(pM,pN);
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->GetSize(pM,pN);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION GetLocalSize{{{*/
		void GetLocalSize(int* pM,int* pN){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->GetLocalSize(pM,pN);
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->GetLocalSize(pM,pN);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION MatMult{{{*/
		void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->MatMult(X->pvector,AX->pvector);
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->MatMult(X->svector,AX->svector);
			}
			else _error_("Matrix type: " << type << " not supported yet!");

		}
		/*}}}*/
		/*FUNCTION Duplicate{{{*/
		Matrix* Duplicate(void){

			Matrix* output=new Matrix();

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				output->pmatrix=this->pmatrix->Duplicate();
				#endif
			}
			else if(type==SeqMatType){
				output->smatrix=this->smatrix->Duplicate();
			}
			else _error_("Matrix type: " << type << " not supported yet!");

			return output;
		}
		/*}}}*/
		/*FUNCTION ToSerial{{{*/
		doubletype* ToSerial(void){

			doubletype* output=NULL;

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				output=this->pmatrix->ToSerial();
				#endif
			}
			else if(type==SeqMatType){
				output=this->smatrix->ToSerial();
			}
			else _error_("Matrix type: " << type << " not supported yet!");

			return output;
		}
		/*}}}*/
		/*FUNCTION SetValues{{{*/
		void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->SetValues(m,idxm,n,idxn,values,mode);
			}
			else _error_("Matrix type: " << type << " not supported yet!");
		}
		/*}}}*/
		/*FUNCTION Convert{{{*/
		void Convert(MatrixType newtype){

			if(type==PetscMatType){
				#ifdef _HAVE_PETSC_
				this->pmatrix->Convert(newtype);
				#endif
			}
			else if(type==SeqMatType){
				this->smatrix->Convert(newtype);
			}
			else{
				_error_("Matrix type: " << type << " not supported yet!");
			}

		}
		/*}}}*/
};

#endif //#ifndef _MATRIX_H_
