Index: /issm/trunk-jpl/src/wrappers/matlab/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/Makefile.am	(revision 15519)
+++ /issm/trunk-jpl/src/wrappers/matlab/Makefile.am	(revision 15520)
@@ -13,22 +13,7 @@
 io_sources= ./include/matlabincludes.h\
 				./io/matlabio.h\
-				./io/MatlabNArrayToNArray.cpp\
 				./io/CheckNumMatlabArguments.cpp\
-				./io/mxGetAssignedField.cpp\
 				./io/WriteMatlabData.cpp\
-				./io/FetchMatlabData.cpp\
-				./io/OptionParse.cpp\
-				./io/MatlabMatrixToMatrix.cpp\
-				./io/MatlabVectorToVector.cpp\
-				./io/MatlabVectorToDoubleVector.cpp\
-				./io/MatlabMatrixToDoubleMatrix.cpp\
-				./io/MatlabMatrixToIssmMat.cpp\
-				./io/MatlabVectorToIssmVec.cpp
-
-				
-if PETSC
-io_sources += ./io/MatlabMatrixToPetscMat.cpp\
-				./io/MatlabVectorToPetscVec.cpp
-endif
+				./io/FetchMatlabData.cpp
 
 ALLCXXFLAGS= -fPIC -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_ -D_WRAPPERS_ $(CXXFLAGS) $(CXXOPTFLAGS) 
@@ -160,5 +145,5 @@
 Chaco_la_SOURCES = ../Chaco/Chaco.cpp\
 						 ../Chaco/Chaco.h
-Chaco_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(CHACOLIB) $(GSLLIB)
+Chaco_la_LIBADD = ${deps} $(MPILIB) $(CHACOLIB) $(GSLLIB) $(PETSCLIB)
 
 ContourToMesh_la_SOURCES = ../ContourToMesh/ContourToMesh.cpp\
Index: /issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp	(revision 15519)
+++ /issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp	(revision 15520)
@@ -10,4 +10,5 @@
 
 #include "./matlabio.h"
+#include <cstring> 
 
 /*Primitive data types*/
@@ -484,49 +485,4 @@
 
 /*ISSM objects*/
-/*FUNCTION FetchData(Matrix<double>** pmatrix,const mxArray* dataref){{{*/
-void FetchData(Matrix<double>** pmatrix,const mxArray* dataref){
-
-	Matrix<double>* outmatrix=NULL;
-	int dummy=0;
-
-	if (mxIsClass(dataref,"double") ){
-
-		/*Convert matlab matrix to matrix: */
-		outmatrix=MatlabMatrixToMatrix(dataref);
-
-	}
-	else{
-		/*This is an error: we don't have the correct input!: */
-		_error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet");
-	}
-
-	/*Assign output pointers:*/
-	*pmatrix=outmatrix;
-}
-/*}}}*/
-/*FUNCTION FetchData(Vector<double>** pvector,const mxArray* dataref){{{*/
-void FetchData(Vector<double>** pvector,const mxArray* dataref){
-
-	Vector<double>* vector=NULL;
-	int dummy;
-
-	if(mxIsEmpty(dataref)){
-		/*Nothing to pick up. Just initialize matrix pointer to NULL: */
-		vector=new Vector<double>(0);
-	}
-	else if (mxIsClass(dataref,"double") ){
-
-		/*Convert matlab vector to petsc vector: */
-		vector=MatlabVectorToVector(dataref);
-	}
-	else{
-		/*This is an error: we don't have the correct input!: */
-		_error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet");
-	}
-
-	/*Assign output pointers:*/
-	*pvector=vector;
-}
-/*}}}*/
 /*FUNCTION FetchData(BamgGeom** pbamggeom,const mxArray* dataref){{{*/
 void FetchData(BamgGeom** pbamggeom,const mxArray* dataref){
@@ -684,2 +640,574 @@
 }
 /*}}}*/
+
+/*Toolkit*/
+/*FUNCTION MatlabMatrixToDoubleMatrix {{{*/
+int MatlabMatrixToDoubleMatrix(double** pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){
+
+	int        i,j,count,rows,cols;
+
+	/*output: */
+	double* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex*    ir=NULL;
+	mwIndex*    jc=NULL;
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		double* pmxmatrix=(double*)mxGetPr(mxmatrix);
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+
+		if(rows*cols){
+			matrix=xNewZeroInit<double>(rows*cols);
+
+			/*Now, get ir,jc and pr: */
+			ir=mxGetIr(mxmatrix);
+			jc=mxGetJc(mxmatrix);
+
+			/*Now, start inserting data into double* matrix: */
+			count=0;
+			for(i=0;i<cols;i++){
+				for(j=0;j<(jc[i+1]-jc[i]);j++){
+					matrix[rows*ir[count]+i]=pmxmatrix[count];
+					count++;
+				}
+			}
+		}
+
+	}
+	else if(mxIsClass(mxmatrix,"double")){
+		/*Dealing with dense matrix: recover pointer and size: */
+		double* pmxmatrix=(double*)mxGetPr(mxmatrix);
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+
+		/*Create serial matrix: */
+		if(rows*cols){
+			matrix=xNewZeroInit<double>(rows*cols);
+
+			for(i=0;i<rows;i++){
+				for(j=0;j<cols;j++){
+					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
+				}
+			}
+		}
+	}
+	else if(mxIsClass(mxmatrix,"single")){
+		/*Dealing with dense matrix: recover pointer and size: */
+		float *pmxmatrix=(float*)mxGetPr(mxmatrix);
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+
+		/*Create serial matrix: */
+		if(rows*cols){
+			matrix=xNewZeroInit<double>(rows*cols);
+
+			for(i=0;i<rows;i++){
+				for(j=0;j<cols;j++){
+					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
+				}
+			}
+		}
+	}
+	else if(mxIsClass(mxmatrix,"int16")){
+		/*Dealing with dense matrix: recover pointer and size: */
+		short int *pmxmatrix=(short*)mxGetPr(mxmatrix);
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+
+		/*Create serial matrix: */
+		if(rows*cols){
+			matrix=xNewZeroInit<double>(rows*cols);
+
+			for(i=0;i<rows;i++){
+				for(j=0;j<cols;j++){
+					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
+				}
+			}
+		}
+	}
+	else if(mxIsClass(mxmatrix,"uint8")){
+		/*Dealing with dense matrix: recover pointer and size: */
+		char *pmxmatrix=(char*)mxGetPr(mxmatrix);
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+
+		/*Create serial matrix: */
+		if(rows*cols){
+			matrix=xNewZeroInit<double>(rows*cols);
+
+			for(i=0;i<rows;i++){
+				for(j=0;j<cols;j++){
+					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
+				}
+			}
+		}
+	}
+	else{
+		_error_("Matlab matrix type Not implemented yet");
+	}
+
+	/*Assign output pointer: */
+	*pmatrix=matrix;
+	*pmatrix_rows=rows;
+	*pmatrix_cols=cols;
+
+	return 1;
+}/*}}}*/
+/*FUNCTION MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{*/
+int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
+
+	int  i,j,rows,cols;
+	int  numel,ndims;
+	int *size,*dims;
+	double* mxmatrix_ptr=NULL;
+	const mwSize* ipt=NULL;
+
+	/*output: */
+	double* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex *ir    = NULL;
+	mwIndex *jc    = NULL;
+	double  *pr    = NULL;
+	int      count;
+
+	/*get Matlab matrix information: */
+	numel=mxGetNumberOfElements(mxmatrix);
+	ndims=mxGetNumberOfDimensions(mxmatrix);
+	ipt  =mxGetDimensions(mxmatrix);
+	size =xNew<int>(ndims);
+	for (i=0;i<ndims;i++) size[i]=(int)ipt[i];
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		rows = mxGetM(mxmatrix);
+		cols = mxGetN(mxmatrix);
+
+		matrix=xNewZeroInit<double>(rows*cols);
+
+		/*Now, get ir,jc and pr: */
+		ir = mxGetIr(mxmatrix);
+		jc = mxGetJc(mxmatrix);
+		pr = mxGetPr(mxmatrix);
+
+		/*Now, start inserting data into double* matrix: */
+		count=0;
+		for(i=0;i<cols;i++){
+			for(j=0;j<(jc[i+1]-jc[i]);j++){
+				*(matrix+rows*ir[count]+i)=pr[count];
+				count++;
+			}
+		}
+
+	}
+	else{
+
+		/*Dealing with dense matrix: recover pointer and size: */
+		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
+
+		/*Create serial matrix: */
+		matrix=xNewZeroInit<double>(numel);
+
+		dims=xNew<int>(ndims);
+		for(i=0;i<numel;i++){
+			ColumnWiseDimsFromIndex(dims,i,size,ndims);
+			j = IndexFromRowWiseDims(dims,size,ndims);
+			matrix[j]=(double)mxmatrix_ptr[i];
+		}
+		xDelete<int>(dims);
+	}
+
+	/*Assign output pointer: */
+	*pmatrix       = matrix;
+	*pmatrix_numel = numel;
+	*pmatrix_ndims = ndims;
+	*pmatrix_size  = size;
+
+	return 1;
+}
+/*}}}*/
+/*FUNCTION MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{*/
+int MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
+
+	int  i,j,rows,cols;
+	int  numel,ndims;
+	int *size,*dims;
+	bool* mxmatrix_ptr=NULL;
+	const mwSize* ipt=NULL;
+
+	/*output: */
+	bool* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex *ir    = NULL;
+	mwIndex *jc    = NULL;
+	bool    *pm    = NULL;
+	int      count;
+
+	/*get Matlab matrix information: */
+	numel = mxGetNumberOfElements(mxmatrix);
+	ndims = mxGetNumberOfDimensions(mxmatrix);
+	ipt   = mxGetDimensions(mxmatrix);
+	size  = xNew<int>(ndims);
+	for (i=0;i<ndims;i++) size[i]=(int)ipt[i];
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		rows=mxGetM(mxmatrix);
+		cols=mxGetN(mxmatrix);
+		matrix=xNewZeroInit<bool>(rows*cols);
+
+		/*Now, get ir,jc and pm: */
+		ir=mxGetIr(mxmatrix);
+		jc=mxGetJc(mxmatrix);
+		pm=(bool*)mxGetData(mxmatrix);
+
+		/*Now, start inserting data into bool* matrix: */
+		count=0;
+		for(i=0;i<cols;i++){
+			for(j=0;j<(jc[i+1]-jc[i]);j++){
+				matrix[rows*ir[count]+i]=pm[count];
+				count++;
+			}
+		}
+	}
+	else{
+
+		/*Dealing with dense matrix: recover pointer and size: */
+		mxmatrix_ptr=(bool*)mxGetData(mxmatrix);
+
+		/*Create serial matrix: */
+		matrix=xNew<bool>(numel);
+		dims=xNew<int>(ndims);
+		for(i=0;i<numel;i++){
+			ColumnWiseDimsFromIndex(dims,i,size,ndims);
+			j=IndexFromRowWiseDims(dims,size,ndims);
+			matrix[j]=(bool)mxmatrix_ptr[i];
+		}
+		xDelete<int>(dims);
+	}
+
+	/*Assign output pointer: */
+	*pmatrix       = matrix;
+	*pmatrix_numel = numel;
+	*pmatrix_ndims = ndims;
+	*pmatrix_size  = size;
+
+	return 1;
+}
+/*}}}*/
+/*FUNCTION MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{*/
+int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
+
+	int           i,j,rows,cols;
+	int           numel,ndims;
+	int          *size , *dims;
+	mxChar       *mxmatrix_ptr = NULL;
+	const mwSize *ipt          = NULL;
+
+	/*output: */
+	char* matrix=NULL;
+
+	/*matlab indices: */
+	mwIndex *ir    = NULL;
+	mwIndex *jc    = NULL;
+	char    *pm    = NULL;
+	int      count;
+
+	/*get Matlab matrix information: */
+	numel = mxGetNumberOfElements(mxmatrix);
+	ndims = mxGetNumberOfDimensions(mxmatrix);
+	ipt   = mxGetDimensions(mxmatrix);
+	size  = xNew<int>(ndims);
+	for (i=0;i<ndims;i++) size[i]=(int)ipt[i];
+
+	/*Ok, first check if we are dealing with a sparse or full matrix: */
+	if (mxIsSparse(mxmatrix)){
+
+		/*Dealing with sparse matrix: recover size first: */
+		rows = mxGetM(mxmatrix);
+		cols = mxGetN(mxmatrix);
+		matrix=xNew<char>(rows*cols);
+
+		/*Now, get ir,jc and pm: */
+		ir = mxGetIr(mxmatrix);
+		jc = mxGetJc(mxmatrix);
+		pm = (char*)mxGetData(mxmatrix);
+
+		/*Now, start inserting data into char* matrix: */
+		count=0;
+		for(i=0;i<cols;i++){
+			for(j=0;j<(jc[i+1]-jc[i]);j++){
+				matrix[rows*ir[count]+i]=(char)pm[count];
+				count++;
+			}
+		}
+	}
+	else{
+		/*Dealing with dense matrix: recover pointer and size: */
+		mxmatrix_ptr=mxGetChars(mxmatrix);
+
+		/*Create serial matrix: */
+		matrix=xNew<char>(numel+1);
+		matrix[numel]='\0';
+
+		/*looping code adapted from Matlab example explore.c: */
+		int elements_per_page = size[0] * size[1];
+		/* total_number_of_pages = size[2] x size[3] x ... x size[N-1] */
+		int total_number_of_pages = 1;
+		for (i=2; i<ndims; i++) {
+			total_number_of_pages *= size[i];
+		}
+
+		i=0;
+		for (int page=0; page < total_number_of_pages; page++) {
+			int row;
+			/* On each page, walk through each row. */
+			for (row=0; row<size[0]; row++)  {
+				int column;
+				j = (page * elements_per_page) + row;
+
+				/* Walk along each column in the current row. */
+				for (column=0; column<size[1]; column++) {
+					*(matrix+i++)=(char)*(mxmatrix_ptr+j);
+					j += size[0];
+				}
+			}
+		}
+	}
+
+	/*Assign output pointer: */
+	*pmatrix       = matrix;
+	*pmatrix_numel = numel;
+	*pmatrix_ndims = ndims;
+	*pmatrix_size  = size;
+
+	return 1;
+}
+/*}}}*/
+/*FUNCTION mxGetAssignedField{{{*/
+mxArray* mxGetAssignedField(const mxArray* pmxa_array,int number,const char* field){
+
+	//output
+	mxArray* mxfield=NULL;
+
+	//input
+	mxArray    *inputs[2];
+	mxArray    *pindex      = NULL;
+	const char *fnames[2];
+	mwSize      ndim        = 2;
+	mwSize      onebyone[2] = {1,1};
+
+	//We want to call the subsasgn method, and get the returned array.This ensures that if we are running 
+	//large sized problems, the data is truly loaded from disk by the model subsasgn class method.
+	inputs[0]=(mxArray*)pmxa_array; //this is the model
+
+	//create index structure used in the assignment (index.type='.' and index.subs='x' for field x for ex)
+	fnames[0] = "type";
+	fnames[1] = "subs";
+	pindex=mxCreateStructArray( ndim,onebyone,2,fnames);
+	mxSetField( pindex, 0, "type",mxCreateString("."));
+	mxSetField( pindex, 0, "subs",mxCreateString(field));
+	inputs[1]=pindex;
+
+	mexCallMATLAB( 1, &mxfield, 2, (mxArray**)inputs, "subsref");
+
+	return mxfield;
+}/*}}}*/
+
+GenericOption<double>* OptionDoubleParse( char* name, const mxArray* prhs[]){ /*{{{*/
+
+	GenericOption<double> *odouble = NULL;
+
+	/*check and parse the name  */
+	odouble=new GenericOption<double>();
+	odouble->name =xNew<char>(strlen(name)+1);
+	memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
+	FetchData(&odouble->value,prhs[0]);
+	odouble->numel=1;
+	odouble->ndims=1;
+	odouble->size=NULL;
+
+	return(odouble);
+}/*}}}*/
+GenericOption<double*>* OptionDoubleArrayParse( char* name, const mxArray* prhs[]){ /*{{{*/
+
+	GenericOption<double*> *odouble = NULL;
+
+	/*check and parse the name  */
+	odouble=new GenericOption<double*>();
+	odouble->name =xNew<char>(strlen(name)+1);
+	memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
+
+	/*check and parse the value  */
+	if (!mxIsClass(prhs[0],"double")){
+		_error_("Value of option \"" << odouble->name  << "\" must be class \"double\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
+	}
+	FetchData(&odouble->value,&odouble->numel,&odouble->ndims,&odouble->size,prhs[0]);
+
+	return(odouble);
+}/*}}}*/
+GenericOption<bool*>* OptionLogicalParse( char* name, const mxArray* prhs[]){ /*{{{*/
+
+	GenericOption<bool*> *ological = NULL;
+
+	/*check and parse the name  */
+	ological=new GenericOption<bool*>();
+	ological->name =xNew<char>(strlen(name)+1);
+	memcpy(ological->name,name,(strlen(name)+1)*sizeof(char));
+
+	/*check and parse the value  */
+	if (!mxIsClass(prhs[0],"logical")){
+		_error_("Value of option \"" << ological->name  << "\" must be class \"logical\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
+	}
+	FetchData(&ological->value,&ological->numel,&ological->ndims,&ological->size,prhs[0]);
+
+	return(ological);
+}/*}}}*/
+GenericOption<char*>* OptionCharParse( char* name, const mxArray* prhs[]){ /*{{{*/
+
+	GenericOption<char*>  *ochar = NULL;
+
+	/*check and parse the name  */
+	ochar=new GenericOption<char*>();
+	ochar->name =xNew<char>(strlen(name)+1);
+	memcpy(ochar->name,name,(strlen(name)+1)*sizeof(char));
+
+	/*check and parse the value  */
+	if (!mxIsClass(prhs[0],"char")){
+		_error_("Value of option \"" << ochar->name  << "\" must be class \"char\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
+	}
+	FetchData(&ochar->value,&ochar->numel,&ochar->ndims,&ochar->size,prhs[0]);
+
+	return(ochar);
+}/*}}}*/
+GenericOption<Options**>* OptionStructParse( char* name, const mxArray* prhs[]){ /*{{{*/
+
+	int            i;
+	char           namei[161];
+	Option*                   option      = NULL;
+	GenericOption<Options**>  *ostruct    = NULL;
+	const mwSize  *ipt        = NULL;
+	const mxArray *structi;
+	mwIndex        sindex;
+
+	/*check and parse the name  */
+	ostruct=new GenericOption<Options**>();
+	ostruct->name =xNew<char>(strlen(name)+1);
+	memcpy(ostruct->name,name,(strlen(name)+1)*sizeof(char));
+
+	/*check and parse the value  */
+	if (!mxIsClass(prhs[0],"struct")){
+		_error_("Value of option \"" << ostruct->name  << "\" must be class \"struct\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
+	}
+	ostruct->numel=mxGetNumberOfElements(prhs[0]);
+	ostruct->ndims=mxGetNumberOfDimensions(prhs[0]);
+	ipt           =mxGetDimensions(prhs[0]);
+	ostruct->size =xNew<int>(ostruct->ndims);
+	for (i=0; i<ostruct->ndims; i++) ostruct->size[i]=(int)ipt[i];
+	if (ostruct->numel) ostruct->value=xNew<Options*>(ostruct->numel);
+
+	/*loop through and process each element of the struct array  */
+	for (sindex=0; sindex<ostruct->numel; sindex++) {
+		ostruct->value[sindex]=new Options;
+
+		/*loop through and process each field for the element  */
+		for (i=0; i<mxGetNumberOfFields(prhs[0]); i++) {
+			sprintf(namei,"%s.%s",name,mxGetFieldNameByNumber(prhs[0],i));
+			structi=mxGetFieldByNumber(prhs[0],sindex,i);
+
+			option=(Option*)OptionParse(namei,&structi);
+			ostruct->value[sindex]->AddObject((Object*)option);
+			option=NULL;
+		}
+	}
+
+	return(ostruct);
+}/*}}}*/
+GenericOption<Options*>* OptionCellParse( char* name, const mxArray* prhs[]){ /*{{{*/
+
+	int            i;
+	int           *dims;
+	char           namei[161];
+	char           cstr[81];
+	GenericOption<Options*> *ocell      = NULL;
+	Option        *option     = NULL;
+	const mwSize  *ipt        = NULL;
+	const mxArray *celli;
+	mwIndex        cindex;
+
+	/*check and parse the name  */
+	ocell=new GenericOption<Options*>();
+	ocell->name =xNew<char>(strlen(name)+1);
+	memcpy(ocell->name,name,(strlen(name)+1)*sizeof(char));
+
+	/*check and parse the value  */
+	if (!mxIsClass(prhs[0],"cell")){
+		_error_("Value of option \"" << ocell->name  << "\" must be class \"cell\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
+	}
+
+	ocell->numel=mxGetNumberOfElements(prhs[0]);
+	ocell->ndims=mxGetNumberOfDimensions(prhs[0]);
+	ipt         =mxGetDimensions(prhs[0]);
+	ocell->size =xNew<int>(ocell->ndims);
+	for (i=0; i<ocell->ndims; i++) ocell->size[i]=(int)ipt[i];
+	ocell->value=new Options;
+
+	/*loop through and process each element of the cell array  */
+	dims=xNew<int>(ocell->ndims);
+	for (cindex=0; cindex<ocell->numel; cindex++) {
+		ColumnWiseDimsFromIndex(dims,(int)cindex,ocell->size,ocell->ndims);
+		StringFromDims(cstr,dims,ocell->ndims);
+		#ifdef _INTEL_WIN_
+			_snprintf(namei,161,"%s%s",name,cstr);
+		#else
+			snprintf(namei,161,"%s%s",name,cstr);
+		#endif
+		celli=mxGetCell(prhs[0],cindex);
+
+		option=(Option*)OptionParse(namei,&celli);
+		ocell->value->AddObject((Object*)option);
+		option=NULL;
+	}
+	xDelete<int>(dims);
+
+	return(ocell);
+}/*}}}*/
+Option* OptionParse(char* name, const mxArray* prhs[]){ /*{{{*/
+
+	Option  *option = NULL;
+	mxArray *lhs[1];
+
+	/*parse the value according to the matlab data type  */
+	if     (mxIsClass(prhs[0],"double")  && (mxGetNumberOfElements(prhs[0])==1))
+	 option=(Option*)OptionDoubleParse(name,prhs);
+	else if(mxIsClass(prhs[0],"double")  && (mxGetNumberOfElements(prhs[0])!=1))
+	 option=(Option*)OptionDoubleArrayParse(name,prhs);
+	else if(mxIsClass(prhs[0],"logical"))
+	 option=(Option*)OptionLogicalParse(name,prhs);
+	else if(mxIsClass(prhs[0],"char"))
+	 option=(Option*)OptionCharParse(name,prhs);
+	else if(mxIsClass(prhs[0],"struct"))
+	 option=(Option*)OptionStructParse(name,prhs);
+	else if(mxIsClass(prhs[0],"cell"))
+	 option=(Option*)OptionCellParse(name,prhs);
+	else {
+		_printf0_("  Converting value of option \"" << name << "\" from unrecognized class \"" << mxGetClassName(prhs[0]) << "\" to class \"" << "struct" << "\".\n");
+		if (!mexCallMATLAB(1,lhs,1,(mxArray**)prhs,"struct")) {
+			option=(Option*)OptionStructParse(name,(const mxArray**)lhs);
+			mxDestroyArray(lhs[0]);
+		}
+		else _error_("Second argument value of option \""<< name <<"\" is of unrecognized class \""<< mxGetClassName(prhs[0]) <<"\".");
+	}
+
+	return(option);
+}/*}}}*/
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToDoubleMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToDoubleMatrix.cpp	(revision 15519)
+++ 	(revision )
@@ -1,129 +1,0 @@
-/* \file MatlabMatrixToDoubleMatrix.cpp
- * \brief: convert a sparse or dense matlab matrix to a double* pointer
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-/*Matlab includes: */
-#include "./matlabio.h"
-
-int MatlabMatrixToDoubleMatrix(double** pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){
-
-	int        i,j,count,rows,cols;
-
-	/*output: */
-	double* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex*    ir=NULL;
-	mwIndex*    jc=NULL;
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		double* pmxmatrix=(double*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-
-		if(rows*cols){
-			matrix=xNewZeroInit<double>(rows*cols);
-
-			/*Now, get ir,jc and pr: */
-			ir=mxGetIr(mxmatrix);
-			jc=mxGetJc(mxmatrix);
-
-			/*Now, start inserting data into double* matrix: */
-			count=0;
-			for(i=0;i<cols;i++){
-				for(j=0;j<(jc[i+1]-jc[i]);j++){
-					matrix[rows*ir[count]+i]=pmxmatrix[count];
-					count++;
-				}
-			}
-		}
-
-	}
-	else if(mxIsClass(mxmatrix,"double")){
-		/*Dealing with dense matrix: recover pointer and size: */
-		double* pmxmatrix=(double*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-
-		/*Create serial matrix: */
-		if(rows*cols){
-			matrix=xNewZeroInit<double>(rows*cols);
-
-			for(i=0;i<rows;i++){
-				for(j=0;j<cols;j++){
-					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
-				}
-			}
-		}
-	}
-	else if(mxIsClass(mxmatrix,"single")){
-		/*Dealing with dense matrix: recover pointer and size: */
-		float *pmxmatrix=(float*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-
-		/*Create serial matrix: */
-		if(rows*cols){
-			matrix=xNewZeroInit<double>(rows*cols);
-
-			for(i=0;i<rows;i++){
-				for(j=0;j<cols;j++){
-					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
-				}
-			}
-		}
-	}
-	else if(mxIsClass(mxmatrix,"int16")){
-		/*Dealing with dense matrix: recover pointer and size: */
-		short int *pmxmatrix=(short*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-
-		/*Create serial matrix: */
-		if(rows*cols){
-			matrix=xNewZeroInit<double>(rows*cols);
-
-			for(i=0;i<rows;i++){
-				for(j=0;j<cols;j++){
-					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
-				}
-			}
-		}
-	}
-	else if(mxIsClass(mxmatrix,"uint8")){
-		/*Dealing with dense matrix: recover pointer and size: */
-		char *pmxmatrix=(char*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-
-		/*Create serial matrix: */
-		if(rows*cols){
-			matrix=xNewZeroInit<double>(rows*cols);
-
-			for(i=0;i<rows;i++){
-				for(j=0;j<cols;j++){
-					matrix[cols*i+j]=(double)pmxmatrix[rows*j+i];
-				}
-			}
-		}
-	}
-	else{
-		_error_("Matlab matrix type Not implemented yet");
-	}
-
-	/*Assign output pointer: */
-	*pmatrix=matrix;
-	*pmatrix_rows=rows;
-	*pmatrix_cols=cols;
-
-	return 1;
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToIssmMat.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToIssmMat.cpp	(revision 15519)
+++ 	(revision )
@@ -1,21 +1,0 @@
-/*!\file MatlabMatrixToIssmDenseMat.cpp
- */
-
-/*Headers:*/
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./matlabio.h"
-
-IssmMat<double>* MatlabMatrixToIssmMat(const mxArray* dataref){
-
-	IssmDenseMat<double>* output=NULL;
-
-	output=new IssmDenseMat<double>();
-	MatlabMatrixToDoubleMatrix(&output->matrix,&output->M,&output->N,dataref);
-	return (IssmMat<double>*)output;
-
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToMatrix.cpp	(revision 15519)
+++ 	(revision )
@@ -1,27 +1,0 @@
-/*!\file MatlabMatrixToMatrix.cpp
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./matlabio.h"
-
-Matrix<double>* MatlabMatrixToMatrix(const mxArray* mxmatrix){
-
-	int dummy;
-	Matrix<double>* matrix=NULL;
-
-	/*allocate matrix object: */
-	matrix=new Matrix<double>();
-
-	#ifdef _HAVE_PETSC_
-	matrix->pmatrix=MatlabMatrixToPetscMat(mxmatrix);
-	#else
-	matrix->imatrix=MatlabMatrixToIssmMat(mxmatrix);
-	#endif
-
-	return matrix;
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToPetscMat.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToPetscMat.cpp	(revision 15519)
+++ 	(revision )
@@ -1,124 +1,0 @@
-/* \file MatlabMatrixToPetscMatrix.cpp
- * \brief: convert a sparse or dense matlab matrix to a serial Petsc matrix:
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-#include "../../c/shared/shared.h"
-
-/*Petsc includes: */
-#ifdef _HAVE_PETSC_
-#include <petscmat.h>
-#include <petscvec.h>
-#include <petscksp.h>
-#endif
-
-/*Matlab includes: */
-#include "./matlabio.h"
-
-PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix){
-
-	int dummy;
-	PetscMat* matrix=new PetscMat();
-
-	MatlabMatrixToPetscMat(&matrix->matrix, &dummy, &dummy, mxmatrix);
-
-	return matrix;
-}
-int MatlabMatrixToPetscMat(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){
-
-	int rows, cols;
-	double *mxmatrix_ptr = NULL;
-	double *tmatrix      = NULL;
-	int ierr;
-	int i,j;
-
-	/*output: */
-	Mat matrix = NULL;
-
-	/*matlab indices: */
-	mwIndex *ir = NULL;
-	mwIndex *jc = NULL;
-	double  *pr = NULL;
-	int     count;
-	int     nnz;
-	int     nz;
-
-	/*petsc indices: */
-	int *idxm = NULL;
-	int *idxn = NULL;
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-		nnz=mxGetNzmax(mxmatrix);
-		if(rows){
-			nz=(int)((double)nnz/(double)rows);
-		}
-		else{
-			nz=0;
-		}
-
-		ierr=MatCreateSeqAIJ(PETSC_COMM_SELF,rows,cols,nz,PETSC_NULL,&matrix);CHKERRQ(ierr);
-
-		/*Now, get ir,jc and pr: */
-		pr=mxGetPr(mxmatrix);
-		ir=mxGetIr(mxmatrix);
-		jc=mxGetJc(mxmatrix);
-
-		/*Now, start inserting data into sparse matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				MatSetValue(matrix,ir[count],i,pr[count],INSERT_VALUES);
-				count++;
-			}
-		}
-	}
-	else{
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-
-		/*transpose, as Petsc now does not allows MAT_COLUMN_ORIENTED matrices in MatSetValues: */
-		tmatrix=xNew<double>(rows*cols);
-		for(i=0;i<cols;i++){
-			for(j=0;j<rows;j++){
-				*(tmatrix+rows*i+j)=*(mxmatrix_ptr+cols*j+i);
-			}
-		}
-
-		/*Create serial matrix: */
-		ierr=MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&matrix);CHKERRQ(ierr);
-
-		/*Insert mxmatrix_ptr values into petsc matrix: */
-		idxm=xNew<int>(rows);
-		idxn=xNew<int>(cols);
-
-		for(i=0;i<rows;i++)idxm[i]=i;
-		for(i=0;i<cols;i++)idxn[i]=i;
-
-		ierr=MatSetValues(matrix,rows,idxm,cols,idxn,tmatrix,INSERT_VALUES); CHKERRQ(ierr);
-
-		xDelete<double>(tmatrix);
-	}
-
-	/*Assemble matrix: */
-	MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY); 
-	MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY);
-
-	/*Assign output pointer: */
-	*pmatrix=matrix;
-	if(pmatrix_rows) *pmatrix_rows=rows;
-	if(pmatrix_cols) *pmatrix_cols=cols;
-
-	return 1;
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabNArrayToNArray.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabNArrayToNArray.cpp	(revision 15519)
+++ 	(revision )
@@ -1,248 +1,0 @@
-/* \file MatlabNArrayToNArray.cpp
- * \brief: convert a sparse or dense matlab n-dimensional array to cpp n-dimensional array
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./matlabio.h"
-
-/*FUNCTION MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{*/
-int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
-
-	int  i,j,rows,cols;
-	int  numel,ndims;
-	int *size,*dims;
-	double* mxmatrix_ptr=NULL;
-	const mwSize* ipt=NULL;
-
-	/*output: */
-	double* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex *ir    = NULL;
-	mwIndex *jc    = NULL;
-	double  *pr    = NULL;
-	int      count;
-
-	/*get Matlab matrix information: */
-	numel=mxGetNumberOfElements(mxmatrix);
-	ndims=mxGetNumberOfDimensions(mxmatrix);
-	ipt  =mxGetDimensions(mxmatrix);
-	size =xNew<int>(ndims);
-	for (i=0;i<ndims;i++) size[i]=(int)ipt[i];
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		rows = mxGetM(mxmatrix);
-		cols = mxGetN(mxmatrix);
-
-		matrix=xNewZeroInit<double>(rows*cols);
-
-		/*Now, get ir,jc and pr: */
-		ir = mxGetIr(mxmatrix);
-		jc = mxGetJc(mxmatrix);
-		pr = mxGetPr(mxmatrix);
-
-		/*Now, start inserting data into double* matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				*(matrix+rows*ir[count]+i)=pr[count];
-				count++;
-			}
-		}
-
-	}
-	else{
-
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
-
-		/*Create serial matrix: */
-		matrix=xNewZeroInit<double>(numel);
-
-		dims=xNew<int>(ndims);
-		for(i=0;i<numel;i++){
-			ColumnWiseDimsFromIndex(dims,i,size,ndims);
-			j = IndexFromRowWiseDims(dims,size,ndims);
-			matrix[j]=(double)mxmatrix_ptr[i];
-		}
-		xDelete<int>(dims);
-	}
-
-	/*Assign output pointer: */
-	*pmatrix       = matrix;
-	*pmatrix_numel = numel;
-	*pmatrix_ndims = ndims;
-	*pmatrix_size  = size;
-
-	return 1;
-}
-/*}}}*/
-/*FUNCTION MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{*/
-int MatlabNArrayToNArray(bool** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
-
-	int  i,j,rows,cols;
-	int  numel,ndims;
-	int *size,*dims;
-	bool* mxmatrix_ptr=NULL;
-	const mwSize* ipt=NULL;
-
-	/*output: */
-	bool* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex *ir    = NULL;
-	mwIndex *jc    = NULL;
-	bool    *pm    = NULL;
-	int      count;
-
-	/*get Matlab matrix information: */
-	numel = mxGetNumberOfElements(mxmatrix);
-	ndims = mxGetNumberOfDimensions(mxmatrix);
-	ipt   = mxGetDimensions(mxmatrix);
-	size  = xNew<int>(ndims);
-	for (i=0;i<ndims;i++) size[i]=(int)ipt[i];
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		rows=mxGetM(mxmatrix);
-		cols=mxGetN(mxmatrix);
-		matrix=xNewZeroInit<bool>(rows*cols);
-
-		/*Now, get ir,jc and pm: */
-		ir=mxGetIr(mxmatrix);
-		jc=mxGetJc(mxmatrix);
-		pm=(bool*)mxGetData(mxmatrix);
-
-		/*Now, start inserting data into bool* matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				matrix[rows*ir[count]+i]=pm[count];
-				count++;
-			}
-		}
-	}
-	else{
-
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=(bool*)mxGetData(mxmatrix);
-
-		/*Create serial matrix: */
-		matrix=xNew<bool>(numel);
-		dims=xNew<int>(ndims);
-		for(i=0;i<numel;i++){
-			ColumnWiseDimsFromIndex(dims,i,size,ndims);
-			j=IndexFromRowWiseDims(dims,size,ndims);
-			matrix[j]=(bool)mxmatrix_ptr[i];
-		}
-		xDelete<int>(dims);
-	}
-
-	/*Assign output pointer: */
-	*pmatrix       = matrix;
-	*pmatrix_numel = numel;
-	*pmatrix_ndims = ndims;
-	*pmatrix_size  = size;
-
-	return 1;
-}
-/*}}}*/
-/*FUNCTION MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){{{*/
-int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix){
-
-	int           i,j,rows,cols;
-	int           numel,ndims;
-	int          *size , *dims;
-	mxChar       *mxmatrix_ptr = NULL;
-	const mwSize *ipt          = NULL;
-
-	/*output: */
-	char* matrix=NULL;
-
-	/*matlab indices: */
-	mwIndex *ir    = NULL;
-	mwIndex *jc    = NULL;
-	char    *pm    = NULL;
-	int      count;
-
-	/*get Matlab matrix information: */
-	numel = mxGetNumberOfElements(mxmatrix);
-	ndims = mxGetNumberOfDimensions(mxmatrix);
-	ipt   = mxGetDimensions(mxmatrix);
-	size  = xNew<int>(ndims);
-	for (i=0;i<ndims;i++) size[i]=(int)ipt[i];
-
-	/*Ok, first check if we are dealing with a sparse or full matrix: */
-	if (mxIsSparse(mxmatrix)){
-
-		/*Dealing with sparse matrix: recover size first: */
-		rows = mxGetM(mxmatrix);
-		cols = mxGetN(mxmatrix);
-		matrix=xNew<char>(rows*cols);
-
-		/*Now, get ir,jc and pm: */
-		ir = mxGetIr(mxmatrix);
-		jc = mxGetJc(mxmatrix);
-		pm = (char*)mxGetData(mxmatrix);
-
-		/*Now, start inserting data into char* matrix: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				matrix[rows*ir[count]+i]=(char)pm[count];
-				count++;
-			}
-		}
-	}
-	else{
-		/*Dealing with dense matrix: recover pointer and size: */
-		mxmatrix_ptr=mxGetChars(mxmatrix);
-
-		/*Create serial matrix: */
-		matrix=xNew<char>(numel+1);
-		matrix[numel]='\0';
-
-		/*looping code adapted from Matlab example explore.c: */
-		int elements_per_page = size[0] * size[1];
-		/* total_number_of_pages = size[2] x size[3] x ... x size[N-1] */
-		int total_number_of_pages = 1;
-		for (i=2; i<ndims; i++) {
-			total_number_of_pages *= size[i];
-		}
-
-		i=0;
-		for (int page=0; page < total_number_of_pages; page++) {
-			int row;
-			/* On each page, walk through each row. */
-			for (row=0; row<size[0]; row++)  {
-				int column;
-				j = (page * elements_per_page) + row;
-
-				/* Walk along each column in the current row. */
-				for (column=0; column<size[1]; column++) {
-					*(matrix+i++)=(char)*(mxmatrix_ptr+j);
-					j += size[0];
-				}
-			}
-		}
-	}
-
-	/*Assign output pointer: */
-	*pmatrix       = matrix;
-	*pmatrix_numel = numel;
-	*pmatrix_ndims = ndims;
-	*pmatrix_size  = size;
-
-	return 1;
-}
-/*}}}*/
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToDoubleVector.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToDoubleVector.cpp	(revision 15519)
+++ 	(revision )
@@ -1,90 +1,0 @@
-/* \file MatlabVectorToDoubleVector.cpp
- * \brief: convert a sparse or dense matlab vector to a serial vector:
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-/*Matlab includes: */
-#include "./matlabio.h"
-
-int MatlabVectorToDoubleVector(double** pvector,int* pvector_rows,const mxArray* mxvector){
-
-	int rows,cols;
-	double* mxvector_ptr=NULL;
-	int ierr;
-	int i,j;
-
-	/*output: */
-	double* vector=NULL;
-
-	/*matlab indices: */
-	mwIndex*    ir=NULL;
-	mwIndex*    jc=NULL;
-	double* pr=NULL;
-	int     count;
-	int     nnz;
-	int     nz;
-
-	/*Ok, first check if we are dealing with a sparse or full vector: */
-	if (mxIsSparse(mxvector)){
-
-		/*Dealing with sparse vector: recover size first: */
-		mxvector_ptr=(double*)mxGetPr(mxvector);
-		rows=mxGetM(mxvector);
-		cols=mxGetN(mxvector);
-		nnz=mxGetNzmax(mxvector);
-
-		/*Check that input is actualy a vector*/
-		if (cols!=1) _error_("input vector of size " << rows << "x" << cols << " should have only one column");
-
-		nz=(int)((double)nnz/(double)rows);
-
-		if(rows){
-			vector=xNewZeroInit<double>(rows);
-
-			/*Now, get ir,jc and pr: */
-			pr=mxGetPr(mxvector);
-			ir=mxGetIr(mxvector);
-			jc=mxGetJc(mxvector);
-
-			/*Now, start inserting data into sparse vector: */
-			count=0;
-			for(i=0;i<cols;i++){
-				for(j=0;j<(jc[i+1]-jc[i]);j++){
-					vector[ir[count]]=pr[count];
-					count++;
-				}
-			}
-		}
-
-	}
-	else{
-
-		/*Dealing with dense vector: recover pointer and size: */
-		mxvector_ptr=(double*)mxGetPr(mxvector);
-		rows=mxGetM(mxvector);
-		cols=mxGetN(mxvector);
-
-		/*Check that input is actualy a vector*/
-		if (cols!=1) _error_("input vector of size " << rows << "x" << cols << " should have only one column");
-
-		/*allocate and memcpy*/
-		if(rows){
-			vector=xNew<double>(rows);
-			memcpy(vector,mxvector_ptr,rows*sizeof(double));
-		}
-		else{
-			vector=NULL;
-		}
-	}
-
-	/*Assign output pointer: */
-	*pvector=vector;
-	*pvector_rows=rows;
-
-	return 1;
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToIssmSeqVec.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToIssmSeqVec.cpp	(revision 15519)
+++ 	(revision )
@@ -1,15 +1,0 @@
-/*!\file MatlabVectorToIssmSeqVec.cpp
- */
-
-/*Headers:*/
-#include "./matlabio.h"
-
-IssmSeqVec<double>* MatlabVectorToIssmSeqVec(const mxArray* dataref){
-
-	IssmSeqVec<double>* output=NULL;
-
-	output=new IssmSeqVec<double>();
-	MatlabVectorToDoubleVector(&output->vector,&output->M,dataref);
-	return output;
-
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToIssmVec.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToIssmVec.cpp	(revision 15519)
+++ 	(revision )
@@ -1,21 +1,0 @@
-/*!\file MatlabVectorToIssmVec.cpp
- */
-
-/*Headers:*/
-#include "./matlabio.h"
-
-IssmVec<double>* MatlabVectorToIssmVec(const mxArray* dataref){
-
-	IssmVec<double>* output=NULL;
-	IssmSeqVec<double>* seqvec=NULL;
-
-	output=new IssmVec<double>();
-
-	seqvec=new IssmSeqVec<double>();
-
-	MatlabVectorToDoubleVector(&seqvec->vector,&seqvec->M,dataref);
-
-	output->vector=(IssmAbsVec<double>*)seqvec;
-	return output;
-
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToPetscVec.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToPetscVec.cpp	(revision 15519)
+++ 	(revision )
@@ -1,100 +1,0 @@
-/* \file MatlabVectorToPetscVector.cpp
- * \brief: convert a sparse or dense matlab vector to a serial Petsc vector:
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-/*Matlab includes: */
-#include "./matlabio.h"
-
-PetscVec* MatlabVectorToPetscVec(const mxArray* mxvector){
-
-	int dummy;
-	PetscVec* vector=new PetscVec();
-
-	MatlabVectorToPetscVec(&vector->vector,&dummy, mxvector);
-
-	return vector;
-}
-
-int MatlabVectorToPetscVec(Vec* pvector,int* pvector_rows,const mxArray* mxvector){
-
-	int rows, cols;
-	double* mxvector_ptr=NULL;
-	int ierr;
-	int i,j;
-
-	/*output: */
-	Vec vector=NULL;
-
-	/*matlab indices: */
-	mwIndex*    ir=NULL;
-	mwIndex*    jc=NULL;
-	double* pr=NULL;
-	int     count;
-	int     nnz;
-	int     nz;
-
-	/*petsc indices: */
-	int* idxm=NULL;
-
-	/*Ok, first check if we are dealing with a sparse or full vector: */
-	if (mxIsSparse(mxvector)){
-
-		/*Dealing with sparse vector: recover size first: */
-		mxvector_ptr=(double*)mxGetPr(mxvector);
-		rows=mxGetM(mxvector);
-		cols=mxGetN(mxvector);
-		nnz=mxGetNzmax(mxvector);
-		nz=(int)((double)nnz/(double)rows);
-
-		ierr=VecCreateSeq(PETSC_COMM_SELF,rows,&vector);CHKERRQ(ierr);
-
-		/*Now, get ir,jc and pr: */
-		pr=mxGetPr(mxvector);
-		ir=mxGetIr(mxvector);
-		jc=mxGetJc(mxvector);
-
-		/*Now, start inserting data into sparse vector: */
-		count=0;
-		for(i=0;i<cols;i++){
-			for(j=0;j<(jc[i+1]-jc[i]);j++){
-				VecSetValue(vector,ir[count],pr[count],INSERT_VALUES);
-				count++;
-			}
-		}
-
-	}
-	else{
-
-		/*Dealing with dense vector: recover pointer and size: */
-		mxvector_ptr=(double*)mxGetPr(mxvector);
-		rows=mxGetM(mxvector);
-		cols=mxGetN(mxvector);
-
-		/*Create serial vector: */
-		ierr=VecCreateSeq(PETSC_COMM_SELF,rows,&vector);CHKERRQ(ierr);
-
-		/*Insert mxvector_ptr values into petsc vector: */
-		idxm=xNew<int>(rows);
-
-		for(i=0;i<rows;i++)idxm[i]=i;
-
-		ierr=VecSetValues(vector,rows,idxm,mxvector_ptr,INSERT_VALUES);CHKERRQ(ierr);
-
-	}
-
-	/*Assemble vector: */
-	VecAssemblyBegin(vector);
-	VecAssemblyEnd(vector);
-
-	/*Assign output pointer: */
-	*pvector=vector;
-	*pvector_rows=rows;
-
-	return 1;
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToVector.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/MatlabVectorToVector.cpp	(revision 15519)
+++ 	(revision )
@@ -1,27 +1,0 @@
-/*!\file MatlabVectorToVector.cpp
- */
-
-/*Headers:*/
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-#include "./matlabio.h"
-
-Vector<double>* MatlabVectorToVector(const mxArray* mxvector){
-
-	int dummy;
-	Vector<double>* vector=NULL;
-
-	/*allocate vector object: */
-	vector=new Vector<double>();
-
-	#ifdef _HAVE_PETSC_
-	vector->pvector=MatlabVectorToPetscVec(mxvector);
-	#else
-	vector->ivector=MatlabVectorToIssmVec(mxvector);
-	#endif
-
-	return vector;
-}
Index: sm/trunk-jpl/src/wrappers/matlab/io/OptionParse.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/OptionParse.cpp	(revision 15519)
+++ 	(revision )
@@ -1,199 +1,0 @@
-/*\file OptionParse.c
- *\brief: functions to parse the mex options.
- */
-#ifdef HAVE_CONFIG_H
-    #include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include <cstring> 
-#include "./matlabio.h"
-
-GenericOption<double>* OptionDoubleParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-	GenericOption<double> *odouble = NULL;
-
-	/*check and parse the name  */
-	odouble=new GenericOption<double>();
-	odouble->name =xNew<char>(strlen(name)+1);
-	memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
-	FetchData(&odouble->value,prhs[0]);
-	odouble->numel=1;
-	odouble->ndims=1;
-	odouble->size=NULL;
-
-	return(odouble);
-}/*}}}*/
-GenericOption<double*>* OptionDoubleArrayParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-	GenericOption<double*> *odouble = NULL;
-
-	/*check and parse the name  */
-	odouble=new GenericOption<double*>();
-	odouble->name =xNew<char>(strlen(name)+1);
-	memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
-
-	/*check and parse the value  */
-	if (!mxIsClass(prhs[0],"double")){
-		_error_("Value of option \"" << odouble->name  << "\" must be class \"double\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-	}
-	FetchData(&odouble->value,&odouble->numel,&odouble->ndims,&odouble->size,prhs[0]);
-
-	return(odouble);
-}/*}}}*/
-GenericOption<bool*>* OptionLogicalParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-	GenericOption<bool*> *ological = NULL;
-
-	/*check and parse the name  */
-	ological=new GenericOption<bool*>();
-	ological->name =xNew<char>(strlen(name)+1);
-	memcpy(ological->name,name,(strlen(name)+1)*sizeof(char));
-
-	/*check and parse the value  */
-	if (!mxIsClass(prhs[0],"logical")){
-		_error_("Value of option \"" << ological->name  << "\" must be class \"logical\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-	}
-	FetchData(&ological->value,&ological->numel,&ological->ndims,&ological->size,prhs[0]);
-
-	return(ological);
-}/*}}}*/
-GenericOption<char*>* OptionCharParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-	GenericOption<char*>  *ochar = NULL;
-
-	/*check and parse the name  */
-	ochar=new GenericOption<char*>();
-	ochar->name =xNew<char>(strlen(name)+1);
-	memcpy(ochar->name,name,(strlen(name)+1)*sizeof(char));
-
-	/*check and parse the value  */
-	if (!mxIsClass(prhs[0],"char")){
-		_error_("Value of option \"" << ochar->name  << "\" must be class \"char\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-	}
-	FetchData(&ochar->value,&ochar->numel,&ochar->ndims,&ochar->size,prhs[0]);
-
-	return(ochar);
-}/*}}}*/
-GenericOption<Options**>* OptionStructParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-	int            i;
-	char           namei[161];
-	Option*                   option      = NULL;
-	GenericOption<Options**>  *ostruct    = NULL;
-	const mwSize  *ipt        = NULL;
-	const mxArray *structi;
-	mwIndex        sindex;
-
-	/*check and parse the name  */
-	ostruct=new GenericOption<Options**>();
-	ostruct->name =xNew<char>(strlen(name)+1);
-	memcpy(ostruct->name,name,(strlen(name)+1)*sizeof(char));
-
-	/*check and parse the value  */
-	if (!mxIsClass(prhs[0],"struct")){
-		_error_("Value of option \"" << ostruct->name  << "\" must be class \"struct\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-	}
-	ostruct->numel=mxGetNumberOfElements(prhs[0]);
-	ostruct->ndims=mxGetNumberOfDimensions(prhs[0]);
-	ipt           =mxGetDimensions(prhs[0]);
-	ostruct->size =xNew<int>(ostruct->ndims);
-	for (i=0; i<ostruct->ndims; i++) ostruct->size[i]=(int)ipt[i];
-	if (ostruct->numel) ostruct->value=xNew<Options*>(ostruct->numel);
-
-	/*loop through and process each element of the struct array  */
-	for (sindex=0; sindex<ostruct->numel; sindex++) {
-		ostruct->value[sindex]=new Options;
-
-		/*loop through and process each field for the element  */
-		for (i=0; i<mxGetNumberOfFields(prhs[0]); i++) {
-			sprintf(namei,"%s.%s",name,mxGetFieldNameByNumber(prhs[0],i));
-			structi=mxGetFieldByNumber(prhs[0],sindex,i);
-
-			option=(Option*)OptionParse(namei,&structi);
-			ostruct->value[sindex]->AddObject((Object*)option);
-			option=NULL;
-		}
-	}
-
-	return(ostruct);
-}/*}}}*/
-GenericOption<Options*>* OptionCellParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-	int            i;
-	int           *dims;
-	char           namei[161];
-	char           cstr[81];
-	GenericOption<Options*> *ocell      = NULL;
-	Option        *option     = NULL;
-	const mwSize  *ipt        = NULL;
-	const mxArray *celli;
-	mwIndex        cindex;
-
-	/*check and parse the name  */
-	ocell=new GenericOption<Options*>();
-	ocell->name =xNew<char>(strlen(name)+1);
-	memcpy(ocell->name,name,(strlen(name)+1)*sizeof(char));
-
-	/*check and parse the value  */
-	if (!mxIsClass(prhs[0],"cell")){
-		_error_("Value of option \"" << ocell->name  << "\" must be class \"cell\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-	}
-
-	ocell->numel=mxGetNumberOfElements(prhs[0]);
-	ocell->ndims=mxGetNumberOfDimensions(prhs[0]);
-	ipt         =mxGetDimensions(prhs[0]);
-	ocell->size =xNew<int>(ocell->ndims);
-	for (i=0; i<ocell->ndims; i++) ocell->size[i]=(int)ipt[i];
-	ocell->value=new Options;
-
-	/*loop through and process each element of the cell array  */
-	dims=xNew<int>(ocell->ndims);
-	for (cindex=0; cindex<ocell->numel; cindex++) {
-		ColumnWiseDimsFromIndex(dims,(int)cindex,ocell->size,ocell->ndims);
-		StringFromDims(cstr,dims,ocell->ndims);
-		#ifdef _INTEL_WIN_
-			_snprintf(namei,161,"%s%s",name,cstr);
-		#else
-			snprintf(namei,161,"%s%s",name,cstr);
-		#endif
-		celli=mxGetCell(prhs[0],cindex);
-
-		option=(Option*)OptionParse(namei,&celli);
-		ocell->value->AddObject((Object*)option);
-		option=NULL;
-	}
-	xDelete<int>(dims);
-
-	return(ocell);
-}/*}}}*/
-Option* OptionParse(char* name, const mxArray* prhs[]){ /*{{{*/
-
-	Option  *option = NULL;
-	mxArray *lhs[1];
-
-	/*parse the value according to the matlab data type  */
-	if     (mxIsClass(prhs[0],"double")  && (mxGetNumberOfElements(prhs[0])==1))
-	 option=(Option*)OptionDoubleParse(name,prhs);
-	else if(mxIsClass(prhs[0],"double")  && (mxGetNumberOfElements(prhs[0])!=1))
-	 option=(Option*)OptionDoubleArrayParse(name,prhs);
-	else if(mxIsClass(prhs[0],"logical"))
-	 option=(Option*)OptionLogicalParse(name,prhs);
-	else if(mxIsClass(prhs[0],"char"))
-	 option=(Option*)OptionCharParse(name,prhs);
-	else if(mxIsClass(prhs[0],"struct"))
-	 option=(Option*)OptionStructParse(name,prhs);
-	else if(mxIsClass(prhs[0],"cell"))
-	 option=(Option*)OptionCellParse(name,prhs);
-	else {
-		_printf0_("  Converting value of option \"" << name << "\" from unrecognized class \"" << mxGetClassName(prhs[0]) << "\" to class \"" << "struct" << "\".\n");
-		if (!mexCallMATLAB(1,lhs,1,(mxArray**)prhs,"struct")) {
-			option=(Option*)OptionStructParse(name,(const mxArray**)lhs);
-			mxDestroyArray(lhs[0]);
-		}
-		else _error_("Second argument value of option \""<< name <<"\" is of unrecognized class \""<< mxGetClassName(prhs[0]) <<"\".");
-	}
-
-	return(option);
-}/*}}}*/
Index: /issm/trunk-jpl/src/wrappers/matlab/io/matlabio.h
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/matlabio.h	(revision 15519)
+++ /issm/trunk-jpl/src/wrappers/matlab/io/matlabio.h	(revision 15520)
@@ -39,10 +39,8 @@
 void FetchData(bool** pmatrix,int* pM,int *pN,const mxArray* dataref);
 void FetchData(bool** pmatrix,int* pnumel,int* pndims,int** psize,const mxArray* dataref);
-void FetchData(Matrix<double>** pmatrix,const mxArray* dataref);
 void FetchData(int** pvector,int* pM,const mxArray* dataref);
 void FetchData(float** pvector,int* pM,const mxArray* dataref);
 void FetchData(double** pvector,int* pM,const mxArray* dataref);
 void FetchData(bool** pvector,int* pM,const mxArray* dataref);
-void FetchData(Vector<double>** pvector,const mxArray* dataref);
 void FetchData(char** pstring,const mxArray* dataref);
 void FetchData(char** pmatrix,int* pnumel,int* pndims,int** psize,const mxArray* dataref);
@@ -73,10 +71,5 @@
 int CheckNumMatlabArguments(int nlhs,int NLHS, int nrhs,int NRHS, const char* THISFUNCTION, void (*function)( void ));
 
-/*Matlab to Matrix routines: */
-Matrix<double>* MatlabMatrixToMatrix(const mxArray* mxmatrix);
-Vector<double>* MatlabVectorToVector(const mxArray* mxvector);
-
 /*Matlab to double* routines: */
-int MatlabVectorToDoubleVector(double** pvector,int* pvector_rows,const mxArray* mxvector);
 int MatlabMatrixToDoubleMatrix(double** pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix);
 int MatlabNArrayToNArray(double** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
@@ -84,16 +77,4 @@
 int MatlabNArrayToNArray(char** pmatrix,int* pmatrix_numel,int* pmatrix_ndims,int** pmatrix_size,const mxArray* mxmatrix);
 
-/*Matlab to IssmDenseMat routines: */
-IssmMat<double>* MatlabMatrixToIssmMat(const mxArray* dataref);
-IssmVec<double>* MatlabVectorToIssmVec(const mxArray* dataref);
-
-/*Matlab to Petsc routines: */
-#ifdef _HAVE_PETSC_
-int MatlabMatrixToPetscMat(Mat* matrix,int* prows,int* pcols, const mxArray* mxmatrix);
-PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix);
-int MatlabVectorToPetscVec(Vec* pvector,int* pvector_rows,const mxArray* mxvector);
-PetscVec* MatlabVectorToPetscVec(const mxArray* mxvector);
-#endif
-
 /*Print*/
 void ApiPrintf(const char* string);
Index: sm/trunk-jpl/src/wrappers/matlab/io/mxGetAssignedField.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/io/mxGetAssignedField.cpp	(revision 15519)
+++ 	(revision )
@@ -1,39 +1,0 @@
-/*!\file: mxGetAssignedField.c: 
- * \brief: abstract interface on parallel side for i/o, so it ressembles the serial i/o.
- *
- * In serial mode, this routine takes care of returning the field coming 
- * from the model. If largesize is 1, we are running out of core models in 
- * matlab, and we need to call the subsref private method from the model object
- * in order to correctly load the data from disk.
- */
-
-#include "./matlabio.h"
-
-mxArray* mxGetAssignedField(const mxArray* pmxa_array,int number,const char* field){
-
-	//output
-	mxArray* mxfield=NULL;
-
-	//input
-	mxArray    *inputs[2];
-	mxArray    *pindex      = NULL;
-	const char *fnames[2];
-	mwSize      ndim        = 2;
-	mwSize      onebyone[2] = {1,1};
-
-	//We want to call the subsasgn method, and get the returned array.This ensures that if we are running 
-	//large sized problems, the data is truly loaded from disk by the model subsasgn class method.
-	inputs[0]=(mxArray*)pmxa_array; //this is the model
-
-	//create index structure used in the assignment (index.type='.' and index.subs='x' for field x for ex)
-	fnames[0] = "type";
-	fnames[1] = "subs";
-	pindex=mxCreateStructArray( ndim,onebyone,2,fnames);
-	mxSetField( pindex, 0, "type",mxCreateString("."));
-	mxSetField( pindex, 0, "subs",mxCreateString(field));
-	inputs[1]=pindex;
-
-	mexCallMATLAB( 1, &mxfield, 2, (mxArray**)inputs, "subsref");
-
-	return mxfield;
-}