/*!\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 "../objects.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(){

	#ifdef _HAVE_PETSC_
	this->matrix=NULL;
	#else
	this->matrix=NULL;
	#endif
	#ifdef _HAVE_ADOLC_
	this->amatrix=NULL;
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N){{{*/
Matrix::Matrix(int M,int N){

	#ifdef _HAVE_PETSC_
	this->matrix=NewMat(M,N);
	#else
	this->matrix=new SeqMat(M,N);
	#endif
	#ifdef _HAVE_ADOLC_
	this->amatrix=xNew<IssmDouble>(M*N);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N,IssmDouble sparsity){{{*/
Matrix::Matrix(int M,int N,IssmDouble sparsity){

	#ifdef _HAVE_PETSC_
	this->matrix=NewMat(M,N,sparsity);
	#else
	this->matrix=new SeqMat(M,N,sparsity);
	#endif
	#ifdef _HAVE_ADOLC_
 	this->amatrix=xNew<IssmDouble>(M*N);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/
Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){

	#ifdef _HAVE_PETSC_
	int     i;


	int* idxm=xNew<int>(M);
	int* idxn=xNew<int>(N);
	for(i=0;i<M;i++)idxm[i]=i;
	for(i=0;i<N;i++)idxn[i]=i;

	this->matrix=NewMat(M,N,sparsity);
	MatSetValues(this->matrix,M,idxm,N,idxn,serial_mat,INSERT_VALUES);
	MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY);

	xDelete<int>(idxm);
	xDelete<int>(idxn);
	#else
	this->matrix=new SeqMat(serial_mat,M,N,sparsity);
	#endif
	#ifdef _HAVE_ADOLC_
	this->amatrix=xNew<IssmDouble>(M*N);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/
Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){
	
	#ifdef _HAVE_PETSC_
	this->matrix=NewMat(M,N,connectivity,numberofdofspernode);
	#else
	this->matrix=new SeqMat(M,N,connectivity,numberofdofspernode);
	#endif
	#ifdef _HAVE_ADOLC_
	this->amatrix=xNew<IssmDouble>(M*N);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::~Matrix(){{{*/
Matrix::~Matrix(){

 	#ifdef _HAVE_PETSC_
	MatFree(&this->matrix);
	#else
	delete this->matrix;
	#endif
	#ifdef _HAVE_ADOLC_
	xDelete<IssmDouble>(this->amatrix);
	#endif
}
/*}}}*/

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

	int i,j;

	#ifdef _HAVE_PETSC_
	MatView(this->matrix,PETSC_VIEWER_STDOUT_WORLD);
	#else
	this->matrix->Echo();
	#endif

	#ifdef _HAVE_ADOLC_
	/*Not sure about that one. Should we use the overloaded operator >>?*/
	_printString_("ADOLC Matrix equivalent:" );
//	for(i=0;i<M;i++){
//		for(j=0;j<N;j++){
//			_printString_(*(amatrix+N*i+j) << " ");
//		}
//		_printLine_("");
//	}
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Assemble{{{*/
void Matrix::Assemble(void){
	#ifdef _HAVE_PETSC_
		_assert_(this->matrix);
		MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY);
		MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY);
		#if _PETSC_MAJOR_ == 2 
			MatCompress(this->matrix);
		#endif
	#else
		this->matrix->Assemble();
	#endif

}
/*}}}*/
/*FUNCTION Matrix::Norm{{{*/
IssmDouble Matrix::Norm(NormMode norm_type){
	
	IssmDouble norm=0;
	#ifdef _HAVE_PETSC_
		_assert_(this->matrix);
		MatNorm(this->matrix,ISSMToPetscNormMode(norm_type),&norm);
	#else
		norm=this->matrix->Norm(norm_type);
	#endif
	return norm;
}
/*}}}*/
/*FUNCTION Matrix::GetSize{{{*/
void Matrix::GetSize(int* pM,int* pN){
	
	#ifdef _HAVE_PETSC_
		_assert_(this->matrix);
		MatGetSize(this->matrix,pM,pN);
	#else
		this->matrix->GetSize(pM,pN);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::GetLocalSize{{{*/
void Matrix::GetLocalSize(int* pM,int* pN){
	
	#ifdef _HAVE_PETSC_
		_assert_(this->matrix);
		MatGetLocalSize(this->matrix,pM,pN);
	#else
		this->matrix->GetLocalSize(pM,pN);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::MatMult{{{*/
void Matrix::MatMult(Vector* X,Vector* AX){

	#ifdef _HAVE_PETSC_
		_assert_(this->matrix);
		_assert_(X->vector);
		MatMultPatch(this->matrix,X->vector,AX->vector);
	#else
		this->matrix->MatMult(X->vector,AX->vector);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Duplicate{{{*/
Matrix* Matrix::Duplicate(void){

	Matrix* output=NULL;

	output=new Matrix();

	#ifdef _HAVE_PETSC_
		_assert_(this->matrix);
		MatDuplicate(this->matrix,MAT_COPY_VALUES,&output->matrix);
	#else
		output->matrix=this->matrix->Duplicate();
	#endif
	
	return output;
}
/*}}}*/
/*FUNCTION Matrix::ToSerial{{{*/
IssmDouble* Matrix::ToSerial(void){

	IssmDouble* output=NULL;

	#ifdef _HAVE_PETSC_
		MatToSerial(&output,this->matrix);
	#else
		output=this->matrix->ToSerial();
	#endif
	return output;
}
/*}}}*/
/*FUNCTION Matrix::SetValues{{{*/
void Matrix::SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){

	#ifdef _HAVE_PETSC_
		MatSetValues(this->matrix,m,idxm,n,idxn,values,ISSMToPetscInsertMode(mode));
	#else
		this->matrix->SetValues(m,idxm,n,idxn,values,mode);
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Convert{{{*/
void Matrix::Convert(MatrixType type){

	#ifdef _HAVE_PETSC_
		MatConvert(this->matrix,ISSMToPetscMatrixType(type),MAT_REUSE_MATRIX,&this->matrix);
	#else
		this->matrix->Convert(type);
	#endif
}
/*}}}*/
