Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 14863)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 14864)
@@ -768,4 +768,9 @@
 
 #}}}
+#Mumps sources  {{{
+mumps_sources=      ./toolkits/mumps\
+					./toolkits/mumps/mumpsincludes.h\
+					./toolkits/mumps/MpiDenseMumpsSolve.cpp
+#}}}
 #Mpi sources  {{{
 mpi_sources= ./toolkits/mpi/mpiincludes.h\
@@ -866,4 +871,8 @@
 endif
 
+if MUMPS
+issm_sources  +=  $(mumps_sources)
+endif
+
 if TRANSIENT
 issm_sources  +=  $(transient_sources)
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h	(revision 14863)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h	(revision 14864)
@@ -208,4 +208,6 @@
 			outvector->vector=this->matrix->Solve(pf->vector,parameters);
 
+			return outvector;
+
 		}/*}}}*/
 
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h	(revision 14863)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h	(revision 14864)
@@ -25,5 +25,6 @@
 #include "../../classes/IssmComm.h"
 #include "../../classes/objects/Bucket.h"
-#include "../../toolkits/toolkits.h"
+#include "../mumps/mumpsincludes.h"
+#include "./IssmMpiVec.h"
 #include <math.h>
 
@@ -69,7 +70,18 @@
 		/*}}}*/
 		/*FUNCTION IssmMpiDenseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
-		IssmMpiDenseMat(int m,int n,int pM,int pN,int* d_nnz,int* o_nnz){
+		IssmMpiDenseMat(int min,int nin,int Min,int Nin,int* d_nnz,int* o_nnz){
 			/*not needed, we are fully dense!: */
-			this->Init(pM,pN);
+			
+			this->buckets=new DataSet();
+
+			this->M=Min;
+			this->N=Nin;
+			this->m=min;
+
+			/*Initialize pointer: */
+			this->matrix=NULL;
+
+			/*Allocate: */
+			if (m*N)this->matrix=xNewZeroInit<doubletype>(this->m*N);
 		}
 		/*}}}*/
@@ -604,13 +616,23 @@
 		IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pfin, Parameters* parameters){
 
+			int i;
+			
+			/*output: */
+			IssmMpiVec<IssmDouble>* uf=NULL;
+			IssmMpiVec<IssmDouble>* pf=NULL;
+
 			/*Assume we are getting an IssmMpiVec in input, downcast: */
-			IssmMpiVec<IssmDouble>* pf=(IssmMpiVec<IssmDouble>*)pfin;
+			pf=(IssmMpiVec<IssmDouble>*)pfin;
+
+			/*Initialize output: */
+			uf=pf->Duplicate();
 
 			/*Let's try and use the MUMPS solver here: */
 			#ifdef _HAVE_MUMPS_
-
+			MpiDenseMumpsSolve(/*output*/ 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);
 			#else
 			_error_("IssmMpiDenseMat solver requires MUMPS solver");
 			#endif
+			return (IssmAbsVec<IssmDouble>*)uf;
 
 		}/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h	(revision 14863)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h	(revision 14864)
@@ -60,4 +60,9 @@
 		}
 		/*}}}*/
+		/*FUNCTION IssmMpiVec(int m,int M){{{*/
+		IssmMpiVec(int min, int Min){
+			this->Init(min,true);
+		}
+		/*}}}*/
 		/*FUNCTION IssmMpiVec(int M,bool fromlocalsize){{{*/
 		IssmMpiVec(int Min, bool fromlocalsize){
@@ -71,4 +76,18 @@
 
 			if(this->M){
+				this->vector=xNew<doubletype>(this->m);
+				xMemCpy<doubletype>(this->vector,buffer,this->m);
+			}
+		}
+		/*}}}*/
+		/*FUNCTION IssmMpiVec(doubletype* serial_vec,int M,int m){{{*/
+		IssmMpiVec(doubletype* buffer,int Min,int min){
+
+			this->vector=NULL;
+			this->buckets=new DataSet();
+			this->M=Min;
+			this->m=min;
+
+			if(this->m){
 				this->vector=xNew<doubletype>(this->m);
 				xMemCpy<doubletype>(this->vector,buffer,this->m);
@@ -119,4 +138,5 @@
 						printf("row %i %g",j,this->vector[j]);
 					}
+					printf("\n");
 				}
 				MPI_Barrier(IssmComm::GetComm());
@@ -437,5 +457,5 @@
 		IssmMpiVec<doubletype>* Duplicate(void){
 
-			return new IssmMpiVec<doubletype>(this->vector,this->M);
+			return new IssmMpiVec<doubletype>(this->vector,this->M,this->m);
 
 		}
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h	(revision 14863)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h	(revision 14864)
@@ -42,20 +42,5 @@
 		/*IssmVec constructors, destructors*/
 		IssmVec(){ /*{{{*/
-
-			switch(IssmVecTypeFromToolkitOptions()){
-
-				case SeqEnum: 
-					this->vector=new IssmSeqVec<doubletype>();
-					break;
-				case MpiEnum:
-					#ifdef _HAVE_MPI_
-					this->vector=new IssmMpiVec<doubletype>();
-					#else
-					_error_("Mpi vector requires compilation of MPI!");
-					#endif
-					break;
-				default:
-					_error_("vector type not supported yet!");
-			}
+			this->vector=NULL;
 		}
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp	(revision 14864)
+++ /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp	(revision 14864)
@@ -0,0 +1,158 @@
+/*!\file MpiDenseMumpsSolve.cpp
+ * \brief: solve dense matrix system with MUMPS
+ */
+
+/*Header files: {{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+#include "../../include/types.h"
+#include "../../include/macros.h"
+#include "../../shared/Alloc/alloc.h"
+#include "../../shared/Exceptions/exceptions.h"
+#include "../../classes/IssmComm.h"
+#include "../mpi/patches/mpipatches.h"
+
+/*Mumps header files: */
+#include "dmumps_c.h"
+
+/*}}}*/
+
+#define JOB_INIT -1 
+#define JOB_SOLVE 6
+#define JOB_END  -2
+
+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){
+
+	/*Variables: */
+
+
+	MPI_Comm   comm;
+	int        my_rank;
+	int        num_procs;
+	int        i;
+	int        j;
+	int         nnz       ,local_nnz;
+	int        *irn_loc = NULL;
+	int        *jcn_loc = NULL;
+	IssmDouble *a_loc   = NULL;
+	int         count;
+	int         lower_row;
+	int         upper_row;
+	IssmDouble* rhs=NULL;
+	int*        recvcounts=NULL;
+	int*        displs=NULL;
+
+	my_rank=IssmComm::GetRank();
+	num_procs=IssmComm::GetSize();
+	comm=IssmComm::GetComm();
+	
+	/*First, some checks:{{{ */
+	if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
+	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 mumps: {{{*/
+	DMUMPS_STRUC_C id;
+	id.job          = JOB_INIT;
+	id.par          = 1;  
+	id.sym          = 1;
+	id.comm_fortran = MPI_Comm_c2f(comm);
+	dmumps_c(&id);
+	/*}}}*/
+	/*Control statements:{{{ */
+	id.icntl[1-1] = 6; //error verbose
+	id.icntl[2-1] = 1; //std verbose
+	id.icntl[4-1] = 4; //verbose everything
+	id.icntl[5-1] = 0;
+	id.icntl[18-1] = 3;
+	
+	id.icntl[20-1] = 0;
+	id.icntl[21-1] = 0;
+	id.icntl[30-1] = 0;
+	/*}}}*/
+	/*Initialize matrix:{{{ */
+	id.n=Kff_M;
+
+	/*figure out number of non-zero entries: */
+	local_nnz=0;
+	for(i=0;i<Kff_m;i++){
+		for(j=0;j<Kff_N;j++){
+			if (Kff[i*Kff_N+j]!=0)local_nnz++;
+		}
+	}
+
+	MPI_Reduce(&local_nnz,&nnz,1,MPI_INT,MPI_SUM,0,comm);
+	MPI_Bcast(&nnz,1,MPI_INT,0,comm);
+	id.nz=nnz;
+	id.nz_loc=local_nnz;
+
+	/*Allocate: */
+	irn_loc=xNew<int>(local_nnz);
+	jcn_loc=xNew<int>(local_nnz);
+	a_loc=xNew<IssmDouble>(local_nnz);
+
+	/*Populate the triplets: */
+	GetOwnershipBoundariesFromRange(&lower_row,&upper_row,Kff_m,comm);
+	count=0;
+	for(i=0;i<Kff_m;i++){
+		for(j=0;j<Kff_N;j++){
+			if (Kff[i*Kff_N+j]!=0){
+				irn_loc[count]=lower_row+i+1; //fortran indexing
+				jcn_loc[count]=j+1; //fortran indexing
+				a_loc[count]=Kff[i*Kff_N+j];
+				count++;
+			}
+		}
+	}
+	id.irn_loc=irn_loc;
+	id.jcn_loc=jcn_loc;
+	id.a_loc=a_loc;
+
+
+	/*Deal with right hand side. We need to MPI_Gather it onto cpu 0: */
+	rhs=xNew<IssmDouble>(pf_M);
+
+	recvcounts=xNew<int>(num_procs);
+	displs=xNew<int>(num_procs);
+
+	/*recvcounts:*/
+	MPI_Allgather(&pf_m,1,MPI_INT,recvcounts,1,MPI_INT,comm);
+	
+	/*displs: */
+	MPI_Allgather(&lower_row,1,MPI_INT,displs,1,MPI_INT,comm);
+
+	/*Gather:*/
+	MPI_Gatherv(pf, pf_m, MPI_DOUBLE, rhs, recvcounts, displs, MPI_DOUBLE,0,comm);
+	id.rhs=rhs;
+	id.nrhs=1;
+	id.lrhs=1;
+
+	/*}}}*/
+	/*Solve system: {{{*/
+	id.job = JOB_SOLVE;
+	dmumps_c (&id);
+	/*}}}*/
+	/*Now scatter from cpu 0 to all other cpus: {{{*/
+	MPI_Scatterv( rhs, recvcounts, displs, MPI_DOUBLE, uf, uf_m, MPI_DOUBLE, 0, comm); 
+
+	/*}}}*/
+	/*}}}*/
+	/*Cleanup: {{{*/
+	id.job = JOB_END; 
+	dmumps_c (&id);
+
+	xDelete<int>(irn_loc);
+	xDelete<int>(jcn_loc);
+	xDelete<IssmDouble>(a_loc);
+	xDelete<IssmDouble>(rhs);
+	xDelete<int>(recvcounts);
+	xDelete<int>(displs);
+
+	/*}}}*/
+}
+
Index: /issm/trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h	(revision 14864)
+++ /issm/trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h	(revision 14864)
@@ -0,0 +1,10 @@
+/* \file mumpsincludes.h
+ * \brief all includes for MUMPS layer
+ */
+
+#ifndef _MUMPS_INCLUDES_H_
+#define _MUMPS_INCLUDES_H_
+
+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);
+
+#endif