/*!\file:  IssmSeqVec.h
 * \brief implementation of an ISSM vector which run serially (1 cpu only), which is made of a fully dense 
 * vector. Internally, this dense vector is just a linear buffer of type doubletype. 
 * This object needs to answer the API defined by the virtual functions in IssmAbsVec, 
 * and the contructors required by IssmVec (see IssmVec.h)
 */ 

#ifndef _ISSM_SEQ_VEC_H_
#define _ISSM_SEQ_VEC_H_

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

#include "../../shared/shared.h"
#include <math.h>

/*}}}*/

/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. 
  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 IssmAbsVec;

template <class doubletype> 
class IssmSeqVec: public IssmAbsVec<doubletype>{

	public:

		doubletype* vector;
		int M;

		/*IssmSeqVec constructors, destructors*/
		IssmSeqVec(){/*{{{*/

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

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

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

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

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

		/*IssmSeqVec specific routines*/
		void Echo(void){/*{{{*/

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

			/*do nothing*/

		}
		/*}}}*/
		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;
			}

		}
		/*}}}*/
		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;
			}
		}
		/*}}}*/
		void GetValue(doubletype* pvalue,int dof){/*{{{*/

			*pvalue=this->vector[dof];

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

			*pM=this->M;

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

			*pM=this->M;

		}
		/*}}}*/
		IssmSeqVec<doubletype>* Duplicate(void){/*{{{*/

			return new IssmSeqVec<doubletype>(this->vector,this->M);

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

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

		}
		/*}}}*/
		void AXPY(IssmAbsVec<doubletype>* Xin, doubletype a){/*{{{*/

			int i;

			/*Assume X is of the correct type, and downcast: */
			IssmSeqVec* X=NULL;

			X=(IssmSeqVec<doubletype>*)Xin;

			/*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];

		}
		/*}}}*/
		void AYPX(IssmAbsVec<doubletype>* Xin, doubletype a){/*{{{*/

			int i;

			/*Assume X is of the correct type, and downcast: */
			IssmSeqVec* X=NULL;

			X=(IssmSeqVec<doubletype>*)Xin;

			/*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];

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

			doubletype* buffer=NULL;

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

		}
		/*}}}*/
		void Copy(IssmAbsVec<doubletype>* toin){/*{{{*/

			int i;

			/*Assume toin is of the correct type, and downcast: */
			IssmSeqVec* to=NULL;

			to=(IssmSeqVec<doubletype>*)toin;

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

		}
		/*}}}*/
		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]));
					norm=0; for(i=0;i<this->M;i++)norm=max(norm,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;
			}
		}
		/*}}}*/
		void Scale(doubletype scale_factor){/*{{{*/

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

		}
		/*}}}*/
		doubletype Dot(IssmAbsVec<doubletype>* inputin){/*{{{*/

			int i;

			/*Assume inputin is of the correct type, and downcast: */
			IssmSeqVec* input=NULL;

			input=(IssmSeqVec<doubletype>*)inputin;

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

		}
		/*}}}*/
		void PointwiseDivide(IssmAbsVec<doubletype>* xin,IssmAbsVec<doubletype>* yin){/*{{{*/

			int i;

			/*Assume xin and yin are of the correct type, and downcast: */
			IssmSeqVec* x=NULL;
			IssmSeqVec* y=NULL;

			x=(IssmSeqVec<doubletype>*)xin;
			y=(IssmSeqVec<doubletype>*)yin;

			/*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 _ISSM_SEQ_VEC_H_
