/*!\file:  SeqMat.h
 * \brief wrapper to SeqMat objects, which are just wrappers to a simple IssmDouble or IssmPDouble* buffer.
 */ 

#ifndef _SEQMAT_H_
#define _SEQMAT_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"
#include "./SeqVec.h"

/*}}}*/

/*We need to template this class, in case we want to create Matrices that hold IssmDouble* matrix or IssmPDouble* matrix. 
  Such matrices 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 SeqMat{

	public:

		int M,N; 
		doubletype* matrix;  /*here, doubletype is either IssmDouble or IssmPDouble*/

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

			this->M=0;
			this->N=0;
			this->matrix=NULL;
		}
		/*}}}*/
		/*FUNCTION SeqMat(int M,int N){{{*/
		SeqMat(int pM,int pN){

			this->M=pM;
			this->N=pN;
			this->matrix=NULL;
			if(M*N) this->matrix=xNewZeroInit<doubletype>(pM*pN);
		}
		/*}}}*/
		/*FUNCTION SeqMat(int M,int N, doubletype sparsity){{{*/
		SeqMat(int pM,int pN, doubletype sparsity){

			this->M=pM;
			this->N=pN;
			this->matrix=NULL;
			if(M*N) this->matrix=xNewZeroInit<doubletype>(pM*pN);
		}
		/*}}}*/
		/*FUNCTION SeqMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
		SeqMat(int m,int n,int pM,int pN,int* d_nnz,int* o_nnz){

			this->M=pM;
			this->N=pN;
			this->matrix=NULL;
			if(pM*pN) this->matrix=xNewZeroInit<doubletype>(pM*pN);
		}
		/*}}}*/
		/*FUNCTION SeqMat(doubletype* serial_mat,int M,int N,doubletype sparsity){{{*/
		SeqMat(doubletype* serial_mat,int pM,int pN,doubletype sparsity){

			this->M=pM;
			this->N=pN;
			this->matrix=NULL;
			if(M*N){
				this->matrix=xNewZeroInit<doubletype>(pM*pN);
				xMemCpy<doubletype>(this->matrix,serial_mat,pM*pN);
			}

		}
		/*}}}*/
		/*FUNCTION SeqMat(int M,int N, int connectivity, int numberofdofspernode){{{*/
		SeqMat(int pM,int pN, int connectivity,int numberofdofspernode){

			this->M=pM;
			this->N=pN;
			this->matrix=NULL;
			if(M*N) this->matrix=xNewZeroInit<doubletype>(pM*pN);
		}
		/*}}}*/
		/*FUNCTION ~SeqMat(){{{*/
		~SeqMat(){

			xDelete<doubletype>(this->matrix);
			M=0;
			N=0;
		}
		/*}}}*/

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

			int i,j;
			_printLine_("SeqMat size " << this->M << "-" << this->N);
			for(i=0;i<M;i++){
				for(j=0;j<N;j++){
					_printString_(this->matrix[N*i+j] << " ");
				}
				_printLine_("");
			}
		}
		/*}}}*/
		/*FUNCTION Assemble{{{*/
		void Assemble(void){

			/*do nothing*/

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

			doubletype norm;
			doubletype absolute;
			int i,j;

			switch(mode){
				case NORM_INF:
					norm=0;
					for(i=0;i<this->M;i++){
						absolute=0;
						for(j=0;j<this->N;j++){
							absolute+=fabs(this->matrix[N*i+j]);
						}
						norm=max(norm,absolute);
					}
					return norm;
					break; 
				default:
					_error_("unknown norm !");
					break;
			}
		}
		/*}}}*/
		/*FUNCTION GetSize{{{*/
		void GetSize(int* pM,int* pN){

			*pM=this->M;
			*pN=this->N;

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

			*pM=this->M;
			*pN=this->N;

		}
		/*}}}*/
		/*FUNCTION MatMult{{{*/
		void MatMult(SeqVec<doubletype>* X,SeqVec<doubletype>* AX){

			int i,j;
			int XM,AXM;
			doubletype dummy;

			X->GetSize(&XM);
			AX->GetSize(&AXM);

			if(M!=AXM)_error_("A and AX should have the same number of rows!");
			if(N!=XM)_error_("A and X should have the same number of columns!");

			for(i=0;i<M;i++){
				dummy=0;
				for(j=0;j<N;j++){
					dummy+= this->matrix[N*i+j]*X->vector[j];
				}
				AX->vector[i]=dummy;
			}

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

			doubletype dummy=0;

			return new SeqMat(this->matrix,this->M,this->N,dummy);

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

			doubletype* buffer=NULL;

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

		}
		/*}}}*/
		/*FUNCTION SetValues{{{*/
		void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode){

			int i,j;
			switch(mode){
				case ADD_VAL:
					for(i=0;i<m;i++) for(j=0;j<n;j++) this->matrix[N*idxm[i]+idxn[j]]+=values[n*i+j];
					break;
				case INS_VAL:
					for(i=0;i<m;i++) for(j=0;j<n;j++) this->matrix[N*idxm[i]+idxn[j]]=values[n*i+j];
					break;
				default:
					_error_("unknown insert mode!");
					break;
			}

		}
		/*}}}*/
		/*FUNCTION Convert{{{*/
		void Convert(MatrixType type){

			/*do nothing*/

		}
		/*}}}*/		

};

#endif //#ifndef _SEQMAT_H_
