/*!\file Vector.cpp
 * \brief: implementation of the Vector 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 "./Vector.h"
/*}}}*/

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

	#ifdef _HAVE_PETSC_
	this->vector=NULL;
	#else
	this->vector=NULL;
	#endif
	#ifdef _HAVE_ADOLC_
	this->avector=NULL;
	#endif
}
/*}}}*/
/*FUNCTION Vector::Vector(int M,bool fromlocalsize){{{*/
Vector::Vector(int pM,bool fromlocalsize){

	#ifdef _HAVE_PETSC_
	this->vector=NewVec(pM,fromlocalsize);
	#else
	this->vector=new SeqVec(pM,fromlocalsize);
	#endif
	#ifdef _HAVE_ADOLC_
	this->avector=xNew<IssmDouble>(pM);
	#endif
}
/*}}}*/
/*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M){{{*/
Vector::Vector(IssmDouble* serial_vec,int M){

	#ifdef _HAVE_PETSC_
		int* idxm=xNew<int>(M);
		for(int i=0;i<M;i++) idxm[i]=i;

		this->vector=NewVec(M);
		VecSetValues(this->vector,M,idxm,serial_vec,INSERT_VALUES);
		VecAssemblyBegin(this->vector);
		VecAssemblyEnd(this->vector);

		xDelete<int>(idxm);
	#else
		this->vector=new SeqVec(serial_vec,M);
	#endif
	#ifdef _HAVE_ADOLC_
		this->avector=xNew<IssmDouble>(M);
	#endif
}
/*}}}*/
#ifdef _HAVE_PETSC_
/*FUNCTION Vector::Vector(Vec petsc_vec){{{*/
Vector::Vector(Vec petsc_vec){

	if(petsc_vec==NULL){
		this->vector=NewVec(0);
	}
	else{
		/*copy vector*/
		VecDuplicate(petsc_vec,&this->vector);
		VecCopy(petsc_vec,this->vector);
	}

}
/*}}}*/
#endif
#if defined(_HAVE_GSL_) && !defined(_HAVE_PETSC_)
/*FUNCTION Vector::Vector(SeqVec* seq_vec){{{*/
Vector::Vector(SeqVec*  seq_vec){

	if(seq_vec==NULL){
		this->vector=NULL;
	}
	else{
		/*copy vector*/
		this->vector=seq_vec->Duplicate();
	}
}
/*}}}*/
#endif

		/*FUNCTION Vector::~Vector(){{{*/
Vector::~Vector(){

 	#ifdef _HAVE_PETSC_
	VecFree(&this->vector);
	#else
	delete this->vector;
	#endif
	#ifdef _HAVE_ADOLC_
	xDelete<IssmDouble>(this->avector);
	#endif
}
/*}}}*/

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

	int i;

	#ifdef _HAVE_PETSC_
	_assert_(this->vector);
	VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD);
	#else
	this->vector->Echo();
	#endif

	#ifdef _HAVE_ADOLC_
	/*do nothing for now: */
	#endif
}
/*}}}*/
/*FUNCTION Vector::Assemble{{{*/
void Vector::Assemble(void){
		
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecAssemblyBegin(this->vector); 
		VecAssemblyEnd(this->vector);
	#else
		this->vector->Assemble();
	#endif

}
/*}}}*/
/*FUNCTION Vector::SetValues{{{*/
void Vector::SetValues(int ssize, int* list, IssmDouble* values, InsMode mode){
		
		
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode));
	#else
		this->vector->SetValues(ssize,list,values,mode);
	#endif

}
/*}}}*/
/*FUNCTION Vector::SetValue{{{*/
void Vector::SetValue(int dof, IssmDouble value, InsMode mode){
		
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode));
	#else
		this->vector->SetValue(dof,value,mode);
	#endif

}
/*}}}*/
/*FUNCTION Vector::GetValue{{{*/
void Vector::GetValue(IssmDouble* pvalue,int dof){
		
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecGetValues(this->vector,1,&dof,pvalue);
	#else
	this->vector->GetValue(pvalue,dof);
	#endif
}
/*}}}*/
/*FUNCTION Vector::GetSize{{{*/
void Vector::GetSize(int* pM){
		
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecGetSize(this->vector,pM);
	#else
		this->vector->GetSize(pM);
	#endif

}
/*}}}*/
/*FUNCTION Vector::IsEmpty{{{*/
bool Vector::IsEmpty(void){

	int M;

	_assert_(this->vector);
	this->GetSize(&M);

	if(M==0) 
	 return true;
	else
	 return false;
}
/*}}}*/
/*FUNCTION Vector::GetLocalSize{{{*/
void Vector::GetLocalSize(int* pM){
		
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecGetLocalSize(this->vector,pM);
	#else
		this->vector->GetLocalSize(pM);
	#endif

}
/*}}}*/
/*FUNCTION Vector::Duplicate{{{*/
Vector* Vector::Duplicate(void){
	
	Vector* output=NULL;
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		Vec vec_output=NULL;
		VecDuplicate(this->vector,&vec_output);
		output=new Vector(vec_output);
		VecFree(&vec_output);
	#else
		output=new Vector();
		output->vector=this->vector->Duplicate();
	#endif
	return output;

}
/*}}}*/
/*FUNCTION Vector::Set{{{*/
void Vector::Set(IssmDouble value){
	
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecSet(this->vector,value);
	#else
		this->vector->Set(value);
	#endif

}
/*}}}*/
/*FUNCTION Vector::AXPY{{{*/
void Vector::AXPY(Vector* X, IssmDouble a){
	
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecAXPY(this->vector,a,X->vector);
	#else
		this->vector->AXPY(X->vector,a);
	#endif
}
/*}}}*/
/*FUNCTION Vector::AYPX{{{*/
void Vector::AYPX(Vector* X, IssmDouble a){
	
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecAYPX(this->vector,a,X->vector);
	#else
		this->vector->AYPX(X->vector,a);
	#endif

}
/*}}}*/
/*FUNCTION Vector::ToMPISerial{{{*/
IssmDouble* Vector::ToMPISerial(void){

	IssmDouble* vec_serial=NULL;

	#ifdef _HAVE_PETSC_
		VecToMPISerial(&vec_serial, this->vector);
	#else
		vec_serial=this->vector->ToMPISerial();
	#endif

	return vec_serial;

}
/*}}}*/
/*FUNCTION Vector::Copy{{{*/
void Vector::Copy(Vector* to){

	#ifdef _HAVE_PETSC_
		if(this->vector) VecCopy(this->vector,to->vector);
	#else
		this->vector->Copy(to->vector);
	#endif

}
/*}}}*/
/*FUNCTION Vector::Norm{{{*/
IssmDouble Vector::Norm(NormMode norm_type){
	
	IssmDouble norm=0;
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecNorm(this->vector,ISSMToPetscNormMode(norm_type),&norm);
	#else
		norm=this->vector->Norm(norm_type);
	#endif
	return norm;
}
/*}}}*/
/*FUNCTION Vector::Scale{{{*/
void Vector::Scale(IssmDouble scale_factor){
	
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecScale(this->vector,scale_factor); 
	#else
		this->vector->Scale(scale_factor);
	#endif
}
/*}}}*/
/*FUNCTION Vector::Dot{{{*/
IssmDouble Vector::Dot(Vector* vector){

	IssmDouble dot;
	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecDot(this->vector,vector->vector,&dot);
	#else
		dot=this->vector->Dot(vector->vector);
	#endif
	return dot;
}
/*}}}*/
/*FUNCTION Vector::PointwiseDivide{{{*/
void Vector::PointwiseDivide(Vector* x,Vector* y){

	#ifdef _HAVE_PETSC_
		_assert_(this->vector);
		VecPointwiseDivide(this->vector,x->vector,y->vector);
	#else
		this->vector->PointwiseDivide(x->vector,y->vector);
	#endif
}
/*}}}*/
