/*!\file Matrix.cpp
 * \brief: implementation of the Matrix object
 */

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

#include <stdio.h>
#include <string.h>
#include "../classes.h"
#include "../../shared/shared.h"
#include "../../Container/Container.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "./Matrix.h"
	
/*}}}*/

/*Matrix constructors and destructor*/
/*FUNCTION Matrix::Matrix(){{{*/
Matrix::Matrix(){

	pmatrix=NULL;
	smatrix=NULL;

	type=PetscMatType; //default
	#ifndef _HAVE_PETSC_
	type=SeqMatType;
	#endif
	
}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N,int type){{{*/
Matrix::Matrix(int M,int N,int in_type){

	pmatrix=NULL;
	smatrix=NULL;
	type=in_type;

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

}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N,IssmDouble sparsity,int type){{{*/
Matrix::Matrix(int M,int N,IssmDouble sparsity,int in_type){

	pmatrix=NULL;
	smatrix=NULL;
	type=in_type;

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

	pmatrix=NULL;
	smatrix=NULL;
	type=in_type;

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

	pmatrix=NULL;
	smatrix=NULL;
	type=in_type;

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

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

}
/*}}}*/

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

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

}
/*}}}*/
/*FUNCTION Matrix::Assemble{{{*/
void Matrix::Assemble(void){

	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		this->pmatrix->Assemble();
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		this->smatrix->Assemble();
	}
	else{
		_error2_("Matrix type: " << type << " not supported yet!");
	}
}
/*}}}*/
/*FUNCTION Matrix::Norm{{{*/
IssmDouble Matrix::Norm(NormMode norm_type){
	
	IssmDouble norm=0;

	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		norm=this->pmatrix->Norm(norm_type);
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		norm=this->smatrix->Norm(norm_type);
	}
	else _error2_("Matrix type: " << type << " not supported yet!");

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

	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		this->pmatrix->GetSize(pM,pN);
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		this->smatrix->GetSize(pM,pN);
	}
	else _error2_("Matrix type: " << type << " not supported yet!");
	
}
/*}}}*/
/*FUNCTION Matrix::GetLocalSize{{{*/
void Matrix::GetLocalSize(int* pM,int* pN){
	
	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		this->pmatrix->GetLocalSize(pM,pN);
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		this->smatrix->GetLocalSize(pM,pN);
	}
	else _error2_("Matrix type: " << type << " not supported yet!");

}
/*}}}*/
/*FUNCTION Matrix::MatMult{{{*/
void Matrix::MatMult(Vector* X,Vector* AX){

	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		this->pmatrix->MatMult(X->pvector,AX->pvector);
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		this->smatrix->MatMult(X->svector,AX->svector);
	}
	else _error2_("Matrix type: " << type << " not supported yet!");

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

	Matrix* output=NULL;

	output=new Matrix();

	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		output->pmatrix=this->pmatrix->Duplicate();
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		output->smatrix=this->smatrix->Duplicate();
	}
	else _error2_("Matrix type: " << type << " not supported yet!");
	
	return output;
}
/*}}}*/
/*FUNCTION Matrix::ToSerial{{{*/
IssmDouble* Matrix::ToSerial(void){

	IssmDouble* output=NULL;
	
	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		output=this->pmatrix->ToSerial();
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		output=this->smatrix->ToSerial();
	}
	else _error2_("Matrix type: " << type << " not supported yet!");


	return output;
}
/*}}}*/
/*FUNCTION Matrix::SetValues{{{*/
void Matrix::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);
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		this->smatrix->SetValues(m,idxm,n,idxn,values,mode);
	}
	else _error2_("Matrix type: " << type << " not supported yet!");
}
/*}}}*/
/*FUNCTION Matrix::Convert{{{*/
void Matrix::Convert(MatrixType newtype){

	if(type==PetscMatType){
		#ifdef _HAVE_PETSC_
		this->pmatrix->Convert(newtype);
		#else
		_error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
		#endif
	}
	else if(type==SeqMatType){
		this->smatrix->Convert(newtype);
	}
	else{
		_error2_("Matrix type: " << type << " not supported yet!");
	}

}
/*}}}*/
