/*!\file:  SeqVec.h
 * \brief wrapper to our SeqVec object, which is just a wrapper to a IssmDouble* 
 */ 

#ifndef _SEQVEC_H_
#define _SEQVEC_H_

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

#include "../../shared/Exceptions/exceptions.h"
#include "../../shared/MemOps/xMemCpy.h"
#include "../../shared/Alloc/alloc.h"
#include "../../include/macros.h"

/*}}}*/

/*We need to template this class, in case we want to create vectors that hold IssmDouble* matrix or IssmPDouble* matrix. 
  Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, 
  but only rely on IssmPDouble types)*/

template <class doubletype> 
class SeqVec{

	public:

		doubletype* vector;
		int M;

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

			this->M=0;
			this->vector=NULL;
		}
		/*}}}*/
		/*FUNCTION SeqVec(int M){{{*/
		SeqVec(int pM){

			this->M=pM;
			this->vector=NULL;
			if(this->M) this->vector=xNewZeroInit<doubletype>(pM);
		}
		/*}}}*/
		/*FUNCTION SeqVec(int m,int M){{{*/
		SeqVec(int pm,int pM){

			this->M=pM;
			this->vector=NULL;
			if(this->M) this->vector=xNewZeroInit<doubletype>(pM);
		}
		/*}}}*/
		/*FUNCTION SeqVec(int M,bool fromlocalsize){{{*/
		SeqVec(int pM,bool fromlocalsize){

			this->M=pM;
			this->vector=NULL;
			if(this->M) this->vector=xNewZeroInit<doubletype>(pM);
		}
		/*}}}*/
		/*FUNCTION SeqVec(doubletype* serial_vec,int M){{{*/
		SeqVec(doubletype* buffer,int pM){

			this->M=pM;
			this->vector=NULL;
			if(this->M){
				this->vector=xNew<doubletype>(pM);
				xMemCpy<doubletype>(this->vector,buffer,pM);
			}
		}
		/*}}}*/
		/*FUNCTION ~SeqVec(){{{*/
		~SeqVec(){
			xDelete<doubletype>(this->vector);
			M=0;
		}
		/*}}}*/

		/*SeqVec specific routines*/
		/*FUNCTION Echo{{{*/
		void Echo(void){

			int i;
			_printLine_("SeqVec size " << this->M);
			for(i=0;i<M;i++){
				_printString_(vector[i] << "\n ");
			}
		}
		/*}}}*/
		/*FUNCTION Assemble{{{*/
		void Assemble(void){

			/*do nothing*/

		}
		/*}}}*/
		/*FUNCTION SetValues{{{*/
		void SetValues(int ssize, int* list, doubletype* values, InsMode mode){

			int i;
			switch(mode){
				case ADD_VAL:
					for(i=0;i<ssize;i++) this->vector[list[i]]+=values[i];
					break;
				case INS_VAL:
					for(i=0;i<ssize;i++) this->vector[list[i]]=values[i];
					break;
				default:
					_error_("unknown insert mode!");
					break;
			}

		}
		/*}}}*/
		/*FUNCTION SetValue{{{*/
		void SetValue(int dof, doubletype value, InsMode mode){

			switch(mode){
				case ADD_VAL:
					this->vector[dof]+=value;
					break;
				case INS_VAL:
					this->vector[dof]=value;
					break;
				default:
					_error_("unknown insert mode!");
					break;
			}
		}
		/*}}}*/
		/*FUNCTION GetValue{{{*/
		void GetValue(doubletype* pvalue,int dof){

			*pvalue=this->vector[dof];

		}
		/*}}}*/
		/*FUNCTION GetSize{{{*/
		void GetSize(int* pM){

			*pM=this->M;

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

			*pM=this->M;

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

			return new SeqVec(this->vector,this->M);

		}
		/*}}}*/
		/*FUNCTION Set{{{*/
		void Set(doubletype value){

			int i;
			for(i=0;i<this->M;i++)this->vector[i]=value;

		}
		/*}}}*/
		/*FUNCTION AXPY{{{*/
		void AXPY(SeqVec* X, doubletype a){

			int i;

			/*y=a*x+y where this->vector is y*/
			for(i=0;i<this->M;i++)this->vector[i]=a*X->vector[i]+this->vector[i];

		}
		/*}}}*/
		/*FUNCTION AYPX{{{*/
		void AYPX(SeqVec* X, doubletype a){

			int i;

			/*y=x+a*y where this->vector is y*/
			for(i=0;i<this->M;i++)this->vector[i]=X->vector[i]+a*this->vector[i];

		}
		/*}}}*/
		/*FUNCTION ToMPISerial{{{*/
		doubletype* ToMPISerial(void){

			doubletype* buffer=NULL;

			if(this->M){
				buffer=xNew<doubletype>(this->M);
				xMemCpy<doubletype>(buffer,this->vector,this->M);
			}
			return buffer;

		}
		/*}}}*/
		/*FUNCTION Copy{{{*/
		void Copy(SeqVec* to){

			int i;

			to->M=this->M;
			for(i=0;i<this->M;i++)to->vector[i]=this->vector[i];

		}
		/*}}}*/
		/*FUNCTION Norm{{{*/
		doubletype Norm(NormMode mode){

			doubletype norm;
			int i;

			switch(mode){
				case NORM_INF:
					norm=0; for(i=0;i<this->M;i++)norm=max(norm,fabs(this->vector[i]));
					return norm;
					break;
				case NORM_TWO:
					norm=0; 
					for(i=0;i<this->M;i++)norm+=pow(this->vector[i],2);
					return sqrt(norm);
					break;
				default:
					_error_("unknown norm !");
					break;
			}
		}
		/*}}}*/
		/*FUNCTION Scale{{{*/
		void Scale(doubletype scale_factor){

			int i;
			for(i=0;i<this->M;i++)this->vector[i]=scale_factor*this->vector[i];

		}
		/*}}}*/
		/*FUNCTION Dot{{{*/
		doubletype Dot(SeqVec* input){

			int i;

			doubletype dot=0;
			for(i=0;i<this->M;i++)dot+=this->vector[i]*input->vector[i];
			return dot;

		}
		/*}}}*/
		/*FUNCTION PointwiseDivide{{{*/
		void PointwiseDivide(SeqVec* x,SeqVec* y){

			int i;
			/*pointwise w=x/y where this->vector is w: */
			for(i=0;i<this->M;i++)this->vector[i]=x->vector[i]/y->vector[i];
		}
		/*}}}*/
};
#endif //#ifndef _SEQVEC_H_
