Index: /issm/trunk-jpl/src/c/toolkits/double/MatlabVectorToDoubleVector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/double/MatlabVectorToDoubleVector.cpp	(revision 11732)
+++ /issm/trunk-jpl/src/c/toolkits/double/MatlabVectorToDoubleVector.cpp	(revision 11733)
@@ -80,6 +80,11 @@
 
 		/*allocate and memcpy*/
-		vector=(double*)xmalloc(rows*sizeof(double));
-		memcpy(vector,mxvector_ptr,rows*sizeof(double));
+		if(rows){
+			vector=(double*)xmalloc(rows*sizeof(double));
+			memcpy(vector,mxvector_ptr,rows*sizeof(double));
+		}
+		else{
+			vector=NULL;
+		}
 	}
 
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.cpp	(revision 11733)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.cpp	(revision 11733)
@@ -0,0 +1,234 @@
+/*!\file SeqMat.cpp
+ * \brief: implementation of the SeqMat 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 "../toolkits.h"
+#include "../../shared/shared.h"
+
+/*}}}*/
+
+/*SeqMat constructors and destructor*/
+/*FUNCTION SeqMat::SeqMat(){{{1*/
+SeqMat::SeqMat(){
+
+	this->M=0;
+	this->N=0;
+	this->matrix=NULL;
+}
+/*}}}*/
+/*FUNCTION SeqMat::SeqMat(int M,int N){{{1*/
+SeqMat::SeqMat(int pM,int pN){
+
+	this->M=pM;
+	this->N=pN;
+	this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
+}
+/*}}}*/
+/*FUNCTION SeqMat::SeqMat(int M,int N, double sparsity){{{1*/
+SeqMat::SeqMat(int pM,int pN, double sparsity){
+
+	this->M=pM;
+	this->N=pN;
+	this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
+}
+/*}}}*/
+/*FUNCTION SeqMat(double* serial_mat,int M,int N,double sparsity){{{1*/
+SeqMat::SeqMat(double* serial_mat,int pM,int pN,double sparsity){
+
+	int i,j;
+
+	this->M=pM;
+	this->N=pN;
+	this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
+	memcpy(this->matrix,serial_mat,pM*pN*sizeof(double));
+
+}
+/*}}}*/
+/*FUNCTION SeqMat::SeqMat(int M,int N, int connectivity, int numberofdofspernode){{{1*/
+SeqMat::SeqMat(int pM,int pN, int connectivity,int numberofdofspernode){
+
+	this->M=pM;
+	this->N=pN;
+	this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
+}
+/*}}}*/
+/*FUNCTION SeqMat::~SeqMat(){{{1*/
+SeqMat::~SeqMat(){
+
+	xfree((void**)&this->matrix);
+	M=0;
+	N=0;
+}
+/*}}}*/
+
+/*SeqMat specific routines: */
+/*FUNCTION SeqMat::Echo{{{1*/
+void SeqMat::Echo(void){
+
+	int i,j;
+	printf("SeqMat size %i-%i\n",this->M,this->N);
+	for(i=0;i<M;i++){
+		for(j=0;j<N;j++){
+			printf("%g ",this->matrix[N*i+j]);
+		}
+		printf("\n");
+	}
+}
+/*}}}*/
+
+#ifdef _SERIAL_
+/*FUNCTION SeqMat::ToMatlabMatrix{{{1*/
+mxArray* SeqMat::ToMatlabMatrix(void){
+
+	double* buffer=NULL;
+	
+	mxArray* dataref=NULL;
+
+	/*copy vector into a new buffer: */
+	buffer=(double*)xmalloc(this->M*this->N*sizeof(double));
+	memcpy(buffer,this->matrix,M*N*sizeof(double));
+
+	dataref = mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(dataref,this->M);
+	mxSetN(dataref,this->N);
+	mxSetPr(dataref,buffer);	
+
+
+	/*do not erase buffer!: */
+	return dataref;
+
+}
+/*}}}*/
+/*FUNCTION MatlabMatrixToSeqMat{{{1*/
+SeqMat* MatlabMatrixToSeqMat(const mxArray* dataref){
+
+	SeqMat* output=NULL;
+
+	output=new SeqMat();
+	MatlabMatrixToDoubleMatrix(&output->matrix,&output->M,&output->N,dataref);
+	return output;
+
+}
+/*}}}*/
+#endif
+/*FUNCTION SeqMat::Assemble{{{1*/
+void SeqMat::Assemble(void){
+		
+	/*do nothing*/
+
+}
+/*}}}*/
+/*FUNCTION SeqMat::Norm{{{1*/
+double SeqMat::Norm(NormMode mode){
+
+	double norm;
+	double 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 SeqMat::GetSize{{{1*/
+void SeqMat::GetSize(int* pM,int* pN){
+
+	*pM=this->M;
+	*pN=this->N;
+
+}
+/*}}}*/
+/*FUNCTION SeqMat::GetLocalSize{{{1*/
+void SeqMat::GetLocalSize(int* pM,int* pN){
+	
+	*pM=this->M;
+	*pN=this->N;
+
+}
+/*}}}*/
+/*FUNCTION SeqMat::MatMult{{{1*/
+void SeqMat::MatMult(SeqVec* X,SeqVec* AX){
+
+	int i,j;
+	int XM,AXM;
+	double 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 SeqMat::Duplicate{{{1*/
+SeqMat* SeqMat::Duplicate(void){
+
+	double dummy=0;
+
+	return new SeqMat(this->matrix,this->M,this->N,dummy);
+
+}
+/*}}}*/
+/*FUNCTION SeqMat::ToSerial{{{1*/
+double* SeqMat::ToSerial(void){
+
+	double* buffer=NULL;
+
+	if(this->M*this->N){
+		buffer=(double*)xmalloc(this->M*this->N*sizeof(double));
+		memcpy(buffer,this->matrix,this->M*this->N*sizeof(double));
+	}
+	return buffer;
+
+}
+/*}}}*/
+/*FUNCTION SeqMat::SetValues{{{1*/
+void SeqMat::SetValues(int m,int* idxm,int n,int* idxn,double* 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;
+	}
+
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.h	(revision 11733)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqMat.h	(revision 11733)
@@ -0,0 +1,62 @@
+/*!\file:  SeqMat.h
+ * \brief wrapper to SeqMat objects, which are just wrappers to a simple double* buffer.
+ */ 
+
+#ifndef _SEQMAT_H_
+#define _SEQMAT_H_
+
+/*Headers:*/
+/*{{{1*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../toolkitsenums.h"
+
+#ifdef _SERIAL_
+#include "mex.h"
+#endif
+
+/*}}}*/
+class SeqVec;
+
+class SeqMat{
+
+	public:
+	
+		int M,N; 
+		double* matrix; 
+
+		/*SeqMat constructors, destructors {{{1*/
+		SeqMat();
+		SeqMat(int M,int N);
+		SeqMat(int M,int N,double sparsity);
+		SeqMat(double* serial_mat,int M,int N,double sparsity);
+		SeqMat(int M,int N,int connectivity,int numberofdofspernode);
+		~SeqMat();
+		/*}}}*/
+		/*SeqMat specific routines {{{1*/
+		void Echo(void);
+		#ifdef _SERIAL_
+		mxArray* ToMatlabMatrix(void);
+		#endif
+		void Assemble(void);
+		double Norm(NormMode norm_type);
+		void GetSize(int* pM,int* pN);
+		void GetLocalSize(int* pM,int* pN);
+		void MatMult(SeqVec* X,SeqVec* AX);
+		SeqMat* Duplicate(void);
+		double* ToSerial(void);
+		void SetValues(int m,int* idxm,int n,int* idxn,double* values,InsMode mode);
+		/*}}}*/
+
+};
+		
+/*API :*/
+#ifdef _SERIAL_
+SeqMat*  MatlabMatrixToSeqMat(const mxArray* dataref);
+#endif
+
+#endif //#ifndef _SEQMAT_H_
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.cpp	(revision 11733)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.cpp	(revision 11733)
@@ -0,0 +1,269 @@
+/*!\file SeqVec.cpp
+ * \brief: implementation of the SeqVec 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 "../toolkits.h"
+#include "../../shared/shared.h"
+
+/*}}}*/
+
+/*SeqVec constructors and destructor*/
+/*FUNCTION SeqVec::SeqVec(){{{1*/
+SeqVec::SeqVec(){
+
+	this->M=0;
+	this->vector=NULL;
+}
+/*}}}*/
+/*FUNCTION SeqVec::SeqVec(int M,bool fromlocalsize){{{1*/
+SeqVec::SeqVec(int pM,bool fromlocalsize){
+
+	this->M=pM;
+	if(this->M) this->vector=(double*)xcalloc(pM,sizeof(double));
+	else this->vector=NULL;
+}
+/*}}}*/
+/*FUNCTION SeqVec::SeqVec(double* serial_vec,int M){{{1*/
+SeqVec::SeqVec(double* buffer,int pM){
+
+	int i,j;
+
+	this->M=pM;
+	this->vector=(double*)xcalloc(pM,sizeof(double));
+	memcpy(this->vector,buffer,pM*sizeof(double));
+}
+/*}}}*/
+		/*FUNCTION SeqVec::~SeqVec(){{{1*/
+SeqVec::~SeqVec(){
+	xfree((void**)&this->vector);
+	M=0;
+}
+/*}}}*/
+
+/*SeqVec specific routines: */
+/*FUNCTION SeqVec::Echo{{{1*/
+void SeqVec::Echo(void){
+
+	int i;
+	printf("SeqVec size %i\n",this->M);
+	for(i=0;i<M;i++){
+		printf("%g\n ",vector[i]);
+	}
+}
+/*}}}*/
+
+#ifdef _SERIAL_
+/*FUNCTION SeqVec::ToMatlabVector{{{1*/
+mxArray* SeqVec::ToMatlabVector(void){
+
+	double* buffer=NULL;
+	
+	mxArray* dataref=NULL;
+
+	/*copy vector into a new buffer: */
+	buffer=(double*)xmalloc(this->M*sizeof(double));
+	memcpy(buffer,vector,M*sizeof(double));
+
+	dataref = mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(dataref,this->M);
+	mxSetN(dataref,1);
+	mxSetPr(dataref,buffer);	
+
+
+	/*do not erase buffer!: */
+	return dataref;
+
+}
+/*}}}*/
+/*FUNCTION MatlabVectorToSeqVec{{{1*/
+SeqVec* MatlabVectorToSeqVec(const mxArray* dataref){
+
+	SeqVec* output=NULL;
+
+	output=new SeqVec();
+	MatlabVectorToDoubleVector(&output->vector,&output->M,dataref);
+	return output;
+
+}
+/*}}}*/
+#endif
+/*FUNCTION SeqVec::Assemble{{{1*/
+void SeqVec::Assemble(void){
+		
+	/*do nothing*/
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::SetValues{{{1*/
+void SeqVec::SetValues(int ssize, int* list, double* 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 SeqVec::SetValue{{{1*/
+void SeqVec::SetValue(int dof, double 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 SeqVec::GetValue{{{1*/
+void SeqVec::GetValue(double* pvalue,int dof){
+
+	*pvalue=this->vector[dof];
+
+}
+/*}}}*/
+		
+/*FUNCTION SeqVec::GetSize{{{1*/
+void SeqVec::GetSize(int* pM){
+
+	*pM=this->M;
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::GetLocalSize{{{1*/
+void SeqVec::GetLocalSize(int* pM){
+	
+	*pM=this->M;
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::Duplicate{{{1*/
+SeqVec* SeqVec::Duplicate(void){
+	
+	return new SeqVec(this->vector,this->M);
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::Set{{{1*/
+void SeqVec::Set(double value){
+
+	int i;
+	for(i=0;i<this->M;i++)this->vector[i]=value;
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::AXPY{{{1*/
+void SeqVec::AXPY(SeqVec* X, double 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 SeqVec::AYPX{{{1*/
+void SeqVec::AYPX(SeqVec* X, double 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 SeqVec::ToMPISerial{{{1*/
+double* SeqVec::ToMPISerial(void){
+
+	double* buffer=NULL;
+
+	if(this->M){
+		buffer=(double*)xmalloc(this->M*sizeof(double));
+		memcpy(buffer,this->vector,this->M*sizeof(double));
+	}
+	return buffer;
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::Copy{{{1*/
+void SeqVec::Copy(SeqVec* to){
+
+	int i;
+
+	to->M=this->M;
+	for(i=0;i<this->M;i++)to->vector[i]=this->vector[i];
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::Norm{{{1*/
+double SeqVec::Norm(NormMode mode){
+
+	double norm;
+	int i;
+
+	switch(mode){
+		case NORM_INF:
+			norm=0; for(i=0;i<this->M;i++)norm=fmax(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 SeqVec::Scale{{{1*/
+void SeqVec::Scale(double scale_factor){
+
+	int i;
+	for(i=0;i<this->M;i++)this->vector[i]=scale_factor*this->vector[i];
+	
+}
+/*}}}*/
+/*FUNCTION SeqVec::Dot{{{1*/
+double SeqVec::Dot(SeqVec* input){
+
+	int i;
+
+	double dot=0;
+	for(i=0;i<this->M;i++)dot+=this->vector[i]*input->vector[i];
+	return dot;
+
+}
+/*}}}*/
+/*FUNCTION SeqVec::PointwiseDivide{{{1*/
+void SeqVec::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];
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.h	(revision 11733)
+++ /issm/trunk-jpl/src/c/toolkits/issm/SeqVec.h	(revision 11733)
@@ -0,0 +1,67 @@
+/*!\file:  SeqVec.h
+ * \brief wrapper to our SeqVec object, which is just a wrapper to a double* 
+ */ 
+
+#ifndef _SEQVEC_H_
+#define _SEQVEC_H_
+
+/*Headers:*/
+/*{{{1*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../toolkitsenums.h"
+
+#ifdef _SERIAL_
+#include "mex.h"
+#endif
+
+/*}}}*/
+
+class SeqVec{
+
+	public:
+	
+		double* vector;
+		int M;
+
+		/*SeqVec constructors, destructors {{{1*/
+		SeqVec();
+		SeqVec(int M,bool fromlocalsize=false);
+		SeqVec(double* buffer, int M);
+		~SeqVec();
+		/*}}}*/
+		/*SeqVec specific routines {{{1*/
+		void Echo(void);
+		#ifdef _SERIAL_
+		mxArray* ToMatlabVector(void);
+		#endif
+		void Assemble(void);
+		void SetValues(int ssize, int* list, double* values, InsMode mode);
+		void SetValue(int dof, double value, InsMode  mode);
+		void GetValue(double* pvalue, int dof);
+		void GetSize(int* pM);
+		void GetLocalSize(int* pM);
+		SeqVec* Duplicate(void);
+		void Set(double value);
+		void AXPY(SeqVec* X, double a);
+		void AYPX(SeqVec* X, double a);
+		double* ToMPISerial(void);
+		void Copy(SeqVec* to);
+		double Norm(NormMode norm_type);
+		void Scale(double scale_factor);
+		void PointwiseDivide(SeqVec* x,SeqVec* y);
+		double Dot(SeqVec* vector);
+		/*}}}*/
+};
+
+
+/*API :*/
+#ifdef _SERIAL_
+SeqVec*  MatlabVectorToSeqVec(const mxArray* dataref);
+#endif
+
+#endif //#ifndef _SEQVEC_H_
Index: /issm/trunk-jpl/src/c/toolkits/issm/issmtoolkit.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/issmtoolkit.h	(revision 11733)
+++ /issm/trunk-jpl/src/c/toolkits/issm/issmtoolkit.h	(revision 11733)
@@ -0,0 +1,11 @@
+/* \file issmtoolkit.h
+ * \brief all includes for MPI layer
+ */
+
+#ifndef _ISSM_TOOLKIT_H_
+#define _ISSM_TOOLKIT_H_
+
+#include "./SeqMat.h"
+#include "./SeqVec.h"
+
+#endif
Index: /issm/trunk-jpl/src/c/toolkits/toolkits.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/toolkits.h	(revision 11732)
+++ /issm/trunk-jpl/src/c/toolkits/toolkits.h	(revision 11733)
@@ -22,4 +22,5 @@
 #include "./double/double.h"
 #include "./toolkitsenums.h"
+#include "./issm/issmtoolkit.h"
 
 #endif