diff --git a/CMakeLists.txt b/CMakeLists.txt
index 13ba82a8fc52be155e031afeacaf994dcb1fe540..aff6ab2830bbb0be1552d87cf97b94f6c792f8e0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -58,6 +58,7 @@ opt(MMG3D "Enable MMG3D 3D anisotropic mesh refinement" ${DEFAULT})
 opt(MPEG_ENCODE "Enable built-in MPEG movie encoder" ${DEFAULT})
 opt(MPI "Enable MPI (mostly for parser and solver - mesh generation is sequential)" OFF)
 opt(MSVC_STATIC_RUNTIME "Enable static Visual C++ runtime" OFF)
+opt(MUMPS "Enable MUMPS sparse direct linear solver" OFF)
 opt(NATIVE_FILE_CHOOSER "Enable native file chooser in GUI" ${DEFAULT})
 opt(NETGEN "Enable Netgen 3D frontal mesh generator" ${DEFAULT})
 opt(OCC "Enable Open CASCADE geometrical models" ${DEFAULT})
@@ -124,8 +125,8 @@ set(GMSH_API
   Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
     Solver/crossConfTerm.h Solver/orthogonalTerm.h
     Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h 
-    Solver/linearSystemFull.h Solver/elasticitySolver.h Solver/sparsityPattern.h
-    Solver/groupOfElements.h Solver/linearSystemPETSc.h
+    Solver/linearSystemFull.h Solver/elasticitySolver.h Solver/sparsityPattern.h 
+    Solver/groupOfElements.h Solver/linearSystemPETSc.h Solver/linearSystemMUMPS.h
   Post/PView.h Post/PViewData.h Plugin/PluginManager.h Post/OctreePost.h 
   Post/PViewDataGModel.h Post/PViewOptions.h Post/ColorTable.h
    Numeric/nodalBasis.h
@@ -837,6 +838,20 @@ if(HAVE_SOLVER)
     endif(HAVE_METIS)
   endif(ENABLE_TAUCS)
 
+  if(ENABLE_MUMPS)
+    find_library(DMUMPS_LIB dmumps PATH_SUFFIXES lib)
+    find_library(ZMUMPS_LIB zmumps PATH_SUFFIXES lib)
+    find_path(DMUMPS_INC "dmumps_c.h" PATH_SUFFIXES src include)
+    find_path(ZMUMPS_INC "zmumps_c.h" PATH_SUFFIXES src include)
+    if(DMUMPS_LIB AND DMUMPS_INC AND ZMUMPS_LIB AND ZMUMPS_INC)
+        set_config_option(HAVE_MUMPS "MUMPS")
+        list(APPEND EXTERNAL_LIBRARIES ${DMUMPS_LIB})
+        list(APPEND EXTERNAL_INCLUDES ${DMUMPS_INC})
+        list(APPEND EXTERNAL_LIBRARIES ${ZMUMPS_LIB})
+        list(APPEND EXTERNAL_INCLUDES ${ZMUMPS_INC})
+      endif(DMUMPS_LIB AND DMUMPS_INC AND ZMUMPS_LIB AND ZMUMPS_INC)
+  endif(ENABLE_MUMPS)
+
   if(ENABLE_PETSC)
     if(PETSC_DIR)
       set(ENV_PETSC_DIR ${PETSC_DIR})
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 3ce389e68cff078ff89b82c060fcb37e4819603e..1edfc746828d87a67a18c8c705c2c5a34d2c463e 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -36,6 +36,7 @@
 #cmakedefine HAVE_MMG3D
 #cmakedefine HAVE_MPEG_ENCODE
 #cmakedefine HAVE_MPI
+#cmakedefine HAVE_MUMPS
 #cmakedefine HAVE_NATIVE_FILE_CHOOSER
 #cmakedefine HAVE_NETGEN
 #cmakedefine HAVE_NO_INTPTR_T
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 85990e92fc31c9d47eb1bd9c05e0a8735892ac5c..262a8bb7a8d34da2ad9caf04b756b6186893c05e 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -82,8 +82,7 @@ bool CellComplex::_insertCells(std::vector<MElement*>& elements,
     MElement* element = elements.at(i);
     int dim = element->getDim();
     int type = element->getType();
-    if(//type == TYPE_PYR || type == TYPE_PRI ||
-       type == TYPE_POLYG || type == TYPE_POLYH) {
+    if(type == TYPE_POLYG || type == TYPE_POLYH) {
       Msg::Error("Mesh element type %d not implemented in homology solver", type);
     }
     if(type == TYPE_QUA || type == TYPE_HEX ||
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 420e6e4df5e61da96123e3e298b55e5c89d7318c..03fc23a71b7ba5a6f77181d330795eb3a9df0fca 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -745,6 +745,15 @@ class fullMatrix
     printf("};\n");
   }
 
+  void binarySave (FILE *f) const
+  {
+    fwrite (_data, sizeof(scalar), _r*_c, f);
+  }
+  void binaryLoad (FILE *f)
+  {
+    if(fread (_data, sizeof(scalar), _r*_c, f) != (size_t)_r) return;
+  }
+
   // specific functions for dgshell
   void mult_naiveBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol,
                        const int alpha, const int beta, fullVector<scalar> &c,
diff --git a/Solver/linearSystem.cpp b/Solver/linearSystem.cpp
index b9fa8177661f6a2cd5af6e1fa2ef11483114a4fc..06dfddd34f7e7f0e946207eefd3174e66189c4a1 100644
--- a/Solver/linearSystem.cpp
+++ b/Solver/linearSystem.cpp
@@ -7,9 +7,17 @@
 #include "linearSystemCSR.h"
 #include "linearSystemGMM.h"
 
-void linearSystemBase::setParameter (std::string key, std::string value) 
+void linearSystemBase::setParameter (std::string key, std::string value)
 {
   if (isAllocated())
     Msg::Error("this system is already allocated, parameters cannot be set");
   _parameters[key] = value;
 }
+
+std::string linearSystemBase::getParameter(std::string key) const
+{
+  std::map<std::string, std::string>::const_iterator it;
+  it = this->_parameters.find(key);
+  if(it == this->_parameters.end()) return "";
+  else return it->second;
+}
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index ad9e493b129d1bbe87ddd98c7d0e1f4c84696f1c..8921b7410eb77d9ce1d1b0792d2bf36cfbbfb372 100644
--- a/Solver/linearSystem.h
+++ b/Solver/linearSystem.h
@@ -24,6 +24,7 @@ class linearSystemBase {
   virtual void zeroSolution() = 0;
   virtual int systemSolve() = 0;
   void setParameter (std::string key, std::string value);
+  std::string getParameter(std::string key) const;
   virtual void insertInSparsityPattern(int _row, int _col){};
   virtual double normInfRightHandSide() const = 0;
   virtual double normInfSolution() const {return 0;};
diff --git a/Solver/linearSystemMUMPS.cpp b/Solver/linearSystemMUMPS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bd6b7b11348e607a4ce350c7ef7b2a155bf274f4
--- /dev/null
+++ b/Solver/linearSystemMUMPS.cpp
@@ -0,0 +1,253 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+
+#include <stdio.h>
+#include <math.h>
+#include "linearSystemMUMPS.h"
+
+#if defined(HAVE_MUMPS)
+
+#define USE_COMM_WORLD -987654
+
+void mumpserror(int id, int subid)
+{
+  if(id < 0) {
+    Msg::Error("MUMPS INFO(1) = %d, INFO(2) = %d", id, subid);
+
+    switch (id) {
+    case -6:
+      Msg::Error("Matrix is singular in structure, structural rank: %d",
+                 subid);
+      break;
+     case -10:
+      Msg::Error("Matrix is numerically singular");
+      break;
+    case -13 :
+      Msg::Error("Not enough memory");
+      break;
+    case -40 :
+      Msg::Error("Matrix is not symmetric positive definite");
+      break;
+    default:
+      Msg::Error("Check MUMPS user's guide");
+      break;
+    }
+  }
+}
+
+linearSystemMUMPS<double>::linearSystemMUMPS()
+{
+  _n = 0;
+  _nz = 0;
+}
+
+bool linearSystemMUMPS<double>::isAllocated() const
+{
+  if (_n > 0) return true;
+  else return false;
+}
+
+void linearSystemMUMPS<double>::allocate(int nbRows)
+{
+  _n = nbRows;
+  _b.reserve(_n);
+  _x.reserve(_n);
+  _a.reserve(10*_n);
+  _irn.reserve(10*_n);
+  _jcn.reserve(10*_n);
+  _ij.reserve(_n);
+}
+
+void linearSystemMUMPS<double>::clear()
+{
+  _x.clear();
+  _a.clear();
+  _b.clear();
+  _n = 0;
+  _nz = 0;
+  _irn.clear();
+  _jcn.clear();
+  _ij.clear();
+}
+
+void linearSystemMUMPS<double>::zeroMatrix()
+{
+  _nz = 0;
+  _a.clear();
+  _ij.clear();
+  _irn.clear();
+  _jcn.clear();
+}
+
+void linearSystemMUMPS<double>::zeroRightHandSide()
+{
+  for(unsigned int i=0; i < _b.size(); i++) _b[i] = 0.;
+}
+
+void linearSystemMUMPS<double>::zeroSolution()
+{
+  for(unsigned int i=0; i < _x.size(); i++) _x[i] = 0.;
+}
+
+int linearSystemMUMPS<double>::systemSolve()
+{
+  // MUMPS will overwrite _b with the solution
+  std::vector<DMUMPS_REAL> b = _b;
+
+  DMUMPS_STRUC_C id;
+
+  id.par = 1;
+  const std::string sym = getParameter("symmetry");
+  if(sym == "spd") id.sym = 1;
+  else if(sym == "symmetric") id.sym = 2;
+  else id.sym = 0;
+
+  id.comm_fortran = USE_COMM_WORLD;
+
+  Msg::Info("MUMPS initialization");
+  id.job=-1;
+  dmumps_c(&id);
+  mumpserror(id.info[0], id.info[1]);
+
+  id.n = _n;
+  id.nz = _nz;
+  // Fortran indices start from 1
+  for(unsigned int i=0; i < _irn.size(); i++) _irn[i]++;
+  for(unsigned int i=0; i < _jcn.size(); i++) _jcn[i]++;
+  id.irn= &*_irn.begin();
+  id.jcn= &*_jcn.begin();
+  id.a = &*_a.begin();
+  id.rhs = &*_b.begin();
+
+  // Fortran indices start from 1
+  id.icntl[1-1]=-1;
+  id.icntl[2-1]=-1;
+  id.icntl[3-1]=-1;
+  id.icntl[4-1]=0;
+  id.icntl[5-1]=0;
+  id.icntl[18-1]=0;
+
+  Msg::Info("MUMPS analysis, LU factorization, and back substitution");
+  id.job = 6;
+  dmumps_c(&id);
+  mumpserror(id.info[0], id.info[1]);
+
+  Msg::Info("MUMPS destroy");
+  id.job=-2;
+  dmumps_c(&id);
+  Msg::Info("MUMPS end");
+  mumpserror(id.info[0], id.info[1]);
+
+  _x.clear();
+  _x = _b;
+  _b = b;
+  for(unsigned int i=0; i < _irn.size(); i++) _irn[i]--;
+  for(unsigned int i=0; i < _jcn.size(); i++) _jcn[i]--;
+
+  return 1;
+}
+
+void linearSystemMUMPS<double>::insertInSparsityPattern(int row, int col)
+{
+
+}
+
+double linearSystemMUMPS<double>::normInfRightHandSide() const
+{
+  DMUMPS_REAL norm = 0.;
+  for(unsigned int i=0; i < _b.size(); i++) {
+    DMUMPS_REAL temp = fabs(_b[i]);
+    if(temp > norm) norm = temp;
+  }
+  return norm;
+}
+
+double linearSystemMUMPS<double>::normInfSolution() const
+{
+  DMUMPS_REAL norm = 0.;
+  for(unsigned int i=0; i < _x.size(); i++) {
+    DMUMPS_REAL temp = fabs(_x[i]);
+    if(temp > norm) norm = temp;
+  }
+  return norm;
+}
+
+void linearSystemMUMPS<double>::addToMatrix(int row, int col, const double &val)
+{
+  // MUMPS will sum entries with duplicate (row, col)
+  /*_a.push_back(val);
+  _irn.push_back(row);
+  _jcn.push_back(col);
+  _nz++;*/
+  if((int)_ij.size() <= row) {
+    _a.push_back(val);
+    _irn.push_back(row);
+    _jcn.push_back(col);
+    _ij.resize(row+1);
+    _ij[row][col] = _a.size()-1;
+    _nz++;
+    return;
+  }
+  std::map<int, int>::iterator it = _ij[row].find(col);
+  if(it == _ij[row].end()) {
+    _a.push_back(val);
+    _irn.push_back(row);
+    _jcn.push_back(col);
+    _ij[row][col] = _a.size()-1;
+    _nz++;
+  }
+  else {
+    _a[it->second] += val;
+  }
+}
+
+void linearSystemMUMPS<double>::getFromMatrix(int row, int col, double &val) const
+{
+  //Msg::Error("getFromMatrix not implemented for linearSystemMUMPS");
+  if((int)_ij.size() <= row) {
+    val = 0.;
+    return;
+  }
+  std::map<int, int>::const_iterator it = _ij[row].find(col);
+  if(it == _ij[row].end()) val = 0.;
+  else val = _a[it->second];
+}
+
+void linearSystemMUMPS<double>::addToRightHandSide(int row, const double &val)
+{
+  if((int)_b.size() <= row) {
+    _b.resize(row+1);
+    _b[row] = val;
+  }
+  else {
+    _b[row] += val;
+  }
+}
+
+void linearSystemMUMPS<double>::getFromRightHandSide(int row, double &val) const
+{
+  if((int)_b.size() <= row) val = 0.;
+  else val = _b[row];
+}
+
+void linearSystemMUMPS<double>::getFromSolution(int row, double &val) const
+{
+  if((int)_x.size() <= row) val = 0.;
+  else val = _x[row];
+}
+
+void linearSystemMUMPS<double>::addToSolution(int row, const double &val)
+{
+  if((int)_x.size() <= row) {
+    _x.resize(row+1);
+    _x[row] = val;
+  }
+  else {
+    _x[row] += val;
+  }
+}
+
+#endif
diff --git a/Solver/linearSystemMUMPS.h b/Solver/linearSystemMUMPS.h
new file mode 100644
index 0000000000000000000000000000000000000000..5acc492a00556e3a25d6c6dfce0e2d2a5fd0d51c
--- /dev/null
+++ b/Solver/linearSystemMUMPS.h
@@ -0,0 +1,92 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _LINEAR_SYSTEM_MUMPS_H_
+#define _LINEAR_SYSTEM_MUMPS_H_
+
+// Interface to a linear system with MUMPS
+
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+#include "linearSystem.h"
+
+#if defined(HAVE_MUMPS)
+#include "dmumps_c.h"
+#include "zmumps_c.h"
+
+template <class scalar>
+class linearSystemMUMPS : public linearSystem<scalar> {
+
+ public:
+  linearSystemMUMPS()
+  {
+    Msg::Info("linearSystemMUMPS not implemented for this element type");
+  }
+
+  virtual bool isAllocated() const {return false;}
+  virtual void allocate(int nbRows) {}
+  virtual void clear() {}
+  virtual void zeroMatrix() {}
+  virtual void zeroRightHandSide() {}
+  virtual void zeroSolution() {}
+  virtual int systemSolve() {return 1;}
+  virtual void insertInSparsityPattern(int row, int col) {}
+  virtual double normInfRightHandSide() const {return 0.;}
+  virtual double normInfSolution() const {return 0.;}
+
+  virtual void addToMatrix(int row, int col, const double &val) {}
+  virtual void getFromMatrix(int row, int col, double &val) const {}
+  virtual void addToRightHandSide(int row, const scalar &val) {}
+  virtual void getFromRightHandSide(int row, scalar &val) const {}
+  virtual void getFromSolution(int row, scalar &val) const {}
+  virtual void addToSolution(int row, const scalar &val) {}
+
+};
+
+template <>
+class linearSystemMUMPS<double> : public linearSystem<double> {
+ private:
+
+  int _n;
+  int _nz;
+
+  std::vector<int> _irn;
+  std::vector<int> _jcn;
+
+  std::vector<DMUMPS_REAL> _x;
+  std::vector<DMUMPS_REAL> _b;
+  std::vector<DMUMPS_REAL> _a;
+
+  // _ij[i][j] is the index of _a that is the (i, j) element of
+  // the system matrix
+  std::vector<std::map<int,int> > _ij;
+
+ public:
+  linearSystemMUMPS();
+
+  virtual bool isAllocated() const;
+  virtual void allocate(int nbRows);
+  virtual void clear();
+  virtual void zeroMatrix();
+
+  virtual void zeroRightHandSide();
+  virtual void zeroSolution();
+  virtual int systemSolve();
+  virtual void insertInSparsityPattern(int row, int col);
+  virtual double normInfRightHandSide() const;
+  virtual double normInfSolution() const;
+
+  virtual void addToMatrix(int row, int col, const double &val);
+  virtual void getFromMatrix(int row, int col, double &val) const;
+  virtual void addToRightHandSide(int row, const double &val);
+  virtual void getFromRightHandSide(int row, double &val) const;
+  virtual void getFromSolution(int row, double &val) const;
+  virtual void addToSolution(int row, const double &val);
+
+};
+
+#endif
+
+#endif