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

/*Headers:*/
/*{{{1*/
#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(){{{1*/
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){{{1*/
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=(adouble*)xmalloc(M*N*sizeof(adouble));
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N,double sparsity){{{1*/
Matrix::Matrix(int M,int N,double sparsity){

	#ifdef _HAVE_PETSC_
	this->matrix=NewMat(M,N,sparsity);
	#else
	this->matrix=new SeqMat(M,N,sparsity);
	#endif
	#ifdef _HAVE_ADOLC_
	this->amatrix=(adouble*)xmalloc(M*N*sizeof(adouble));
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(double* serial_mat, int M,int N,double sparsity){{{1*/
Matrix::Matrix(double* serial_mat, int M,int N,double sparsity){

	#ifdef _HAVE_PETSC_
	int     i;


	int* idxm=(int*)xmalloc(M*sizeof(int));
	int* idxn=(int*)xmalloc(N*sizeof(int));
	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);

	xfree((void**)&idxm);
	xfree((void**)&idxn);
	#else
	this->matrix=new SeqMat(serial_mat,M,N,sparsity);
	#endif
	#ifdef _HAVE_ADOLC_
	this->amatrix=(adouble*)xmalloc(M*N*sizeof(adouble));
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){{{1*/
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=(adouble*)xmalloc(M*N*sizeof(adouble));
	#endif
}
/*}}}*/
/*FUNCTION Matrix::~Matrix(){{{1*/
Matrix::~Matrix(){

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

/*Matrix specific routines: */
/*FUNCTION Matrix::Echo{{{1*/
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 >>?*/
	printf("ADOLC Matrix equivalent:" );
	for(i=0;i<M;i++){
		for(j=0;j<N;j++){
			printf("%g ",*(amatrix+N*i+j));
		}
		printf("\n");
	}
	#endif
}
/*}}}*/
/*FUNCTION Matrix::Assemble{{{1*/
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{{{1*/
double Matrix::Norm(NormMode norm_type){
	
	double 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{{{1*/
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{{{1*/
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{{{1*/
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{{{1*/
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{{{1*/
double* Matrix::ToSerial(void){

	double* output=NULL;

	#ifdef _HAVE_PETSC_
		MatToSerial(&output,this->matrix);
	#else
		output=this->matrix->ToSerial();
	#endif
	return output;
}
/*}}}*/
/*FUNCTION Matrix::SetValues{{{1*/
void Matrix::SetValues(int m,int* idxm,int n,int* idxn,double* 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{{{1*/
void Matrix::Convert(MatrixType type){

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