Index: ../trunk-jpl/src/c/toolkits/issm/IssmMat.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmMat.h	(revision 23043)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmMat.h	(revision 23044)
@@ -16,7 +16,7 @@
 #include "../../shared/shared.h"
 #include "../ToolkitOptions.h"
 #include "./IssmToolkitUtils.h"
-
+#include "../mumps/mumpsincludes.h"
 /*}}}*/
 
 /*We need to template this class, in case we want to create Matrices that hold
Index: ../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h	(revision 23043)
+++ ../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h	(revision 23044)
@@ -18,8 +18,11 @@
 class Parameters;
 template <class doubletype> class SparseRow;
 
+#ifdef _HAVE_MPI_
 void MpiDenseMumpsSolve(IssmDouble* uf,int uf_M,int uf_n, IssmDouble* Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters);
 void MpiSparseMumpsSolve(IssmDouble* uf,int uf_M,int uf_n, SparseRow<IssmDouble>** Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters);
+#endif
+void SeqDenseMumpsSolve(IssmDouble* uf,int uf_M,int uf_n, IssmDouble* Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters);
 
 #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
 // call back functions:
Index: ../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h	(revision 23043)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h	(revision 23044)
@@ -271,15 +271,23 @@
 			/*Assume we are getting an IssmMpiVec in input, downcast: */
 			pf=(IssmSeqVec<IssmDouble>*)pfin;
 
+			#ifdef _HAVE_MUMPS_
+			/*Assume we have a sequential vec, downcast*/
+			uf=((IssmSeqVec<IssmDouble>*)pfin)->Duplicate();
+			SeqDenseMumpsSolve(uf->vector,uf->M,uf->M, /*stiffness matrix:*/ this->matrix,this->M,this->N,this->M, /*right hand side load vector: */ pf->vector,pf->M,pf->M,parameters);
+			return uf;
+			#endif
+
 			#ifdef _HAVE_GSL_
 			DenseGslSolve(/*output*/ &x,/*stiffness matrix:*/ this->matrix,this->M,this->N, /*right hand side load vector: */ pf->vector,pf->M,parameters);
 
 			uf=new IssmSeqVec<IssmDouble>(x,this->N); xDelete(x);
 			return uf;
-			#else
-				_error_("GSL support not compiled in!");
 			#endif
 
+			_error_("No solver available");
+			return NULL;
+
 		}/*}}}*/
 		#endif
 };
Index: ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp
===================================================================
--- ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp	(revision 23043)
+++ ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp	(revision 23044)
@@ -26,7 +26,9 @@
 void MumpsInit(DMUMPS_STRUC_C &theMumpsStruc){ 
 	theMumpsStruc.par          = 1;  
 	theMumpsStruc.sym          = 0;
+	#ifdef _HAVE_MPI_
 	theMumpsStruc.comm_fortran = MPI_Comm_c2f(IssmComm::GetComm());
+	#endif
 	theMumpsStruc.job          = -1;
 	dmumps_c(&theMumpsStruc);
 }
@@ -34,7 +36,7 @@
 // must be preceded by a call to MumpsInit
 void MumpsSettings(DMUMPS_STRUC_C &theMumpsStruc) { 
 	/*Control statements:{{{ */
-	theMumpsStruc.icntl[1-1] = 6; //error verbose: default 6, 0 or negative -> suppressed
+	theMumpsStruc.icntl[1-1] = 0; //error verbose: default 6, 0 or negative -> suppressed
 	theMumpsStruc.icntl[2-1] = 0; //std verbose: default 1, 0 or negative -> suppressed
 	theMumpsStruc.icntl[3-1] = 0; //global information verbose: default 6, 0 or negative -> suppressed
 	theMumpsStruc.icntl[4-1] = 0; //verbose everything: default is 4
@@ -107,6 +109,7 @@
 		Parameters* parameters);
 #endif 
 
+#ifdef _HAVE_MPI_
 void MpiDenseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ IssmDouble* Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){ /*{{{*/
 
 	/*Variables*/
@@ -200,7 +203,6 @@
 	xDelete<int>(recvcounts);
 	xDelete<int>(displs);
 } /*}}}*/
-
 void MpiSparseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ SparseRow<IssmDouble>** Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){ /*{{{*/
 
 	/*Variables*/
@@ -286,7 +288,72 @@
 	xDelete<int>(recvcounts);
 	xDelete<int>(displs);
 } /*}}}*/
+#endif
+void SeqDenseMumpsSolve(IssmDouble* uf,int uf_M,int uf_m, IssmDouble* Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){/*{{{*/
 
+	/*Variables*/
+	int        *irn_loc = NULL;
+	int        *jcn_loc = NULL;
+	IssmDouble *a_loc   = NULL;
+	IssmDouble *rhs     = NULL;
+
+	/*First, some checks*/
+	if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
+	if (Kff_M!=Kff_m)_error_("stiffness matrix Kff is not serial");
+	if (uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf");
+	if (uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf");
+
+	/*Initialize matrix */
+	/*figure out number of non-zero entries: */
+	int nnz = 0;
+	for(int i=0;i<Kff_M;i++){
+		for(int j=0;j<Kff_N;j++){
+			if(Kff[i*Kff_N+j]!=0)nnz++;
+		}
+	}
+
+	/*Allocate: */
+	if(nnz){
+		irn_loc = xNew<int>(nnz);
+		jcn_loc = xNew<int>(nnz);
+		a_loc   = xNew<IssmDouble>(nnz);
+	}
+
+	/*Populate the triplets: */
+	int count=0;
+	for(int i=0;i<Kff_M;i++){
+		for(int j=0;j<Kff_N;j++){
+			if(Kff[i*Kff_N+j]!=0){
+				irn_loc[count] = i+1; //fortran indexing
+				jcn_loc[count] = j+1; //fortran indexing
+				a_loc[count]   = Kff[i*Kff_N+j];
+				count++;
+			}
+		}
+	}
+
+	/*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
+// AD performance is sensitive to calls to ensurecontiguous.
+// Providing "t" will cause ensurecontiguous to be called.
+#ifdef _HAVE_AD_
+	rhs=xNew<IssmDouble>(pf_M,"t");
+#else
+	rhs=xNew<IssmDouble>(pf_M);
+#endif
+	xMemCpy<IssmDouble>(rhs,pf,Kff_M);
+
+	MumpsSolve(Kff_M,nnz,nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
+
+	/*Now scatter from cpu 0 to all other cpus*/
+	xMemCpy<IssmDouble>(uf,rhs,Kff_M);
+
+	/*Cleanup*/
+	xDelete<int>(irn_loc);
+	xDelete<int>(jcn_loc);
+	xDelete<IssmDouble>(a_loc);
+	xDelete<IssmDouble>(rhs);
+}/*}}}*/
+
 #ifdef _HAVE_ADOLC_
 
 int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y) {