Index: /issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp	(revision 11742)
+++ /issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp	(revision 11743)
@@ -26,7 +26,4 @@
 Matrix::Matrix(){
 
-	this->M=0;
-	this->N=0;
-
 	#ifdef _HAVE_PETSC_
 	this->matrix=NULL;
@@ -40,55 +37,42 @@
 /*}}}*/
 /*FUNCTION Matrix::Matrix(int M,int N){{{1*/
-Matrix::Matrix(int pM,int pN){
-
-	this->M=pM;
-	this->N=pN;
-
-	#ifdef _HAVE_PETSC_
-	this->matrix=NewMat(pM,pN);
-	#else
-	this->matrix=new SeqMat(pM,pN);
-	#endif
-	#ifdef _HAVE_ADOLC_
-	this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble));
+Matrix::Matrix(int M,int N){
+
+	#ifdef _HAVE_PETSC_
+	this->matrix=NewMat(M,N);
+	#else
+	this->matrix=new SeqMat(M,N);
+	#endif
+	#ifdef _HAVE_ADOLC_
+	this->amatrix=(adouble*)xmalloc(M*N*sizeof(adouble));
 	#endif
 }
 /*}}}*/
 /*FUNCTION Matrix::Matrix(int M,int N,double sparsity){{{1*/
-Matrix::Matrix(int pM,int pN,double sparsity){
-
-	this->M=pM;
-	this->N=pN;
-
-	#ifdef _HAVE_PETSC_
-	this->matrix=NewMat(pM,pN,sparsity);
-	#else
-	this->matrix=new SeqMat(pM,pN,sparsity);
-	#endif
-	#ifdef _HAVE_ADOLC_
-	this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble));
+Matrix::Matrix(int M,int N,double sparsity){
+
+	#ifdef _HAVE_PETSC_
+	this->matrix=NewMat(M,N,sparsity);
+	#else
+	this->matrix=new SeqMat(M,N,sparsity);
+	#endif
+	#ifdef _HAVE_ADOLC_
+	this->amatrix=(adouble*)xmalloc(M*N*sizeof(adouble));
 	#endif
 }
 /*}}}*/
 /*FUNCTION Matrix::Matrix(double* serial_mat, int M,int N,double sparsity){{{1*/
-Matrix::Matrix(double* serial_mat, int pM,int pN,double sparsity){
-
-	this->M=pM;
-	this->N=pN;
+Matrix::Matrix(double* serial_mat, int M,int N,double sparsity){
 
 	#ifdef _HAVE_PETSC_
 	int     i;
-	int*    idxm=NULL;
-	int*    idxn=NULL;
-
-	this->matrix=NewMat(pM,pN,sparsity);
-
-	idxm=(int*)xmalloc(M*sizeof(int));
-	idxn=(int*)xmalloc(N*sizeof(int));
+
+	int* idxm=(int*)xmalloc(M*sizeof(int));
+	int* idxn=(int*)xmalloc(N*sizeof(int));
 	for(i=0;i<M;i++)idxm[i]=i;
 	for(i=0;i<N;i++)idxn[i]=i;
 
+	this->matrix=NewMat(M,N,sparsity);
 	MatSetValues(this->matrix,M,idxm,N,idxn,serial_mat,INSERT_VALUES);
-
 	MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY);
 	MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY);
@@ -97,24 +81,21 @@
 	xfree((void**)&idxn);
 	#else
-	this->matrix=new SeqMat(serial_mat,pM,pN,sparsity);
-	#endif
-	#ifdef _HAVE_ADOLC_
-	this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble));
+	this->matrix=new SeqMat(serial_mat,M,N,sparsity);
+	#endif
+	#ifdef _HAVE_ADOLC_
+	this->amatrix=(adouble*)xmalloc(M*N*sizeof(adouble));
 	#endif
 }
 /*}}}*/
 /*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){{{1*/
-Matrix::Matrix(int pM,int pN,int connectivity,int numberofdofspernode){
-	
-	this->M=pM;
-	this->N=pN;
-
-	#ifdef _HAVE_PETSC_
-	this->matrix=NewMat(pM,pN,connectivity,numberofdofspernode);
-	#else
-	this->matrix=new SeqMat(pM,pN,connectivity,numberofdofspernode);
-	#endif
-	#ifdef _HAVE_ADOLC_
-	this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble));
+Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){
+	
+	#ifdef _HAVE_PETSC_
+	this->matrix=NewMat(M,N,connectivity,numberofdofspernode);
+	#else
+	this->matrix=new SeqMat(M,N,connectivity,numberofdofspernode);
+	#endif
+	#ifdef _HAVE_ADOLC_
+	this->amatrix=(adouble*)xmalloc(M*N*sizeof(adouble));
 	#endif
 }
@@ -183,5 +164,5 @@
 
 	#ifdef _HAVE_PETSC_
-	MatlabMatrixToPetscMatrix(&matrix->matrix,&matrix->M,&matrix->N,mxmatrix);
+	MatlabMatrixToPetscMatrix(&matrix->matrix,NULL,NULL,mxmatrix);
 	#else
 	matrix->matrix=MatlabMatrixToSeqMat(mxmatrix);
Index: /issm/trunk-jpl/src/c/objects/Numerics/Matrix.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/Matrix.h	(revision 11742)
+++ /issm/trunk-jpl/src/c/objects/Numerics/Matrix.h	(revision 11743)
@@ -1,4 +1,4 @@
 /*!\file:  Matrix.h
- * \brief wrapper to matrix objects. The goal is to control which API (PETSC,Scalpack, Plapack?) 
+ * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
  * implements our underlying matrix format.
  */ 
@@ -32,6 +32,4 @@
 	public:
 	
-		int M,N; 
-
 		#ifdef _HAVE_PETSC_
 		Mat matrix; 
Index: /issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp	(revision 11742)
+++ /issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp	(revision 11743)
@@ -51,20 +51,14 @@
 Vector::Vector(double* serial_vec,int M){
 
-	int i,j;
-
-	#ifdef _HAVE_PETSC_
-		int* idxm=NULL;
+	#ifdef _HAVE_PETSC_
+		int* idxm=(int*)xmalloc(M*sizeof(int));
+		for(int i=0;i<M;i++) idxm[i]=i;
 
 		this->vector=NewVec(M);
-		
-		idxm=(int*)xmalloc(M*sizeof(int));
-		for(i=0;i<M;i++)idxm[i]=i;
 		VecSetValues(this->vector,M,idxm,serial_vec,INSERT_VALUES);
-
 		VecAssemblyBegin(this->vector);
 		VecAssemblyEnd(this->vector);
 
 		xfree((void**)&idxm);
-
 	#else
 		this->vector=new SeqVec(serial_vec,M);
@@ -112,8 +106,6 @@
 
 	#ifdef _HAVE_PETSC_
-	if(!this->vector){
-		printf("Vector size: 0\n");
-	}
-	else VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD);
+	_assert_(this->vector);
+	VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD);
 	#else
 	this->vector->Echo();
@@ -212,5 +204,5 @@
 		
 	#ifdef _HAVE_PETSC_
-		if(!this->vector)*pM=0;
+		_assert_(this->vector);
 		else VecGetSize(this->vector,pM);
 	#else
@@ -218,4 +210,16 @@
 	#endif
 
+}
+/*}}}*/
+/*FUNCTION Vector::IsEmpty{{{1*/
+bool Vector::IsEmpty(void){
+
+	int M;
+	this->GetSize(&M);
+
+	if(M==0) 
+	 return true;
+	else
+	 return false;
 }
 /*}}}*/
@@ -224,6 +228,6 @@
 		
 	#ifdef _HAVE_PETSC_
-		if(!this->vector)*pM=0;
-		else VecGetLocalSize(this->vector,pM);
+		_assert_(this->vector);
+		VecGetLocalSize(this->vector,pM);
 	#else
 		this->vector->GetLocalSize(pM);
@@ -236,12 +240,11 @@
 	
 	Vector* output=NULL;
-	output=new Vector();
-	#ifdef _HAVE_PETSC_
+	#ifdef _HAVE_PETSC_
+		_assert_(this->vector);
 		Vec vec_output=NULL;
-		if(!this->vector){
-			VecDuplicate(this->vector,&vec_output);
-			output->vector=vec_output;
-		}
-	#else
+		VecDuplicate(this->vector,&vec_output);
+		output=new Vector(vec_output);
+	#else
+		output=new Vector();
 		output->vector=this->vector->Duplicate();
 	#endif
Index: /issm/trunk-jpl/src/c/objects/Numerics/Vector.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/Vector.h	(revision 11742)
+++ /issm/trunk-jpl/src/c/objects/Numerics/Vector.h	(revision 11743)
@@ -1,4 +1,4 @@
 /*!\file:  Vector.h
- * \brief wrapper to vector objects. The goal is to control which API (PETSC,Scalpack, Plapack?) 
+ * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
  * implements our underlying vector format.
  */ 
@@ -32,10 +32,10 @@
 	
 		#ifdef _HAVE_PETSC_
-		Vec vector; 
+			Vec vector; 
 		#else
-		SeqVec* vector;
+			SeqVec* vector;
 		#endif
 		#ifdef _HAVE_ADOLC_
-		adouble* avector;
+			adouble* avector;
 		#endif
 
@@ -54,20 +54,21 @@
 		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);
-		Vector* Duplicate(void);
-		void Set(double value);
-		void AXPY(Vector* X, double a);
-		void AYPX(Vector* X, double a);
-		double* ToMPISerial(void);
-		void Copy(Vector* to);
-		double Norm(NormMode norm_type);
-		void Scale(double scale_factor);
-		void PointwiseDivide(Vector* x,Vector* y);
-		double Dot(Vector* vector);
+		void    AXPY(Vector *X, double a);
+		void    AYPX(Vector *X, double a);
+		void    Assemble(void);
+		void    Copy(Vector *to);
+		double  Dot(Vector *vector);
+		Vector *Duplicate(void);
+		void    GetValue(double *pvalue, int dof);
+		void    GetSize(int *pM);
+		void    GetLocalSize(int *pM);
+		bool    IsEmpty(void);
+		double  Norm(NormMode norm_type);
+		void    PointwiseDivide(Vector  *x,Vector*y);
+		void    Scale(double scale_factor);
+		void    Set(double value);
+		void    SetValues(int ssize, int *list, double*values, InsMode mode);
+		void    SetValue(int dof, double value, InsMode mode);
+		double *ToMPISerial(void);
 		/*}}}*/
 };
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp	(revision 11742)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp	(revision 11743)
@@ -116,6 +116,6 @@
 	/*Assign output pointer: */
 	*pmatrix=matrix;
-	*pmatrix_rows=rows;
-	*pmatrix_cols=cols;
+	if(pmatrix_rows) *pmatrix_rows=rows;
+	if(pmatrix_cols) *pmatrix_cols=cols;
 
 	return 1;
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/MatlabVectorToPetscVector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/MatlabVectorToPetscVector.cpp	(revision 11742)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/MatlabVectorToPetscVector.cpp	(revision 11743)
@@ -94,5 +94,4 @@
 	VecAssemblyEnd(vector);
 
-
 	/*Assign output pointer: */
 	*pvector=vector;