diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 9d4f5e845db5f1a8841c8d4bb455eeaeb29bbf95..de4f0f8d66a09e76d38e2f40e813b5ebb92a9562 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -115,9 +115,6 @@ const char *Get_GmshShortLicense(){ return GMSH_SHORT_LICENSE; }
 std::string Get_GmshBuildOptions()
 {
   std::string opt;
-#if defined(HAVE_GSL)
-  opt += "GSL ";
-#endif
 #if defined(HAVE_NETGEN)
   opt += "NETGEN ";
 #endif
diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index 0741d6f2ea3ecd3afcf73f673798c45a2bead67d..9d3ab7bdf04e3948fc4775b91707cfa866ee32bf 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -7,6 +7,7 @@
 #include <time.h>
 #include "GmshConfig.h"
 #include "GmshDefines.h"
+#include "GmshPredicates.h"
 #include "GModel.h"
 #include "GmshMessage.h"
 #include "OpenFile.h"
@@ -14,7 +15,6 @@
 #include "Options.h"
 #include "CommandLine.h"
 #include "OS.h"
-#include "Numeric.h"
 #include "Generator.h"
 #include "Field.h"
 #include "Context.h"
@@ -48,8 +48,8 @@ int GmshInitialize(int argc, char **argv)
   PluginManager::instance()->registerDefaultPlugins();
 #endif
 
-  // Initialize numeric library (gsl, robust predicates)
-  Init_Numeric();
+  // Initialize robust predicates
+  gmsh::exactinit();
 
   if(dummy) delete dummy;
   return 1;
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index dd637bc3108e4fea5801a4268da01a5d5ca30ec7..1a1e3ca2d3ea65a96ddc4b645d2dd0861162b363 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -9,12 +9,10 @@
 #undef HAVE_64BIT_SIZE_T
 #undef HAVE_ANN
 #undef HAVE_BLAS
-#undef HAVE_CBLAS
 #undef HAVE_CHACO
 #undef HAVE_FLTK
 #undef HAVE_FOURIER_MODEL
 #undef HAVE_GMM
-#undef HAVE_GSL
 #undef HAVE_LAPACK
 #undef HAVE_LIBCGNS
 #undef HAVE_LIBJPEG
diff --git a/Common/Makefile b/Common/Makefile
index a1654365c3bd4ae1f6a5dd9cc61ca757c2f821f6..6478ed92a819b5073a327a75190bb8a3b85d3efd 100644
--- a/Common/Makefile
+++ b/Common/Makefile
@@ -57,21 +57,19 @@ depend:
 	rm -f Makefile.new
 
 # DO NOT DELETE THIS LINE
-Gmsh${OBJEXT}: Gmsh.cpp GmshConfig.h GmshDefines.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
-  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
-  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+Gmsh${OBJEXT}: Gmsh.cpp GmshConfig.h GmshDefines.h ../Numeric/GmshPredicates.h \
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
+  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h \
+  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h GmshMessage.h OpenFile.h CreateFile.h Options.h \
-  ../Post/ColorTable.h CommandLine.h OS.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Mesh/Generator.h \
-  ../Mesh/Field.h ../Post/PView.h Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h ../Mesh/meshPartition.h GmshDaemon.h \
-  ../Plugin/PluginManager.h
+  ../Post/ColorTable.h CommandLine.h OS.h ../Mesh/Generator.h \
+  ../Mesh/Field.h ../Common/GmshConfig.h ../Post/PView.h Context.h \
+  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
+  ../Mesh/meshPartition.h GmshDaemon.h ../Plugin/PluginManager.h
 GmshMessage${OBJEXT}: GmshMessage.cpp GmshConfig.h GmshMessage.h Gmsh.h \
   Options.h ../Post/ColorTable.h Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h ../Common/GmshConfig.h OS.h \
@@ -92,9 +90,9 @@ Options${OBJEXT}: Options.cpp GmshConfig.h GmshDefines.h ../Geo/GModel.h \
   ../Mesh/meshPartitionOptions.h ../Common/GmshConfig.h Options.h \
   ../Post/ColorTable.h DefaultOptions.h ../Mesh/Field.h ../Post/PView.h \
   ../Mesh/BackgroundMesh.h ../Post/PViewOptions.h ../Post/ColorTable.h \
-  ../Post/PViewData.h ../Numeric/GmshMatrix.h ../Common/GmshMessage.h \
-  ../Post/adaptiveData.h ../Plugin/PluginManager.h ../Plugin/Plugin.h \
-  ../Common/Options.h ../Post/PViewDataList.h ../Post/PViewData.h \
+  ../Post/PViewData.h ../Numeric/GmshMatrix.h ../Post/adaptiveData.h \
+  ../Plugin/PluginManager.h ../Plugin/Plugin.h ../Common/Options.h \
+  ../Common/GmshMessage.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Common/ListUtils.h ../Fltk/GUI.h ../Fltk/Solvers.h \
   ../Fltk/menuWindow.h ../Fltk/popupButton.h ../Fltk/graphicWindow.h \
   ../Fltk/openglWindow.h ../Graphics/drawContext.h ../Fltk/optionWindow.h \
@@ -118,8 +116,7 @@ OpenFile${OBJEXT}: OpenFile.cpp GmshConfig.h GmshMessage.h ../Geo/Geo.h \
   ../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
   ../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
   ../Common/ListUtils.h ../Geo/SPoint2.h ../Geo/ExtrudeParams.h \
   ../Common/SmoothData.h ../Geo/GModel.h ../Geo/GVertex.h \
@@ -155,11 +152,9 @@ CreateFile${OBJEXT}: CreateFile.cpp GmshConfig.h GmshMessage.h ../Geo/GModel.h \
 VertexArray${OBJEXT}: VertexArray.cpp VertexArray.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h ../Common/GmshConfig.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
 SmoothData${OBJEXT}: SmoothData.cpp SmoothData.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h
 Octree${OBJEXT}: Octree.cpp Octree.h OctreeInternals.h
 OctreeInternals${OBJEXT}: OctreeInternals.cpp GmshMessage.h OctreeInternals.h
 StringUtils${OBJEXT}: StringUtils.cpp StringUtils.h GmshMessage.h
diff --git a/Fltk/Makefile b/Fltk/Makefile
index 4dba34befbefbdf3341477ae5d1285816263215e..95f4b0a66dc708d9e8691eab17bde945b6158970 100644
--- a/Fltk/Makefile
+++ b/Fltk/Makefile
@@ -111,9 +111,9 @@ graphicWindow${OBJEXT}: graphicWindow.cpp GUI.h graphicWindow.h openglWindow.h \
   ../Graphics/drawContext.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
   paletteWindow.h mainWindow.h menuWindow.h popupButton.h messageWindow.h \
   manipWindow.h extraDialogs.h Draw.h ../Post/PView.h ../Post/PViewData.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
-  ../Common/OS.h ../Common/Options.h ../Post/ColorTable.h \
-  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/OS.h \
+  ../Common/Options.h ../Post/ColorTable.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 openglWindow${OBJEXT}: openglWindow.cpp openglWindow.h ../Graphics/drawContext.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h graphicWindow.h manipWindow.h \
   contextWindow.h visibilityWindow.h ../Common/GmshConfig.h \
@@ -128,9 +128,8 @@ openglWindow${OBJEXT}: openglWindow.cpp openglWindow.h ../Graphics/drawContext.h
   ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
   ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
-  ../Numeric/Gauss.h Draw.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h GUI.h \
-  ../Common/VertexArray.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Numeric/Gauss.h Draw.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  GUI.h ../Common/VertexArray.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h
 menuWindow${OBJEXT}: menuWindow.cpp ../Common/GmshConfig.h \
   ../Common/GmshMessage.h GUI.h Draw.h menuWindow.h popupButton.h \
@@ -342,12 +341,12 @@ classificationEditor${OBJEXT}: classificationEditor.cpp GUI.h \
   ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
   ../Post/ColorTable.h paletteWindow.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h Draw.h \
-  ../Common/Options.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h ../Mesh/meshGFaceDelaunayInsertion.h \
-  ../Mesh/meshGFaceOptimize.h ../Mesh/meshGFaceDelaunayInsertion.h \
-  ../Geo/discreteEdge.h ../Geo/GModel.h ../Geo/GEdge.h \
-  ../Geo/discreteFace.h ../Geo/GModel.h ../Geo/GFace.h
+  ../Numeric/GmshMatrix.h Draw.h ../Common/Options.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
+  ../Mesh/meshGFaceDelaunayInsertion.h ../Mesh/meshGFaceOptimize.h \
+  ../Mesh/meshGFaceDelaunayInsertion.h ../Geo/discreteEdge.h \
+  ../Geo/GModel.h ../Geo/GEdge.h ../Geo/discreteFace.h ../Geo/GModel.h \
+  ../Geo/GFace.h
 partitionDialog${OBJEXT}: partitionDialog.cpp ../Common/GmshConfig.h \
   ../Common/GmshDefines.h ../Common/GmshMessage.h GUI.h paletteWindow.h \
   ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 4a44b0c977b613a030d136b2dc2b883f36e433f6..08549a695a5dd96b9ef71d770c08990206f6970d 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -254,8 +254,8 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
   ym /= (double)ndata;
   zm /= (double)ndata;
 
-  Double_Matrix U(ndata, na), V(na, na);
-  Double_Vector sigma(na);
+  gmshMatrix<double> U(ndata, na), V(na, na);
+  gmshVector<double> sigma(na);
   for(int i = 0; i < ndata; i++) {
     U(i, 0) = points[i].x() - xm;
     U(i, 1) = points[i].y() - ym;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 89b4e9a0591d3d64858995d8926d5ed603936b79..31e0fc7c92442e00eb6d024cb19d674f207ad987 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2882,7 +2882,7 @@ struct PointSurface{
   Surface *s;
 };
 
-static void projectPS(Double_Vector &x, Double_Vector &res, void *data)
+static void projectPS(gmshVector<double> &x, gmshVector<double> &res, void *data)
 {
   PointSurface *ps = (PointSurface*)data;
   Vertex c = InterpolateSurface(ps->s, x(0), x(1), 0, 0);
@@ -2900,7 +2900,7 @@ static void projectPS(Double_Vector &x, Double_Vector &res, void *data)
 
 bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2])
 {
-  Double_Vector x(2);
+  gmshVector<double> x(2);
   x(0) = uv[0];
   x(1) = uv[1];
   PointSurface ps = {&p, s};
@@ -3031,7 +3031,7 @@ struct CurveSurface{
   Surface *s;
 };
 
-static void intersectCS(Double_Vector &uvt, Double_Vector &res, void *data)
+static void intersectCS(gmshVector<double> &uvt, gmshVector<double> &res, void *data)
 {
   CurveSurface *cs = (CurveSurface*)data;
   Vertex vs = InterpolateSurface(cs->s, uvt(0), uvt(1), 0, 0);
@@ -3054,7 +3054,7 @@ bool IntersectCurvesWithSurface(List_T *curve_ids, int surface_id, List_T *shape
     Curve *c = FindCurve((int)curve_id);
     if(c){
       CurveSurface cs = {c, s};
-      Double_Vector uvt(3);
+      gmshVector<double> uvt(3);
       uvt(0) = 0.5;
       uvt(1) = 0.5;
       uvt(2) = 0.5;
diff --git a/Geo/Makefile b/Geo/Makefile
index 93ea9eec425b7b7b5b5efdb2d6abb5d4c948928f..2bbabfe354224a2e94366101f3bdf3245dd93fb7 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -91,27 +91,25 @@ GFace${OBJEXT}: GFace.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   GRegion.h MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
   ../Common/VertexArray.h ../Geo/SVector3.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Numeric/GaussLegendre1D.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h
+  ../Numeric/GmshMatrix.h ../Numeric/GaussLegendre1D.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 GFaceCompound${OBJEXT}: GFaceCompound.cpp ../Common/GmshConfig.h GFaceCompound.h \
   GFace.h GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h \
   GEdgeLoop.h GEdge.h GVertex.h SPoint2.h SVector3.h Pair.h \
   ../Numeric/gmshAssembler.h ../Numeric/gmshLinearSystem.h \
   ../Numeric/gmshLaplace.h ../Numeric/gmshTermOfFormulation.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h \
-  ../Numeric/gmshFunction.h ../Common/Gmsh.h ../Common/GmshMessage.h \
-  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+  ../Numeric/GmshMatrix.h ../Numeric/gmshFunction.h ../Common/Gmsh.h \
+  ../Common/GmshMessage.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h \
+  ../Geo/GFace.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
   ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
   ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
   ../Numeric/Gauss.h ../Numeric/GmshMatrix.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h ../Common/Octree.h \
-  ../Common/OctreeInternals.h ../Numeric/gmshLinearSystemGmm.h \
-  ../Numeric/gmshLinearSystem.h ../Numeric/gmshLinearSystemFull.h \
-  ../Numeric/gmshLinearSystem.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h ../Common/Octree.h ../Common/OctreeInternals.h \
+  ../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
+  ../Numeric/gmshLinearSystemFull.h ../Numeric/gmshLinearSystem.h \
+  ../Numeric/GmshMatrix.h
 GRegion${OBJEXT}: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
@@ -121,42 +119,40 @@ GRegion${OBJEXT}: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
 gmshVertex${OBJEXT}: gmshVertex.cpp GFace.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h SPoint2.h \
   SVector3.h Pair.h gmshVertex.h Geo.h ../Common/GmshDefines.h \
-  gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
-  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
-  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
-  GeoInterpolation.h MVertex.h MElement.h MEdge.h MFace.h \
-  ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h
+  gmshSurface.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshConfig.h ../Common/ListUtils.h ../Common/TreeUtils.h \
+  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
+  ../Common/SmoothData.h GeoInterpolation.h ../Common/GmshMessage.h \
+  MVertex.h MElement.h MEdge.h MFace.h ../Numeric/FunctionSpace.h \
+  ../Numeric/GmshMatrix.h ../Numeric/Gauss.h
 gmshEdge${OBJEXT}: gmshEdge.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h GFaceCompound.h gmshEdge.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/ListUtils.h \
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
   ExtrudeParams.h ../Common/SmoothData.h GeoInterpolation.h \
-  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+  ../Common/GmshMessage.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/meshPartitionOptions.h
 gmshFace${OBJEXT}: gmshFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h gmshFace.h Geo.h ../Common/GmshDefines.h \
-  gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
-  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
-  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
-  GeoInterpolation.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h
+  gmshSurface.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshConfig.h ../Common/ListUtils.h ../Common/TreeUtils.h \
+  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
+  ../Common/SmoothData.h GeoInterpolation.h ../Common/GmshMessage.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 gmshRegion${OBJEXT}: gmshRegion.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h gmshRegion.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/ListUtils.h \
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
-  ExtrudeParams.h ../Common/SmoothData.h
+  ExtrudeParams.h ../Common/SmoothData.h ../Common/GmshMessage.h
 gmshSurface${OBJEXT}: gmshSurface.cpp ../Common/GmshConfig.h \
   ../Common/GmshMessage.h gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 OCCVertex${OBJEXT}: OCCVertex.cpp ../Common/GmshConfig.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h MVertex.h MElement.h \
@@ -172,9 +168,9 @@ OCCFace${OBJEXT}: OCCFace.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   GModel.h GVertex.h GEntity.h Range.h SPoint3.h SBoundingBox3d.h \
   GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h GEdgeLoop.h Pair.h \
   GRegion.h OCCVertex.h OCCIncludes.h OCCEdge.h OCCFace.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/VertexArray.h ../Geo/SVector3.h \
-  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/VertexArray.h \
+  ../Geo/SVector3.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/meshPartitionOptions.h
 OCCRegion${OBJEXT}: OCCRegion.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   GModel.h GVertex.h GEntity.h Range.h SPoint3.h SBoundingBox3d.h \
   GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h GEdgeLoop.h Pair.h \
@@ -184,25 +180,24 @@ discreteEdge${OBJEXT}: discreteEdge.cpp ../Common/GmshConfig.h \
   Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
-  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h
+  ../Numeric/GmshMatrix.h ../Common/ListUtils.h ../Common/TreeUtils.h \
+  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
+  ../Common/SmoothData.h
 discreteFace${OBJEXT}: discreteFace.cpp ../Common/GmshConfig.h \
   ../Common/GmshMessage.h discreteFace.h GModel.h GVertex.h GEntity.h \
   Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
-  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h
+  ../Numeric/GmshMatrix.h ../Common/ListUtils.h ../Common/TreeUtils.h \
+  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
+  ../Common/SmoothData.h
 discreteRegion${OBJEXT}: discreteRegion.cpp ../Common/GmshConfig.h \
   discreteRegion.h GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h Geo.h ../Common/GmshDefines.h \
-  gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Common/ListUtils.h \
-  ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
-  ExtrudeParams.h ../Common/SmoothData.h
+  gmshSurface.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
+  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h
 fourierEdge${OBJEXT}: fourierEdge.cpp ../Common/GmshConfig.h fourierEdge.h \
   GEdge.h GEntity.h Range.h SPoint3.h SBoundingBox3d.h GVertex.h GPoint.h \
   SPoint2.h SVector3.h GModel.h GFace.h GEdgeLoop.h Pair.h GRegion.h \
@@ -223,21 +218,20 @@ GModel${OBJEXT}: GModel.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   GRegion.h MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
   discreteRegion.h discreteFace.h discreteEdge.h discreteVertex.h \
-  gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/Octree.h ../Common/OctreeInternals.h \
-  ../Common/SmoothData.h ../Mesh/Field.h ../Post/PView.h ../Geo/SPoint3.h \
-  ../Mesh/Generator.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h
+  gmshSurface.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/Octree.h ../Common/OctreeInternals.h ../Common/SmoothData.h \
+  ../Mesh/Field.h ../Post/PView.h ../Geo/SPoint3.h ../Mesh/Generator.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 GModelIO_Geo${OBJEXT}: GModelIO_Geo.cpp ../Common/GmshConfig.h \
   ../Common/GmshMessage.h GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h Geo.h ../Common/GmshDefines.h \
-  gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/ListUtils.h ../Common/TreeUtils.h \
-  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
-  ../Common/SmoothData.h ../Common/OpenFile.h gmshVertex.h gmshFace.h \
-  GFaceCompound.h gmshEdge.h gmshRegion.h ../Mesh/Field.h ../Post/PView.h \
-  ../Geo/SPoint3.h ../Parser/Parser.h
+  gmshSurface.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
+  ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
+  ../Common/OpenFile.h gmshVertex.h gmshFace.h GFaceCompound.h gmshEdge.h \
+  gmshRegion.h ../Mesh/Field.h ../Post/PView.h ../Geo/SPoint3.h \
+  ../Parser/Parser.h
 GModelIO_Mesh${OBJEXT}: GModelIO_Mesh.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/GmshDefines.h MElement.h \
@@ -270,14 +264,12 @@ GModelIO_MED${OBJEXT}: GModelIO_MED.cpp ../Common/GmshConfig.h \
 ExtrudeParams${OBJEXT}: ExtrudeParams.cpp ../Common/GmshMessage.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/ListUtils.h ../Common/TreeUtils.h \
-  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
-  ../Common/SmoothData.h
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
+  ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
+  ExtrudeParams.h ../Common/SmoothData.h
 Geo${OBJEXT}: Geo.cpp ../Common/GmshMessage.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/MallocUtils.h Geo.h \
-  ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/MallocUtils.h \
+  Geo.h ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
   ExtrudeParams.h ../Common/SmoothData.h GModel.h GVertex.h GEntity.h \
@@ -285,8 +277,7 @@ Geo${OBJEXT}: Geo.cpp ../Common/GmshMessage.h ../Numeric/Numeric.h \
   ../Mesh/Field.h ../Post/PView.h ../Geo/SPoint3.h ../Common/Context.h \
   ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 GeoStringInterface${OBJEXT}: GeoStringInterface.cpp ../Common/GmshConfig.h \
-  ../Common/GmshMessage.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshMessage.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
   ../Common/StringUtils.h Geo.h ../Common/GmshDefines.h gmshSurface.h \
   Pair.h Range.h SPoint2.h SPoint3.h SVector3.h SBoundingBox3d.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
@@ -298,10 +289,10 @@ GeoStringInterface${OBJEXT}: GeoStringInterface.cpp ../Common/GmshConfig.h \
 GeoInterpolation${OBJEXT}: GeoInterpolation.cpp ../Common/GmshMessage.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/ListUtils.h ../Common/TreeUtils.h \
-  ../Common/avl.h ../Common/ListUtils.h ExtrudeParams.h \
-  ../Common/SmoothData.h GeoInterpolation.h GeoStringInterface.h
+  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
+  ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
+  ExtrudeParams.h ../Common/SmoothData.h GeoInterpolation.h \
+  GeoStringInterface.h
 findLinks${OBJEXT}: findLinks.cpp ../Common/GmshMessage.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/TreeUtils.h \
@@ -311,16 +302,14 @@ MVertex${OBJEXT}: MVertex.cpp MVertex.h SPoint2.h SPoint3.h GVertex.h GEntity.h
   GEdgeLoop.h Pair.h GFaceCompound.h ../Common/GmshMessage.h \
   ../Common/StringUtils.h
 MFace${OBJEXT}: MFace.cpp ../Common/GmshConfig.h MFace.h MVertex.h SPoint2.h \
-  SPoint3.h SVector3.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+  SPoint3.h SVector3.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 MElement${OBJEXT}: MElement.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   MElement.h ../Common/GmshDefines.h MVertex.h SPoint2.h SPoint3.h \
   MEdge.h SVector3.h MFace.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Numeric/Gauss.h GEntity.h Range.h \
   SBoundingBox3d.h GFace.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h Pair.h \
-  ../Common/StringUtils.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h \
+  ../Common/StringUtils.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
   ../Numeric/GaussLegendre1D.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h ../Mesh/qualityMeasures.h \
   ../Mesh/meshGFaceDelaunayInsertion.h ../Geo/MElement.h \
diff --git a/Graphics/Makefile b/Graphics/Makefile
index fb34d92fcb9614d328d7af6451f40e931e198369..fe7c1f0299287ff58cdc18e9fee96c515a36f9bb 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -59,8 +59,8 @@ depend:
 
 # DO NOT DELETE THIS LINE
 Trackball${OBJEXT}: Trackball.cpp Trackball.h
-Iso${OBJEXT}: Iso.cpp ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h
+Iso${OBJEXT}: Iso.cpp ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshConfig.h
 ReadImg${OBJEXT}: ReadImg.cpp ReadImg.h ../Common/GmshMessage.h ../Post/PView.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
@@ -68,16 +68,16 @@ ReadImg${OBJEXT}: ReadImg.cpp ReadImg.h ../Common/GmshMessage.h ../Post/PView.h
 drawContext${OBJEXT}: drawContext.cpp ../Common/GmshMessage.h drawContext.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h Trackball.h \
   ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
-  ../Common/GmshConfig.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
-  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
-  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
-  ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Post/PView.h \
-  ../Post/PViewOptions.h ../Post/ColorTable.h
+  ../Common/GmshConfig.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/GPoint.h \
+  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
+  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Post/PView.h ../Post/PViewOptions.h \
+  ../Post/ColorTable.h
 drawMesh${OBJEXT}: drawMesh.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Common/GmshMessage.h ../Common/GmshDefines.h \
   ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
@@ -108,12 +108,11 @@ drawGeom${OBJEXT}: drawGeom.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h
 drawPost${OBJEXT}: drawPost.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   drawContext.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h Iso.h ../Post/PView.h ../Post/PViewOptions.h \
-  ../Post/ColorTable.h ../Post/PViewData.h ../Common/VertexArray.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Common/SmoothData.h \
-  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
-  gl2ps.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h Iso.h ../Post/PView.h \
+  ../Post/PViewOptions.h ../Post/ColorTable.h ../Post/PViewData.h \
+  ../Common/VertexArray.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+  ../Common/SmoothData.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/meshPartitionOptions.h gl2ps.h
 drawAxes${OBJEXT}: drawAxes.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h \
   ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/GPoint.h \
@@ -124,24 +123,21 @@ drawAxes${OBJEXT}: drawAxes.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h ../Common/GmshConfig.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h gl2ps.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h gl2ps.h
 drawScales${OBJEXT}: drawScales.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Post/PView.h ../Post/PViewOptions.h \
   ../Post/ColorTable.h ../Post/PViewData.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h gl2ps.h
+  ../Common/GmshConfig.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/meshPartitionOptions.h gl2ps.h
 drawGraph2d${OBJEXT}: drawGraph2d.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Post/PView.h ../Post/PViewOptions.h \
   ../Post/ColorTable.h ../Post/PViewData.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h gl2ps.h \
-  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
+  ../Common/GmshConfig.h gl2ps.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/meshPartitionOptions.h ../Numeric/Numeric.h \
   ../Numeric/GmshMatrix.h
 drawGlyph${OBJEXT}: drawGlyph.cpp drawContext.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Fltk/Draw.h ../Common/GmshDefines.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
   ../Common/StringUtils.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h gl2ps.h
 gl2ps${OBJEXT}: gl2ps.cpp gl2ps.h ../Common/GmshConfig.h
diff --git a/Makefile b/Makefile
index 08d94dcdf4b71d1454850937265b964ce31bc532..84050841bff781f0f72822b4758b8cd010dc2006 100644
--- a/Makefile
+++ b/Makefile
@@ -38,7 +38,7 @@ GMSH_EMBEDDED = ${GMSH_API} Geo/discrete*.cpp\
                 Geo/GEdgeLoop.cpp Geo/GFace.cpp Geo/GRegion.cpp\
                 Geo/MElement.cpp Geo/MFace.cpp Geo/MVertex.cpp\
                 Common/StringUtils.{cpp,h}\
-                Numeric/NumericEmbedded.{cpp,h} Numeric/FunctionSpace.cpp\
+                Numeric/Numeric.{cpp,h} Numeric/FunctionSpace.cpp\
                 utils/embed/GmshEmbedded.{cpp,h} utils/embed/Makefile
 
 # Main building rules
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 02353ae25bb0da6c4eee45b6f7cb62551ffa490a..8e8b718777b6aa92d9965a5e67f64d4992f2cb78 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -5,9 +5,10 @@
 
 #include <math.h>
 #include <stdio.h>
+#include "GmshMessage.h"
+#include "GmshPredicates.h"
 #include "Numeric.h"
 #include "BDS.h"
-#include "GmshMessage.h"
 #include "GFace.h"
 #include "meshGFaceDelaunayInsertion.h"
 #include "qualityMeasures.h"
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index 98046d129a0b5dfc950b7cc8a5afb5a531edd4dc..95c47d1930e7bc7ecc9f23e55813068383d2d81d 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -16,6 +16,7 @@
 // value to avoid 3 aligned points or 4 cocyclical points!
 
 #include "GmshMessage.h"
+#include "GmshPredicates.h"
 #include "Numeric.h"
 #include "DivideAndConquer.h"
 #include "MallocUtils.h"
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index e73cbe252bb1bafb8a17bad6a860cbaf6a302fc0..ce58393c77739d8d3d005890c4f330a7ad61f360 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -40,7 +40,7 @@ static bool mappingIsInvertible(MTetrahedron *e)
     if (det0 * detN <= 0.) return false;
   }
 
-  const Double_Matrix &points = e->getFunctionSpace()->points;
+  const gmshMatrix<double> &points = e->getFunctionSpace()->points;
 
   for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
     const double u = points(i,0);
@@ -62,7 +62,7 @@ static double mylength(GEdge *ge, int i, double *u)
   return ge->length(u[i], u[i+1], 10);
 }
 
-static void myresid(int N, GEdge *ge, double *u, Double_Vector &r)
+static void myresid(int N, GEdge *ge, double *u, gmshVector<double> &r)
 {
   double L[100];
   for (int i = 0; i < N - 1; i++) L[i] = mylength(ge, i, u);
@@ -89,10 +89,10 @@ static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
 
   // create the tangent matrix
   const int M = N - 2;
-  Double_Matrix J(M, M);
-  Double_Vector DU(M);
-  Double_Vector R(M);
-  Double_Vector Rp(M);
+  gmshMatrix<double> J(M, M);
+  gmshVector<double> DU(M);
+  gmshVector<double> R(M);
+  gmshVector<double> Rp(M);
   
   int iter = 1;
 
@@ -134,7 +134,7 @@ static double mylength(GFace *gf, int i, double *u, double *v)
   return gf->length(SPoint2(u[i], v[i]), SPoint2(u[i + 1], v[i + 1]), 10);
 }
 
-static void myresid(int N, GFace *gf, double *u, double *v, Double_Vector &r)
+static void myresid(int N, GFace *gf, double *u, double *v, gmshVector<double> &r)
 {
   double L[100];
   for (int i = 0; i < N - 1; i++) L[i] = mylength(gf, i, u, v);  
@@ -168,10 +168,10 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
 
   // create the tangent matrix
   const int M = N - 2;
-  Double_Matrix J(M, M);
-  Double_Vector DU(M);
-  Double_Vector R(M);
-  Double_Vector Rp(M);
+  gmshMatrix<double> J(M, M);
+  gmshVector<double> DU(M);
+  gmshVector<double> R(M);
+  gmshVector<double> Rp(M);
   
   int iter = 1;
 
@@ -401,7 +401,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
      gf->geomType() == GEntity::BoundaryLayerSurface)
     linear = true;
 
-  Double_Matrix points;
+  gmshMatrix<double> points;
   int start = 0;
 
   switch (nPts){
@@ -550,7 +550,7 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
                             faceContainer &faceVertices, edgeContainer &edgeVertices,
                             bool linear, int nPts = 1)
 {
-  Double_Matrix points;
+  gmshMatrix<double> points;
   int start = 0;
   
   switch (nPts){
@@ -643,7 +643,7 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
 static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele, 
                               std::vector<MVertex*> &vr, bool linear, int nPts = 1)
 {
-  Double_Matrix points;
+  gmshMatrix<double> points;
   int start = 0;
 
   switch (nPts){
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 27a17868cb576b16086bff97140a2be093e22846..06475aaf827d9b9890d048b7dfdf552f4c758b31 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -76,16 +76,15 @@ depend:
 
 # DO NOT DELETE THIS LINE
 Generator${OBJEXT}: Generator.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h ../Common/OS.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
-  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
-  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h ../Common/OS.h \
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
+  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h \
+  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
   ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
   ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
@@ -100,8 +99,7 @@ Field${OBJEXT}: Field.cpp ../Common/GmshConfig.h \
   ../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
   ../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Common/ListUtils.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
   ../Geo/SPoint2.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
   ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
@@ -111,9 +109,9 @@ Field${OBJEXT}: Field.cpp ../Common/GmshConfig.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
   ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h \
   ../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
-  ../Post/OctreePost.h ../Common/Octree.h ../Common/OctreeInternals.h \
-  ../Post/PViewDataList.h ../Post/PViewData.h ../Geo/MVertex.h \
-  ../Geo/SPoint2.h ../Geo/SPoint3.h
+  ../Common/GmshMessage.h ../Post/OctreePost.h ../Common/Octree.h \
+  ../Common/OctreeInternals.h ../Post/PViewDataList.h ../Post/PViewData.h \
+  ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h
 gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp HighOrder.h \
   ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
   ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
@@ -137,7 +135,7 @@ gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp HighOrder.h \
   ../Numeric/GmshMatrix.h ../Numeric/gmshLinearSystemGmm.h \
   ../Numeric/gmshLinearSystem.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/meshPartitionOptions.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 meshGEdge${OBJEXT}: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -152,9 +150,8 @@ meshGEdge${OBJEXT}: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
   ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
-  BackgroundMesh.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h
+  BackgroundMesh.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 meshGEdgeExtruded${OBJEXT}: meshGEdgeExtruded.cpp ../Geo/GModel.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -186,9 +183,8 @@ meshGFace${OBJEXT}: meshGFace.cpp meshGFace.h meshGFaceBDS.h \
   ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h \
   ../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
   ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h BDS.h qualityMeasures.h Field.h ../Post/PView.h \
-  ../Common/OS.h HighOrder.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h BDS.h qualityMeasures.h \
+  Field.h ../Post/PView.h ../Common/OS.h HighOrder.h
 meshGFaceTransfinite${OBJEXT}: meshGFaceTransfinite.cpp meshGFace.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -202,8 +198,7 @@ meshGFaceTransfinite${OBJEXT}: meshGFaceTransfinite.cpp meshGFace.h \
   ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
   ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
 meshGFaceExtruded${OBJEXT}: meshGFaceExtruded.cpp ../Geo/GModel.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -219,31 +214,30 @@ meshGFaceExtruded${OBJEXT}: meshGFaceExtruded.cpp ../Geo/GModel.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
   ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Common/Context.h \
   ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
-meshGFaceBDS${OBJEXT}: meshGFaceBDS.cpp meshGFace.h meshGFaceOptimize.h \
+meshGFaceBDS${OBJEXT}: meshGFaceBDS.cpp ../Common/GmshMessage.h \
+  ../Numeric/GmshPredicates.h meshGFace.h meshGFaceOptimize.h \
   ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
   ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
-  meshGFaceDelaunayInsertion.h BackgroundMesh.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
-  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h \
-  ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
-  ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h \
-  ../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h BDS.h qualityMeasures.h Field.h ../Post/PView.h \
-  ../Common/OS.h
+  ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshConfig.h ../Numeric/Gauss.h meshGFaceDelaunayInsertion.h \
+  BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
+  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
+  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  BDS.h qualityMeasures.h Field.h ../Post/PView.h ../Common/OS.h
 meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp \
-  BackgroundMesh.h meshGFaceDelaunayInsertion.h ../Geo/MElement.h \
-  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
-  ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+  ../Common/GmshMessage.h ../Numeric/GmshPredicates.h BackgroundMesh.h \
+  meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
+  ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
   meshGFaceOptimize.h meshGFace.h ../Geo/GFace.h ../Geo/GEntity.h \
   ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
@@ -251,7 +245,7 @@ meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp \
   ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/GPoint.h \
   ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
   ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
   ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
   ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
@@ -267,9 +261,8 @@ meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
   ../Geo/Pair.h BackgroundMesh.h Generator.h
 meshGFaceQuadrilateralize${OBJEXT}: meshGFaceQuadrilateralize.cpp \
   meshGFaceQuadrilateralize.h ../Common/GmshMessage.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Geo/GFace.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
+  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
   ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
@@ -286,31 +279,30 @@ meshGRegion${OBJEXT}: meshGRegion.cpp ../Common/GmshConfig.h \
   ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
-  ../Numeric/Gauss.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h BackgroundMesh.h qualityMeasures.h \
-  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
-  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
-  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/gmshRegion.h ../Geo/Geo.h \
-  ../Geo/gmshSurface.h ../Geo/Pair.h ../Geo/Range.h ../Geo/SPoint2.h \
-  ../Geo/SPoint3.h ../Geo/SVector3.h ../Geo/SBoundingBox3d.h \
-  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
-  ../Common/ListUtils.h ../Geo/SPoint2.h ../Geo/ExtrudeParams.h \
-  ../Common/SmoothData.h ../Geo/GRegion.h BDS.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+  ../Numeric/Gauss.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  BackgroundMesh.h qualityMeasures.h ../Geo/GModel.h ../Geo/GVertex.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
+  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h \
+  ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
+  ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h \
+  ../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
+  ../Geo/gmshRegion.h ../Geo/Geo.h ../Geo/gmshSurface.h ../Geo/Pair.h \
+  ../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
+  ../Geo/SBoundingBox3d.h ../Common/ListUtils.h ../Common/TreeUtils.h \
+  ../Common/avl.h ../Common/ListUtils.h ../Geo/SPoint2.h \
+  ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Geo/GRegion.h BDS.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 meshGRegionDelaunayInsertion${OBJEXT}: meshGRegionDelaunayInsertion.cpp \
-  ../Common/OS.h BackgroundMesh.h meshGRegion.h meshGRegionLocalMeshMod.h \
+  ../Common/GmshMessage.h ../Numeric/GmshPredicates.h ../Common/OS.h \
+  BackgroundMesh.h meshGRegion.h meshGRegionLocalMeshMod.h \
   meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
   ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
   ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
+  ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshConfig.h ../Numeric/Gauss.h ../Numeric/Numeric.h \
   ../Numeric/GmshMatrix.h qualityMeasures.h ../Geo/GModel.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -368,28 +360,26 @@ meshGRegionLocalMeshMod${OBJEXT}: meshGRegionLocalMeshMod.cpp \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h BackgroundMesh.h qualityMeasures.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h BackgroundMesh.h \
+  qualityMeasures.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GRegion.h \
   ../Geo/GEntity.h
 DivideAndConquer${OBJEXT}: DivideAndConquer.cpp ../Common/GmshMessage.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
+  ../Numeric/GmshPredicates.h ../Numeric/Numeric.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h DivideAndConquer.h \
   ../Common/MallocUtils.h
 BackgroundMesh${OBJEXT}: BackgroundMesh.cpp ../Common/GmshMessage.h \
-  BackgroundMesh.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
-  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
-  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h Field.h ../Post/PView.h
+  BackgroundMesh.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshConfig.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  ../Mesh/meshPartitionOptions.h ../Geo/GVertex.h ../Geo/GEntity.h \
+  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
+  ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
+  ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
+  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GModel.h ../Geo/GVertex.h \
+  ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h ../Geo/GEntity.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h Field.h ../Post/PView.h
 qualityMeasures${OBJEXT}: qualityMeasures.cpp qualityMeasures.h BDS.h \
   ../Common/GmshMessage.h ../Geo/MVertex.h ../Geo/SPoint2.h \
   ../Geo/SPoint3.h ../Geo/MElement.h ../Common/GmshDefines.h \
@@ -397,7 +387,7 @@ qualityMeasures${OBJEXT}: qualityMeasures.cpp qualityMeasures.h BDS.h \
   ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Numeric/Gauss.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 BoundaryLayers${OBJEXT}: BoundaryLayers.cpp ../Geo/GModel.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -413,8 +403,8 @@ BoundaryLayers${OBJEXT}: BoundaryLayers.cpp ../Geo/GModel.h ../Geo/GVertex.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
   BoundaryLayers.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
   meshGEdge.h meshGFace.h
-BDS${OBJEXT}: BDS.cpp ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
+BDS${OBJEXT}: BDS.cpp ../Common/GmshMessage.h ../Numeric/GmshPredicates.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
   BDS.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
   ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
@@ -440,9 +430,8 @@ HighOrder${OBJEXT}: HighOrder.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
   ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Geo/MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Numeric/Gauss.h \
-  ../Common/OS.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h
+  ../Common/OS.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 meshPartition${OBJEXT}: meshPartition.cpp ../Common/GmshConfig.h ../Geo/GModel.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index c839a20c8ee9ef91be537b1963352040dadd1aad..a26ff2ba29e2c583831e3a019baf3407ea0b76e2 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -377,10 +377,11 @@ void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
 
   double F = -smooth_obj_HighOrderN(uv, &vdN);
   if (F < eps){
-    double val;
-    minimize_N(2 * vdN.v.size(), smooth_obj_HighOrderN, 
-               deriv_smoothing_objective_function_HighOrderN, 
-               &vdN, 1, uv, val);
+    double val = 0.;
+    Msg::Error("Fletcher-Reeves minimizer routine must be reimplemented");
+    //minimize_N(2 * vdN.v.size(), smooth_obj_HighOrderN, 
+    //          deriv_smoothing_objective_function_HighOrderN, 
+    //          &vdN, 1, uv, val);
     double Fafter = -smooth_obj_HighOrderN(uv, &vdN);
     printf("%12.5E %12.5E\n", F, Fafter);
     if (F < Fafter){
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index ef3a02be2555f208587730f63cf4a098d2a42f53..272b7a634976ae7ca23ad80fc4eca765cf569de9 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -4,6 +4,8 @@
 // bugs and problems to <gmsh@geuz.org>.
 
 #include <stdlib.h>
+#include "GmshMessage.h"
+#include "GmshPredicates.h"
 #include "meshGFace.h"
 #include "meshGFaceOptimize.h"
 #include "BackgroundMesh.h"
@@ -15,7 +17,6 @@
 #include "Context.h"
 #include "GPoint.h"
 #include "GModel.h"
-#include "GmshMessage.h"
 #include "Numeric.h"
 #include "BDS.h"
 #include "qualityMeasures.h"
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index ae57c85d7a9ee2eacda33198053bc3d0a6bf33c1..f96d7d54a0d3b4e4a07ba14693f341780c7f5e4a 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -6,13 +6,14 @@
 #include <set>
 #include <map>
 #include <algorithm>
+#include "GmshMessage.h"
+#include "GmshPredicates.h"
 #include "BackgroundMesh.h"
 #include "meshGFaceDelaunayInsertion.h"
 #include "meshGFaceOptimize.h"
 #include "meshGFace.h"
 #include "GFace.h"
 #include "Numeric.h"
-#include "GmshMessage.h"
 
 const double LIMIT_ = 0.5 * sqrt(2.0);
 
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 1bdfa916a37dac1e0443a18c9f167adae5ca305b..dc453f0cf0b6983d43a83366ea39fe0297a096db 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -3,6 +3,11 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <set>
+#include <map>
+#include <algorithm>
+#include "GmshMessage.h"
+#include "GmshPredicates.h"
 #include "OS.h"
 #include "BackgroundMesh.h"
 #include "meshGRegion.h"
@@ -11,10 +16,6 @@
 #include "GModel.h"
 #include "GRegion.h"
 #include "Numeric.h"
-#include "GmshMessage.h"
-#include <set>
-#include <map>
-#include <algorithm>
 
 int MTet4::inCircumSphere(const double *p) const
 {
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 29f87ea1b345cdac4298c630dc8f390c8593af00..ba34b8ae9db70a69092ad9c3b7bec0021e5cd8c7 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -791,9 +791,10 @@ bool gmshSmoothVertexOptimize(MTet4 *t,
 
   double xyzopti[3] = {vd.v->x(), vd.v->y(), vd.v->z()};
 
-  double val;
-  minimize_N(3, smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4,
-             xyzopti, val);
+  double val = 0.;
+  Msg::Error("Fletcher-Reeves minimizer routine must be reimplemented");
+  //minimize_N(3, smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4,
+  //         xyzopti, val);
 
   double vTot = 0;
 
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 055f8e0798f1e0553412b9d998869476513858f5..96c9640c7834da267308de3936e387a33499db7b 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -223,7 +223,7 @@ double qmDistorsionOfMapping (MTriangle *e)
     const double di  = mesh_functional_distorsion (e,u,v);
     dmin = (i==0)? di : std::min(dmin,di);
   }
-  const Double_Matrix& points = e->getFunctionSpace()->points;
+  const gmshMatrix<double>& points = e->getFunctionSpace()->points;
 
   for (int i=0;i<e->getNumPrimaryVertices();i++) {
     const double u = points(i,0);
@@ -271,7 +271,7 @@ double qmDistorsionOfMapping (MTetrahedron *e)
     dmin = (i==0)? di : std::min(dmin,di);
   }
   
-  const Double_Matrix& points = e->getFunctionSpace()->points;
+  const gmshMatrix<double>& points = e->getFunctionSpace()->points;
 
   for (int i=0;i<e->getNumPrimaryVertices();i++) {
     const double u = points(i,0);
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index 6cab3d5030333820b27fecbbd6d54422aca42bb4..eb18b8e322155b3630a2f00238836f2d5d291cff 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -11,16 +11,16 @@
 #include "GmshDefines.h"
 #include "GmshMessage.h"
 
-Double_Matrix generate1DMonomials(int order)
+gmshMatrix<double> generate1DMonomials(int order)
 {
-  Double_Matrix monomials(order + 1, 1);
+  gmshMatrix<double> monomials(order + 1, 1);
   for (int i = 0; i < order + 1; i++) monomials(i, 0) = i;
   return monomials;
 }
 
-Double_Matrix generate1DPoints(int order)
+gmshMatrix<double> generate1DPoints(int order)
 {
-  Double_Matrix line(order + 1, 1);
+  gmshMatrix<double> line(order + 1, 1);
   line(0, 0) = -1.;
   line(1, 0) =  1.;
   double dd = 2. / order;
@@ -28,9 +28,9 @@ Double_Matrix generate1DPoints(int order)
   return line;
 }
 
-Double_Matrix generatePascalTriangle(int order)
+gmshMatrix<double> generatePascalTriangle(int order)
 {
-  Double_Matrix monomials((order + 1) * (order + 2) / 2, 2);
+  gmshMatrix<double> monomials((order + 1) * (order + 2) / 2, 2);
   int index = 0;
   for (int i = 0; i <= order; i++) {
     for (int j = 0; j <= i; j++) {
@@ -44,9 +44,9 @@ Double_Matrix generatePascalTriangle(int order)
 
 // generate the exterior hull of the Pascal triangle for Serendipity element
 
-Double_Matrix generatePascalSerendipityTriangle(int order)
+gmshMatrix<double> generatePascalSerendipityTriangle(int order)
 {
-  Double_Matrix monomials(3 * order, 2);
+  gmshMatrix<double> monomials(3 * order, 2);
 
   monomials(0, 0) = 0;
   monomials(0, 1) = 0;
@@ -74,24 +74,24 @@ Double_Matrix generatePascalSerendipityTriangle(int order)
 
 // generate the monomials subspace of all monomials of order exactly == p
 
-Double_Matrix generateMonomialSubspace(int dim, int p)
+gmshMatrix<double> generateMonomialSubspace(int dim, int p)
 {
-  Double_Matrix monomials;
+  gmshMatrix<double> monomials;
   
   switch (dim) {
   case 1:
-    monomials = Double_Matrix(1, 1);
+    monomials = gmshMatrix<double>(1, 1);
     monomials(0, 0) = p;
     break;
   case 2:
-    monomials = Double_Matrix(p + 1, 2);
+    monomials = gmshMatrix<double>(p + 1, 2);
     for (int k = 0; k <= p; k++) {
       monomials(k, 0) = p - k;
       monomials(k, 1) = k;
     }
     break;
   case 3:
-    monomials = Double_Matrix((p + 1) * (p + 2) / 2, 3);
+    monomials = gmshMatrix<double>((p + 1) * (p + 2) / 2, 3);
     int index = 0;
     for (int i = 0; i <= p; i++) {
       for (int k = 0; k <= p - i; k++) {
@@ -108,11 +108,11 @@ Double_Matrix generateMonomialSubspace(int dim, int p)
 
 // generate external hull of the Pascal tetrahedron
 
-Double_Matrix generatePascalSerendipityTetrahedron(int order)
+gmshMatrix<double> generatePascalSerendipityTetrahedron(int order)
 {
   int nbMonomials = 4 + 6 * std::max(0, order - 1) + 
     4 * std::max(0, (order - 2) * (order - 1) / 2);
-  Double_Matrix monomials(nbMonomials, 3);
+  gmshMatrix<double> monomials(nbMonomials, 3);
 
   monomials.set_all(0);
   int index = 1;
@@ -128,7 +128,7 @@ Double_Matrix generatePascalSerendipityTetrahedron(int order)
       }
     }
   }
-  Double_Matrix monomialsMaxOrder = generateMonomialSubspace(3, order);
+  gmshMatrix<double> monomialsMaxOrder = generateMonomialSubspace(3, order);
   int nbMaxOrder = monomialsMaxOrder.size1();
   monomials.copy(monomialsMaxOrder, 0, nbMaxOrder, 0, 3, index, 0);
   return monomials;
@@ -136,15 +136,15 @@ Double_Matrix generatePascalSerendipityTetrahedron(int order)
 
 // generate Pascal tetrahedron
 
-Double_Matrix generatePascalTetrahedron(int order)
+gmshMatrix<double> generatePascalTetrahedron(int order)
 {
   int nbMonomials = (order + 1) * (order + 2) * (order + 3) / 6;
 
-  Double_Matrix monomials(nbMonomials, 3);
+  gmshMatrix<double> monomials(nbMonomials, 3);
 
   int index = 0;
   for (int p = 0; p <= order; p++) {
-    Double_Matrix monOrder = generateMonomialSubspace(3, p);
+    gmshMatrix<double> monOrder = generateMonomialSubspace(3, p);
     int nb = monOrder.size1();
     monomials.copy(monOrder, 0, nb, 0, 3, index, 0);
     index += nb;
@@ -323,14 +323,14 @@ void nodepositionface3(int order,  double *u,  double *v,  double *w)
    }
 }
 
-Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip) 
+gmshMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip) 
 {
   int nbPoints = 
     (serendip ?
      4 +  6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
      (order + 1) * (order + 2) * (order + 3) / 6);
   
-  Double_Matrix point(nbPoints, 3);
+  gmshMatrix<double> point(nbPoints, 3);
 
   double overOrder = (order == 0 ? 1. : 1. / order);
 
@@ -444,7 +444,7 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
         
         if (!serendip && order > 3) {
   
-          Double_Matrix interior = gmshGeneratePointsTetrahedron(order - 4, false);
+          gmshMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
           for (int k = 0; k < interior.size1(); k++) {
             point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
             point(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
@@ -459,10 +459,10 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
   return point;
 }
 
-Double_Matrix gmshGeneratePointsTriangle(int order, bool serendip) 
+gmshMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip) 
 {
   int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
-  Double_Matrix point(nbPoints, 2);
+  gmshMatrix<double> point(nbPoints, 2);
   
   point(0, 0) = 0;
   point(0, 1) = 0;
@@ -507,7 +507,7 @@ Double_Matrix gmshGeneratePointsTriangle(int order, bool serendip)
       } 
 
       if (order > 2 && !serendip) {
-        Double_Matrix inner = gmshGeneratePointsTriangle(order - 3, serendip);
+        gmshMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
         inner.scale(1. - 3. * dd);
         inner.add(dd);
         point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
@@ -517,18 +517,18 @@ Double_Matrix gmshGeneratePointsTriangle(int order, bool serendip)
   return point;  
 }
 
-Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial,
-                                                   const Double_Matrix& point) 
+gmshMatrix<double> generateLagrangeMonomialCoefficients(const gmshMatrix<double>& monomial,
+                                                   const gmshMatrix<double>& point) 
 {
   if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
     Msg::Fatal("Wrong sizes for Lagrange coefficients generation");
-    return Double_Matrix(1, 1);
+    return gmshMatrix<double>(1, 1);
   }
   
   int ndofs = monomial.size1();
   int dim = monomial.size2();
   
-  Double_Matrix Vandermonde(ndofs, ndofs);
+  gmshMatrix<double> Vandermonde(ndofs, ndofs);
   for (int i = 0; i < ndofs; i++) {
     for (int j = 0; j < ndofs; j++) {
       double dd = 1.;
@@ -541,13 +541,13 @@ Double_Matrix generateLagrangeMonomialCoefficients(const Double_Matrix& monomial
 
   if (det == 0.){
     Msg::Fatal("Vandermonde matrix has zero determinant!?");
-    return Double_Matrix(1, 1);
+    return gmshMatrix<double>(1, 1);
   }
-  Double_Matrix coefficient(ndofs, ndofs);
+  gmshMatrix<double> coefficient(ndofs, ndofs);
   for (int i = 0; i < ndofs; i++) {
     for (int j = 0; j < ndofs; j++) {
       int f = (i + j) % 2 == 0 ? 1 : -1;
-      Double_Matrix cofactor = Vandermonde.cofactor(i, j);
+      gmshMatrix<double> cofactor = Vandermonde.cofactor(i, j);
       coefficient(i, j) = f * cofactor.determinant() / det;
     }
   }
@@ -659,18 +659,18 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
   return fs[tag];
 }
 
-std::map<std::pair<int, int>, Double_Matrix> gmshFunctionSpaces::injector;
+std::map<std::pair<int, int>, gmshMatrix<double> > gmshFunctionSpaces::injector;
 
-const Double_Matrix &gmshFunctionSpaces::findInjector(int tag1, int tag2)
+const gmshMatrix<double> &gmshFunctionSpaces::findInjector(int tag1, int tag2)
 {
   std::pair<int,int> key(tag1,tag2);
-  std::map<std::pair<int, int>, Double_Matrix>::const_iterator it = injector.find(key);
+  std::map<std::pair<int, int>, gmshMatrix<double> >::const_iterator it = injector.find(key);
   if (it != injector.end()) return it->second;
 
   const gmshFunctionSpace& fs1 = find(tag1);
   const gmshFunctionSpace& fs2 = find(tag2);
 
-  Double_Matrix inj(fs1.points.size1(),fs2.points.size1());
+  gmshMatrix<double> inj(fs1.points.size1(),fs2.points.size1());
   
   double sf[256];
   
diff --git a/Numeric/FunctionSpace.h b/Numeric/FunctionSpace.h
index 242f023cc930aa8dd9d9b4e9a5d8bf0d11e5cc31..d701cf2fadc5e95dfb2dc88b59b6b28bc3919e78 100644
--- a/Numeric/FunctionSpace.h
+++ b/Numeric/FunctionSpace.h
@@ -12,9 +12,9 @@
 
 struct gmshFunctionSpace 
 {
-  Double_Matrix points;
-  Double_Matrix monomials;
-  Double_Matrix coefficients;
+  gmshMatrix<double> points;
+  gmshMatrix<double> monomials;
+  gmshMatrix<double> coefficients;
   inline void evaluateMonomials(double u, double v, double w, double p[]) const 
   {
     for (int j = 0; j < monomials.size1(); j++) {
@@ -98,10 +98,10 @@ class gmshFunctionSpaces
 {
  private:
   static std::map<int, gmshFunctionSpace> fs;
-  static std::map<std::pair<int, int>, Double_Matrix> injector;
+  static std::map<std::pair<int, int>, gmshMatrix<double> > injector;
  public :
   static const gmshFunctionSpace &find(int);
-  static const Double_Matrix &findInjector(int, int);
+  static const gmshMatrix<double> &findInjector(int, int);
 };
 
 #endif
diff --git a/Numeric/GmshMatrix.cpp b/Numeric/GmshMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0923e0ef85782d439551a9428c1d8e9d13f7022a
--- /dev/null
+++ b/Numeric/GmshMatrix.cpp
@@ -0,0 +1,205 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include <complex>
+#include "GmshConfig.h"
+#include "GmshMatrix.h"
+#include "GmshMessage.h"
+
+#if defined(HAVE_BLAS)
+
+extern "C" {
+  void dgemm_(const char *transa, const char *transb, int *m, int *n, int *k, 
+              double *alpha, double *a, int *lda, 
+              double *b, int *ldb, double *beta, 
+              double *c, int *ldc);
+  void zgemm_(const char *transa, const char *transb, int *m, int *n, int *k, 
+              std::complex<double> *alpha, std::complex<double> *a, int *lda, 
+              std::complex<double> *b, int *ldb, std::complex<double> *beta, 
+              std::complex<double> *c, int *ldc);
+  void dgemv_(const char *trans, int *m, int *n, 
+              double *alpha, double *a, int *lda, 
+              double *x, int *incx, double *beta, 
+              double *y, int *incy);
+  void zgemv_(const char *trans, int *m, int *n, 
+              std::complex<double> *alpha, std::complex<double> *a, int *lda, 
+              std::complex<double> *x, int *incx, std::complex<double> *beta, 
+              std::complex<double> *y, int *incy);
+}
+
+template<> 
+void gmshMatrix<double>::mult(const gmshMatrix<double> &b, gmshMatrix<double> &c)
+{
+  int M = c.size1(), N = c.size2(), K = _c;
+  int LDA = _r, LDB = b.size1(), LDC = c.size1();
+  double alpha = 1., beta = 0.;
+  dgemm_("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB, 
+         &beta, c._data, &LDC);
+}
+
+template<> 
+void gmshMatrix<std::complex<double> >::mult(const gmshMatrix<std::complex<double> > &b, 
+                                             gmshMatrix<std::complex<double> > &c)
+{
+  int M = c.size1(), N = c.size2(), K = _c;
+  int LDA = _r, LDB = b.size1(), LDC = c.size1();
+  std::complex<double> alpha = 1., beta = 0.;
+  zgemm_("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB, 
+         &beta, c._data, &LDC);
+}
+
+template<> 
+void gmshMatrix<double>::gemm(gmshMatrix<double> &a, gmshMatrix<double> &b, 
+                              double alpha, double beta)
+{
+  int M = size1(), N = size2(), K = a.size2();
+  int LDA = a.size1(), LDB = b.size1(), LDC = size1();
+  dgemm_("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB, 
+         &beta, _data, &LDC);
+}
+
+template<> 
+void gmshMatrix<std::complex<double> >::gemm(gmshMatrix<std::complex<double> > &a, 
+                                             gmshMatrix<std::complex<double> > &b, 
+                                             std::complex<double> alpha, 
+                                             std::complex<double> beta)
+{
+  int M = size1(), N = size2(), K = a.size2();
+  int LDA = a.size1(), LDB = b.size1(), LDC = size1();
+  zgemm_("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB, 
+         &beta, _data, &LDC);
+}
+
+template<> 
+void gmshMatrix<double>::mult(const gmshVector<double> &x, gmshVector<double> &y)
+{
+  int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
+  double alpha = 1., beta = 0.;
+  dgemv_("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
+         &beta, y._data, &INCY);
+}
+
+template<> 
+void gmshMatrix<std::complex<double> >::mult(const gmshVector<std::complex<double> > &x, 
+                                             gmshVector<std::complex<double> > &y)
+{
+  int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
+  std::complex<double> alpha = 1., beta = 0.;
+  zgemv_("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
+         &beta, y._data, &INCY);
+}
+
+#else
+
+template<class scalar> 
+void gmshMatrix<scalar>::mult(const gmshMatrix<scalar> &b, gmshMatrix<scalar> &c)
+{
+  c.scale(0.);
+  for(int i = 0; i < _r; i++)
+    for(int j = 0; j < b.size2(); j++)
+      for(int k = 0; k < _c; k++)
+        c._data[i + _r * j] += (*this)(i, k) * b(k, j);
+}
+
+template<class scalar> 
+void gmshMatrix<scalar>::gemm(gmshMatrix<scalar> &a, gmshMatrix<scalar> &b, 
+                              scalar alpha, scalar beta)
+{
+  gmshMatrix<scalar> temp(a.size1(), b.size2());
+  a.mult(b, temp);
+  temp.scale(alpha);
+  scale(beta);
+  add(temp);
+}
+
+template<class scalar> 
+void gmshMatrix<scalar>::mult(const gmshVector<scalar> &x, gmshVector<scalar> &y)
+{
+  y.scale(0.);
+  for(int i = 0; i < _r; i++)
+    for(int j = 0; j < _c; j++)
+      y._data[i] += (*this)(i, j) * x(j);
+}
+
+#endif
+
+#if defined(HAVE_LAPACK)
+
+extern "C" {
+  void dgesv_(int *N, int *nrhs, double *A, int *lda, int *ipiv, 
+              double *b, int *ldb, int *info);
+  void dgetrf_(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
+  void dgesvd_(const char* jobu, const char *jobvt, int *M, int *N,
+               double *A, int *lda, double *S, double* U, int *ldu,
+               double *VT, int *ldvt, double *work, int *lwork, int *info);
+}
+
+template<> 
+bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<double> &result)
+{
+  int N = size1(), nrhs = 1, lda = N, ldb = N, info;
+  int *ipiv = new int[N];
+  for(int i = 0; i < N; i++) result(i) = rhs(i);
+  dgesv_(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
+  delete [] ipiv;
+  if(info == 0) return true;
+  Msg::Error("Problem in LAPACK LU (info=%d)", info);
+  return false;
+}
+
+template<> 
+double gmshMatrix<double>::determinant() const
+{
+  gmshMatrix<double> tmp(*this);
+  int M = size1(), N = size2(), lda = size1(), info;
+  int *ipiv = new int[std::min(M, N)];
+  dgetrf_(&M, &N, tmp._data, &lda, ipiv, &info);
+  double det = 1.;
+  for(int i = 0; i < size1(); i++){
+    det *= tmp(i, i);
+    if(ipiv[i] != i + 1) det = -det;
+  }
+  return det;
+}
+
+template<> 
+bool gmshMatrix<double>::svd(gmshMatrix<double> &V, gmshVector<double> &S)
+{
+  gmshMatrix<double> VT(V.size2(), V.size1());
+  int M = size1(), N = size2(), LDA = size1(), LDVT = VT.size1(), info;
+  int LWORK = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
+  gmshVector<double> WORK(LWORK);
+  dgesvd_("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
+          VT._data, &LDVT, WORK._data, &LWORK, &info);
+  V = VT.transpose();
+  if(info == 0) return true;
+  Msg::Error("Problem in LAPACK SVD (info=%d)", info);
+  return false;
+}
+
+#else
+
+template<class scalar> 
+bool gmshMatrix<scalar>::lu_solve(const gmshVector<scalar> &rhs, gmshVector<scalar> &result)
+{
+  Msg::Error("LU factorization requires LAPACK");
+  return false;
+}
+
+template<class scalar> 
+scalar gmshMatrix<scalar>::determinant() const
+{
+  Msg::Error("Determinant computation requires LAPACK");
+  return 0.;
+}
+
+template<class scalar> 
+bool gmshMatrix<scalar>::svd(gmshMatrix<scalar> &V, gmshVector<scalar> &S)
+{
+  Msg::Error("Singular value decomposition requires LAPACK");
+  return false;
+}
+
+#endif
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index c88554b11f86f2cd578ab75b3d7891e1d08cb58c..1c2ec429e82a4be5df98e3efdeae71723c2d956a 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -7,68 +7,44 @@
 #define _GMSH_MATRIX_H_
 
 #include <math.h>
-#include <algorithm>
-#include "GmshConfig.h"
-#include "GmshMessage.h"
 
-#if defined(HAVE_BLAS)
-extern "C" {
-  void dgemm_(const char *transa, const char *transb, int *m, int *n, int *k, 
-              double *alpha, double *a, int *lda, double *b, int *ldb, 
-              double *beta, double *c, int *ldc);
-  void dgemv_(const char *trans, int *m, int *n, double *alpha, double *a, 
-              int *lda, double *x, int *incx, double *beta, double *y, int *incy);
-}
-#endif
-
-#if defined(HAVE_LAPACK)
-extern "C" {
-  void dgesv_(int *N, int *nrhs, double *A, int *lda, int *ipiv, 
-              double *b, int *ldb, int *info);
-  void dgetrf_(int *M, int *N, double *A, int *lda, int *ipiv, int *info);
-  void dgesvd_(const char* jobu, const char *jobvt, int *M, int *N,
-               double *A, int *lda, double *S, double* U, int *ldu,
-               double *VT, int *ldvt, double *work, int *lwork, int *info);
-}
-#endif
-
-template <class SCALAR> class Gmsh_Matrix;
+template <class scalar> class gmshMatrix;
 
-template <class SCALAR>
-class Gmsh_Vector
+template <class scalar>
+class gmshVector
 {
  private:
   int _r;
-  SCALAR *_data;
-  friend class Gmsh_Matrix<SCALAR>;
+  scalar *_data;
+  friend class gmshMatrix<scalar>;
  public:
-  Gmsh_Vector(int r) : _r(r)
+  gmshVector(int r) : _r(r)
   {
-    _data = new SCALAR[_r];
+    _data = new scalar[_r];
     scale(0.);
   }
-  Gmsh_Vector(const Gmsh_Vector<SCALAR> &other) : _r(other._r)
+  gmshVector(const gmshVector<scalar> &other) : _r(other._r)
   {
-    _data = new SCALAR[_r];
+    _data = new scalar[_r];
     for(int i = 0; i < _r; ++i) _data[i] = other._data[i];
   }
-  ~Gmsh_Vector() { if(_data) delete [] _data; }
-  inline SCALAR operator () (int i) const
+  ~gmshVector() { if(_data) delete [] _data; }
+  inline scalar operator () (int i) const
   {
     return _data[i];
   }
   inline int size() const { return _r; }
-  inline SCALAR & operator () (int i)
+  inline scalar & operator () (int i)
   {
     return _data[i];
   }
-  inline SCALAR norm()
+  inline scalar norm()
   {
-    SCALAR n = 0.;
+    scalar n = 0.;
     for(int i = 0; i < _r; ++i) n += _data[i] * _data[i];
     return sqrt(n);
   }
-  inline void scale(const SCALAR s)
+  inline void scale(const scalar s)
   {
     if(s == 0.) 
       for(int i = 0; i < _r; ++i) _data[i] = 0.;
@@ -77,89 +53,60 @@ class Gmsh_Vector
   }
 };
 
-template <class SCALAR>
-class Gmsh_Matrix
+template <class scalar>
+class gmshMatrix
 {
  private:
   int _r, _c;
-  SCALAR *_data;
+  scalar *_data;
  public:
-  Gmsh_Matrix(int r, int c) : _r(r), _c(c)
+  gmshMatrix(int r, int c) : _r(r), _c(c)
   {
-    _data = new SCALAR[_r * _c];
+    _data = new scalar[_r * _c];
     scale(0.);
   }
-  Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other) : _r(other._r), _c(other._c)
+  gmshMatrix(const gmshMatrix<scalar> &other) : _r(other._r), _c(other._c)
   {
-    _data = new SCALAR[_r * _c];
+    _data = new scalar[_r * _c];
     memcpy(other);
   }
-  Gmsh_Matrix() : _r(0), _c(0), _data(0) {}
-  ~Gmsh_Matrix() { if(_data) delete [] _data; }
+  gmshMatrix() : _r(0), _c(0), _data(0) {}
+  ~gmshMatrix() { if(_data) delete [] _data; }
   inline int size1() const { return _r; }
   inline int size2() const { return _c; }
-  Gmsh_Matrix<SCALAR> & operator = (const Gmsh_Matrix<SCALAR> &other)
+  gmshMatrix<scalar> & operator = (const gmshMatrix<scalar> &other)
   {
     if(this != &other){
       _r = other._r; 
       _c = other._c;
-      _data = new SCALAR[_r * _c];
+      _data = new scalar[_r * _c];
       memcpy(other);
     }
     return *this;
   }
-  void memcpy(const Gmsh_Matrix<SCALAR> &other)
+  void memcpy(const gmshMatrix<scalar> &other)
   {
     for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i];
   }
-  inline SCALAR operator () (int i, int j) const
+  inline scalar operator () (int i, int j) const
   {
     return _data[i + _r * j];
   }
-  inline SCALAR & operator () (int i, int j)
+  inline scalar & operator () (int i, int j)
   {
     return _data[i + _r * j];
   }
-  void copy(const Gmsh_Matrix<SCALAR> &a, int i0, int ni, int j0, int nj, 
+  void copy(const gmshMatrix<scalar> &a, int i0, int ni, int j0, int nj, 
             int desti0, int destj0)
   {
     for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++)
       for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
         (*this)(desti, destj) = a(i, j);
   }
-  inline void mult(const Gmsh_Matrix<SCALAR> &b, Gmsh_Matrix<SCALAR> &c)
-  {
-#if defined(HAVE_BLAS)
-    int M = c.size1(), N = c.size2(), K = _c;
-    int LDA = _r, LDB = b.size1(), LDC = c.size1();
-    double alpha = 1., beta = 0.;
-    dgemm_("N", "N", &M, &N, &K, &alpha, _data, &LDA, b._data, &LDB, 
-           &beta, c._data, &LDC);
-#else
-    c.scale(0.);
-    for(int i = 0; i < _r; i++)
-      for(int j = 0; j < b.size2(); j++)
-	for(int k = 0; k < _c; k++)
-	  c._data[i + _r * j] += (*this)(i, k) * b(k, j);
-#endif
-  }
-  inline void blas_dgemm(Gmsh_Matrix<SCALAR> &a, Gmsh_Matrix<SCALAR> &b, 
-			 SCALAR alpha=1., SCALAR beta=1.)
-  {
-#if defined(HAVE_BLAS)
-    int M = size1(), N = size2(), K = a.size2();
-    int LDA = a.size1(), LDB = b.size1(), LDC = size1();
-    dgemm_("N", "N", &M, &N, &K, &alpha, a._data, &LDA, b._data, &LDB, 
-           &beta, _data, &LDC);
-#else
-    Gmsh_Matrix<SCALAR> temp(a.size1(), b.size2());
-    a.mult(b, temp);
-    temp.scale(alpha);
-    scale(beta);
-    add(temp);
-#endif
-  }
-  inline void set_all(const SCALAR &m) 
+  void mult(const gmshMatrix<scalar> &b, gmshMatrix<scalar> &c);
+  void gemm(gmshMatrix<scalar> &a, gmshMatrix<scalar> &b, 
+            scalar alpha=1., scalar beta=1.);
+  inline void set_all(const scalar &m) 
   {
     for(int i = 0; i < _r * _c; i++) _data[i] = m;
   }
@@ -174,241 +121,27 @@ class Gmsh_Matrix
   {
     for(int i = 0; i < _r * _c; ++i) _data[i] += a;
   }
-  inline void add(const Gmsh_Matrix<SCALAR> &m) 
+  inline void add(const gmshMatrix<scalar> &m) 
   {
     for(int i = 0; i < size1(); i++)
       for(int j = 0; j < size2(); j++)
 	(*this)(i, j) += m(i, j);
   }
-  inline void mult(const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &y)
-  {
-#if defined(HAVE_BLAS)
-    int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
-    double alpha = 1., beta = 0.;
-    dgemv_("N", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
-           &beta, y._data, &INCY);
-#else
-    y.scale(0.);
-    for(int i = 0; i < _r; i++)
-      for(int j = 0; j < _c; j++)
-	y._data[i] += (*this)(i, j) * x(j);
-#endif
-  }
-  inline Gmsh_Matrix<SCALAR> transpose()
-  {
-    Gmsh_Matrix<SCALAR> T(size2(), size1());
-    for(int i = 0; i < size1(); i++)
-      for(int j = 0; j < size2(); j++)
-        T(j, i) = (*this)(i, j);
-    return T;
-  }
-  inline bool lu_solve(const Gmsh_Vector<SCALAR> &rhs, Gmsh_Vector<SCALAR> &result)
-  {
-#if defined(HAVE_LAPACK)
-    int N = size1(), nrhs = 1, lda = N, ldb = N, info;
-    int *ipiv = new int[N];
-    for(int i = 0; i < N; i++) result(i) = rhs(i);
-    dgesv_(&N, &nrhs, _data, &lda, ipiv, result._data, &ldb, &info);
-    delete [] ipiv;
-    if(info == 0) return true;
-    Msg::Error("Problem in LAPACK LU (info=%d)", info);
-#else
-    Msg::Error("LU factorization requires LAPACK");
-#endif
-    return false;
-  }
-  Gmsh_Matrix<SCALAR> cofactor(int i, int j) const 
-  {
-    int ni = size1();
-    int nj = size2();
-    Gmsh_Matrix<SCALAR> cof(ni - 1, nj - 1);
-    for(int I = 0; I < ni; I++){
-      for(int J = 0; J < nj; J++){
-	if(J != j && I != i)
-	  cof(I < i ? I : I - 1, J < j ? J : J - 1) = (*this)(I, J);
-      }
-    }
-    return cof;
-  }
-  SCALAR determinant() const
-  {
-#if defined(HAVE_LAPACK)
-    Gmsh_Matrix<SCALAR> tmp(*this);
-    int M = size1(), N = size2(), lda = size1(), info;
-    int *ipiv = new int[std::min(M, N)];
-    dgetrf_(&M, &N, tmp._data, &lda, ipiv, &info);
-    SCALAR det = 1.;
-    for(int i = 0; i < size1(); i++){
-      det *= tmp(i, i);
-      if(ipiv[i] != i + 1) det = -det;
-    }
-    return det;
-#else
-    Msg::Error("Determinant computation requires LAPACK");
-    return 0.;
-#endif
-  }
-  bool svd(Gmsh_Matrix<SCALAR> &V, Gmsh_Vector<SCALAR> &S)
-  {
-#if defined(HAVE_LAPACK)
-    Gmsh_Matrix<SCALAR> VT(V.size2(), V.size1());
-    int M = size1(), N = size2(), LDA = size1(), LDVT = VT.size1(), info;
-    int LWORK = std::max(3 * std::min(M, N) + std::max(M, N), 5 * std::min(M, N));
-    Gmsh_Vector<SCALAR> WORK(LWORK);
-    dgesvd_("O", "A", &M, &N, _data, &LDA, S._data, _data, &LDA,
-            VT._data, &LDVT, WORK._data, &LWORK, &info);
-    V = VT.transpose();
-    if(info == 0) return true;
-    Msg::Error("Problem in LAPACK SVD (info=%d)", info);
-#else
-    Msg::Error("Singular value decomposition requires LAPACK");
-#endif
-    return false;
-  }
-};
-
-#if defined(HAVE_GSL)
-
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
-
-class GSL_Matrix;
-
-class GSL_Vector
-{
- private:
-  int _r;
-  gsl_vector *_data;
-  friend class GSL_Matrix;
- public:
-  GSL_Vector(int r) : _r(r)
-  {
-    _data = gsl_vector_calloc(_r);
-  }
-  GSL_Vector(const GSL_Vector &other) : _r(other._r)
-  {
-    _data = gsl_vector_calloc(_r);
-    gsl_vector_memcpy(_data, other._data);
-  }
-  ~GSL_Vector() { gsl_vector_free(_data); }
-  inline int size() const { return _r; }
-  inline double operator () (int i) const
-  {
-    return gsl_vector_get(_data, i);
-  }
-  inline double & operator () (int i)
-  {
-    return *gsl_vector_ptr(_data, i);
-  }
-  inline double norm()
-  {
-    return gsl_blas_dnrm2(_data);
-  }
-  inline void scale(const double s)
-  {
-    if(s == 0.) gsl_vector_set_zero(_data);
-    else gsl_vector_scale(_data, s);
-  }
-};
-
-class GSL_Matrix
-{
- private:
-  gsl_matrix *_data;
- public:
-  GSL_Matrix(int r, int  c) { _data = gsl_matrix_calloc(r, c); }
-  GSL_Matrix(const GSL_Matrix &other) : _data(0)
-  {
-    if(_data) gsl_matrix_free(_data);
-    _data = gsl_matrix_calloc(other._data->size1, other._data->size2);
-    gsl_matrix_memcpy(_data, other._data);
-  }
-  GSL_Matrix() : _data(0) {}
-  ~GSL_Matrix() { if(_data && _data->owner == 1) gsl_matrix_free(_data); }
-  inline int size1() const { return _data->size1; }
-  inline int size2() const { return _data->size2; }
-  GSL_Matrix & operator = (const GSL_Matrix &other)
-  {
-    if(&other != this){
-      if(_data) gsl_matrix_free(_data);
-      _data = gsl_matrix_calloc(other._data->size1, other._data->size2);
-      gsl_matrix_memcpy(_data, other._data);
-    }
-    return *this;
-  }
-  void memcpy(const GSL_Matrix &other)
-  {
-    gsl_matrix_memcpy(_data, other._data);
-  }
-  inline double operator () (int i, int j) const
+  void mult(const gmshVector<scalar> &x, gmshVector<scalar> &y);
+  inline gmshMatrix<scalar> transpose()
   {
-    return gsl_matrix_get(_data, i, j);
-  }
-  inline double & operator () (int i, int j)
-  {
-    return *gsl_matrix_ptr(_data, i, j);
-  }
-  void copy(const GSL_Matrix &a, int i0, int ni, int j0, int nj, 
-            int desti0, int destj0)
-  {
-    for(int i = i0, desti = desti0; i < i0 + ni; i++, desti++)
-      for(int j = j0, destj = destj0; j < j0 + nj; j++, destj++)
-        (*this)(desti, destj) = a(i, j);
-  }
-  inline void mult(const GSL_Matrix &b, GSL_Matrix &c)
-  {
-    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., _data, b._data, 0., c._data);
-  }
-  inline void blas_dgemm(GSL_Matrix &a, GSL_Matrix &b, 
-			 double alpha=1., double beta=1.)
-  {      
-    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, alpha, a._data, b._data, beta, _data);
-  }
-  inline void set_all(const double &m) 
-  {
-    gsl_matrix_set_all(_data, m);
-  }
-  inline void scale(const double s) 
-  {
-    if(s == 0.) gsl_matrix_set_zero(_data);
-    else gsl_matrix_scale(_data, s);
-  }
-  inline void add(const double &a) 
-  {
-    gsl_matrix_add_constant(_data, a);
-  }
-  inline void add(const GSL_Matrix &m) 
-  {
-    gsl_matrix_add(_data, m._data);
-  }
-  inline void mult(const GSL_Vector &x, GSL_Vector &y)
-  {
-    gsl_blas_dgemv(CblasNoTrans, 1., _data, x._data, 0., y._data);
-  }
-  inline GSL_Matrix transpose()
-  {
-    GSL_Matrix T(size2(), size1());
+    gmshMatrix<scalar> T(size2(), size1());
     for(int i = 0; i < size1(); i++)
       for(int j = 0; j < size2(); j++)
         T(j, i) = (*this)(i, j);
     return T;
   }
-  inline bool lu_solve(const GSL_Vector &rhs, GSL_Vector &result)
-  {
-    int s;
-    gsl_permutation *p = gsl_permutation_alloc(size1());
-    gsl_linalg_LU_decomp(_data, p, &s);
-    gsl_linalg_LU_solve(_data, p, rhs._data, result._data);
-    gsl_permutation_free(p);
-    return true;
-  }
-  GSL_Matrix cofactor(int i, int j) const 
+  bool lu_solve(const gmshVector<scalar> &rhs, gmshVector<scalar> &result);
+  gmshMatrix<scalar> cofactor(int i, int j) const 
   {
     int ni = size1();
     int nj = size2();
-    GSL_Matrix cof(ni - 1, nj - 1);
+    gmshMatrix<scalar> cof(ni - 1, nj - 1);
     for(int I = 0; I < ni; I++){
       for(int J = 0; J < nj; J++){
 	if(J != j && I != i)
@@ -417,31 +150,8 @@ class GSL_Matrix
     }
     return cof;
   }
-  double determinant() const 
-  {
-    GSL_Matrix tmp = *this;
-    gsl_permutation *p = gsl_permutation_alloc(size1());
-    int s;
-    gsl_linalg_LU_decomp(tmp._data, p, &s);
-    gsl_permutation_free(p);
-    return gsl_linalg_LU_det(tmp._data, s);
-  } 
-  bool svd(GSL_Matrix &V, GSL_Vector &S)
-  {
-    GSL_Vector tmp(S.size());
-    gsl_linalg_SV_decomp(_data, V._data, S._data, tmp._data);
-    return true;
-  }
+  scalar determinant() const;
+  bool svd(gmshMatrix<scalar> &V, gmshVector<scalar> &S);
 };
 
-typedef GSL_Matrix Double_Matrix;
-typedef GSL_Vector Double_Vector;
-
-#else
-
-typedef Gmsh_Matrix<double> Double_Matrix;
-typedef Gmsh_Vector<double> Double_Vector;
-
-#endif
-
 #endif
diff --git a/Numeric/GmshPredicates.h b/Numeric/GmshPredicates.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4f6535de73b96ec3290e6ea25a04b77fc9f32eb
--- /dev/null
+++ b/Numeric/GmshPredicates.h
@@ -0,0 +1,18 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _GMSH_PREDICATES_H_
+#define _GMSH_PREDICATES_H_
+
+// namespace necessary to avoid conflicts with predicates used by Tetgen
+namespace gmsh{
+double exactinit();
+double incircle(double *pa, double *pb, double *pc, double *pd);
+double insphere(double *pa, double *pb, double *pc, double *pd, double *pe);
+double orient2d(double *pa, double *pb, double *pc);
+double orient3d(double *pa, double *pb, double *pc, double *pd);
+}
+
+#endif
diff --git a/Numeric/Makefile b/Numeric/Makefile
index f88949c1e2e8dcb2eba6d0fb1908ad82f98eb599..0f773c137409f3090a38820ace7fb5acfcd526c6 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -12,7 +12,7 @@ INC = ${DASH}I../Common ${DASH}I../Numeric ${DASH}I../Geo
 CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
 
 SRC = Numeric.cpp\
-      NumericEmbedded.cpp\
+      GmshMatrix.cpp\
       gmshAssembler.cpp\
       gmshTermOfFormulation.cpp\
       gmshLaplace.cpp\
@@ -54,9 +54,9 @@ depend:
 
 # DO NOT DELETE THIS LINE
 Numeric${OBJEXT}: Numeric.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
-  Numeric.h NumericEmbedded.h GmshMatrix.h
-NumericEmbedded${OBJEXT}: NumericEmbedded.cpp NumericEmbedded.h GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h
+  Numeric.h GmshMatrix.h
+GmshMatrix${OBJEXT}: GmshMatrix.cpp ../Common/GmshConfig.h GmshMatrix.h \
+  ../Common/GmshMessage.h
 gmshAssembler${OBJEXT}: gmshAssembler.cpp ../Geo/MVertex.h ../Geo/SPoint2.h \
   ../Geo/SPoint3.h gmshAssembler.h gmshLinearSystem.h
 gmshTermOfFormulation${OBJEXT}: gmshTermOfFormulation.cpp ../Common/Gmsh.h \
@@ -76,9 +76,9 @@ gmshTermOfFormulation${OBJEXT}: gmshTermOfFormulation.cpp ../Common/Gmsh.h \
   gmshTermOfFormulation.h gmshFunction.h gmshLinearSystem.h \
   gmshAssembler.h
 gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
-  GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
-  gmshFunction.h ../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  GmshMatrix.h ../Common/GmshConfig.h gmshFunction.h ../Common/Gmsh.h \
+  ../Common/GmshMessage.h ../Geo/GModel.h ../Geo/GVertex.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
   ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
@@ -89,25 +89,25 @@ gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
   ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
   ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
-  ../Numeric/Gauss.h Numeric.h NumericEmbedded.h
+  ../Numeric/Gauss.h Numeric.h
 gmshElasticity${OBJEXT}: gmshElasticity.cpp gmshElasticity.h \
   gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshConfig.h \
-  ../Common/GmshMessage.h ../Common/Gmsh.h ../Common/GmshMessage.h \
-  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
-  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h \
-  ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h \
-  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+  ../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/GModel.h \
+  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
+  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
+  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
   ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
   ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
-  ../Numeric/Gauss.h Numeric.h NumericEmbedded.h
+  ../Numeric/Gauss.h Numeric.h
 EigSolve${OBJEXT}: EigSolve.cpp
 FunctionSpace${OBJEXT}: FunctionSpace.cpp FunctionSpace.h GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/GmshDefines.h
+  ../Common/GmshConfig.h ../Common/GmshDefines.h ../Common/GmshMessage.h
 GaussQuadratureLin${OBJEXT}: GaussQuadratureLin.cpp Gauss.h GaussLegendre1D.h
 GaussQuadratureTri${OBJEXT}: GaussQuadratureTri.cpp Gauss.h GaussLegendre1D.h
 GaussQuadratureQuad${OBJEXT}: GaussQuadratureQuad.cpp Gauss.h GaussLegendre1D.h
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 25b917a560b7e7c8856b71a7b9d811335977433c..335e79ea7e7d4a48a99c1c2324482e62ffff516b 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -7,113 +7,515 @@
 #include "GmshMessage.h"
 #include "Numeric.h"
 
-#if defined(HAVE_GSL)
+#define SQU(a)      ((a)*(a))
 
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_min.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_multiroots.h>
-#include <gsl/gsl_linalg.h>
+double myatan2(double a, double b)
+{
+  if(a == 0.0 && b == 0)
+    return 0.0;
+  return atan2(a, b);
+}
+
+double myasin(double a)
+{
+  if(a <= -1.)
+    return -M_PI / 2.;
+  else if(a >= 1.)
+    return M_PI / 2.;
+  else
+    return asin(a);
+}
 
-void my_gsl_msg(const char *reason, const char *file, int line,
-		int gsl_errno)
+double myacos(double a)
 {
-  Msg::Error("GSL: %s (%s, line %d)", reason, file, line);
+  if(a <= -1.)
+    return M_PI;
+  else if(a >= 1.)
+    return 0.;
+  else
+    return acos(a);
 }
 
-int Init_Numeric()
+void matvec(double mat[3][3], double vec[3], double res[3])
 {
-  // set new error handler
-  gsl_set_error_handler(&my_gsl_msg);
+  res[0] = mat[0][0] * vec[0] + mat[0][1] * vec[1] + mat[0][2] * vec[2];
+  res[1] = mat[1][0] * vec[0] + mat[1][1] * vec[1] + mat[1][2] * vec[2];
+  res[2] = mat[2][0] * vec[0] + mat[2][1] * vec[1] + mat[2][2] * vec[2];
+}
 
-  // initilize robust geometric predicates
-  gmsh::exactinit() ;
-  return 1;
+void normal3points(double x0, double y0, double z0,
+                   double x1, double y1, double z1,
+                   double x2, double y2, double z2,
+                   double n[3])
+{
+  double t1[3] = {x1 - x0, y1 - y0, z1 - z0};
+  double t2[3] = {x2 - x0, y2 - y0, z2 - z0};
+  prodve(t1, t2, n);
+  norme(n);
 }
 
-static double (*f_statN) (double*, void *data); 
+int sys2x2(double mat[2][2], double b[2], double res[2])
+{
+  double det, ud, norm;
+  int i;
 
-static void (*df_statN) (double*, double*, double &, void *data);
+  norm = SQU(mat[0][0]) + SQU(mat[1][1]) + SQU(mat[0][1]) + SQU(mat[1][0]);
+  det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
 
-static double fobjN (const gsl_vector *x, void *data)
+  // TOLERANCE ! WARNING WARNING
+  if(norm == 0.0 || fabs(det) / norm < 1.e-12) {
+    if(norm)
+      Msg::Debug("Assuming 2x2 matrix is singular (det/norm == %.16g)",
+          fabs(det) / norm);
+    res[0] = res[1] = 0.0;
+    return 0;
+  }
+  ud = 1. / det;
+
+  res[0] = b[0] * mat[1][1] - mat[0][1] * b[1];
+  res[1] = mat[0][0] * b[1] - mat[1][0] * b[0];
+
+  for(i = 0; i < 2; i++)
+    res[i] *= ud;
+
+  return (1);
+}
+
+double det3x3(double mat[3][3])
 {
-  return f_statN(x->data, data);
+  return (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
+          mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
+          mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]));
 }
 
-static void dfobjN(const gsl_vector *x, void *params, gsl_vector *g)
+double trace3x3(double mat[3][3])
 {
-  double f;
-  df_statN(x->data, g->data, f, params);
+  return mat[0][0] + mat[1][1] + mat[2][2];
+}
+
+double trace2 (double mat[3][3])
+{
+  double a00 =  mat[0][0] * mat[0][0] + mat[1][0] * mat[0][1] + mat[2][0] * mat[0][2]; 
+  double a11 =  mat[1][0] * mat[0][1] + mat[1][1] * mat[1][1] + mat[1][2] * mat[2][1]; 
+  double a22 =  mat[2][0] * mat[0][2] + mat[2][1] * mat[1][2] + mat[2][2] * mat[2][2];
+  
+  return a00 + a11 + a22;
+}
+
+int sys3x3(double mat[3][3], double b[3], double res[3], double *det)
+{
+  double ud;
+  int i;
+
+  *det = det3x3(mat);
+
+  if(*det == 0.0) {
+    res[0] = res[1] = res[2] = 0.0;
+    return (0);
+  }
+
+  ud = 1. / (*det);
+
+  res[0] = b[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
+    mat[0][1] * (b[1] * mat[2][2] - mat[1][2] * b[2]) +
+    mat[0][2] * (b[1] * mat[2][1] - mat[1][1] * b[2]);
+
+  res[1] = mat[0][0] * (b[1] * mat[2][2] - mat[1][2] * b[2]) -
+    b[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
+    mat[0][2] * (mat[1][0] * b[2] - b[1] * mat[2][0]);
+
+  res[2] = mat[0][0] * (mat[1][1] * b[2] - b[1] * mat[2][1]) -
+    mat[0][1] * (mat[1][0] * b[2] - b[1] * mat[2][0]) +
+    b[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
+
+  for(i = 0; i < 3; i++)
+    res[i] *= ud;
+  return (1);
 }
 
-static void fdfobjN(const gsl_vector *x, void *params, double *f, gsl_vector *g)
+int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
 {
-  df_statN(x->data, g->data, *f, params);
+  int out;
+  double norm;
+
+  out = sys3x3(mat, b, res, det);
+  norm =
+    SQU(mat[0][0]) + SQU(mat[0][1]) + SQU(mat[0][2]) +
+    SQU(mat[1][0]) + SQU(mat[1][1]) + SQU(mat[1][2]) +
+    SQU(mat[2][0]) + SQU(mat[2][1]) + SQU(mat[2][2]);
+
+  // TOLERANCE ! WARNING WARNING
+  if(norm == 0.0 || fabs(*det) / norm < 1.e-12) {
+    if(norm)
+      Msg::Debug("Assuming 3x3 matrix is singular (det/norm == %.16g)",
+          fabs(*det) / norm);
+    res[0] = res[1] = res[2] = 0.0;
+    return 0;
+  }
+
+  return out;
 }
 
-void minimize_N(int N, double (*f)(double*, void *data), 
-                void (*df)(double*, double*, double &, void *data),
-                void *data, int niter, double *U, double &res)
+double det2x2(double mat[2][2])
 {
-  f_statN = f;
-  df_statN = df;
+  return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
+}
 
-  int iter = 0;
-  int status;
+double det2x3(double mat[2][3])
+{
+  double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
+  double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
+  double n[3];
   
-  const gsl_multimin_fdfminimizer_type *T;
-  gsl_multimin_fdfminimizer *s;
-  gsl_vector *x;
-  gsl_multimin_function_fdf my_func;
+  prodve(v1, v2, n);
+  return norm3(n);
+}
 
-  my_func.f = &fobjN;
-  my_func.df = &dfobjN;
-  my_func.fdf = &fdfobjN;
-  my_func.n = N;
-  my_func.params = data;
+double inv2x2(double mat[2][2], double inv[2][2])
+{
+  const double det = det2x2(mat);
+  if(det){
+    double ud = 1. / det;
+    inv[0][0] =  mat[1][1] * ud;
+    inv[1][0] = -mat[1][0] * ud;
+    inv[0][1] = -mat[0][1] * ud;
+    inv[1][1] =  mat[0][0] * ud;
+  }
+  else{
+    Msg::Error("Singular matrix");
+    for(int i = 0; i < 2; i++)
+      for(int j = 0; j < 2; j++)
+        inv[i][j] = 0.;
+  }
+  return det;
+}
 
-  x = gsl_vector_alloc(N);
-  for (int i = 0; i < N; i++)
-    gsl_vector_set(x, i, U[i]);
+double inv3x3(double mat[3][3], double inv[3][3])
+{
+  double det = det3x3(mat);
+  if(det){
+    double ud = 1. / det;
+    inv[0][0] =  (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) * ud;
+    inv[1][0] = -(mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) * ud;
+    inv[2][0] =  (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) * ud;
+    inv[0][1] = -(mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) * ud;
+    inv[1][1] =  (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) * ud;
+    inv[2][1] = -(mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) * ud;
+    inv[0][2] =  (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]) * ud;
+    inv[1][2] = -(mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]) * ud;
+    inv[2][2] =  (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]) * ud;
+  }
+  else{
+    Msg::Error("Singular matrix");
+    for(int i = 0; i < 3; i++)
+      for(int j = 0; j < 3; j++)
+        inv[i][j] = 0.;
+  }
+  return det;
+}
 
-  T = gsl_multimin_fdfminimizer_conjugate_fr;
-  s = gsl_multimin_fdfminimizer_alloc(T, N);
+double angle_02pi(double A3)
+{
+  double DP = 2 * M_PI;
+  while(A3 > DP || A3 < 0.) {
+    if(A3 > 0)
+      A3 -= DP;
+    else
+      A3 += DP;
+  }
+  return A3;
+}
 
-  gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 1e-4);
+double angle_plan(double V[3], double P1[3], double P2[3], double n[3])
+{
+  double PA[3], PB[3], angplan;
+  double cosc, sinc, c[3];
 
-  do{
-    iter++;
-    status = gsl_multimin_fdfminimizer_iterate(s);
-    
-    if (status)
-      break;
-    
-    status = gsl_multimin_test_gradient(s->gradient, 1e-3);
+  PA[0] = P1[0] - V[0];
+  PA[1] = P1[1] - V[1];
+  PA[2] = P1[2] - V[2];
+
+  PB[0] = P2[0] - V[0];
+  PB[1] = P2[1] - V[1];
+  PB[2] = P2[2] - V[2];
+
+  norme(PA);
+  norme(PB);
+
+  prodve(PA, PB, c);
+
+  prosca(PA, PB, &cosc);
+  prosca(c, n, &sinc);
+  angplan = myatan2(sinc, cosc);
+
+  return angplan;
+}
+
+double triangle_area(double p0[3], double p1[3], double p2[3])
+{
+  double a[3], b[3], c[3];
+  
+  a[0] = p2[0] - p1[0];
+  a[1] = p2[1] - p1[1];
+  a[2] = p2[2] - p1[2];
+  
+  b[0] = p0[0] - p1[0];
+  b[1] = p0[1] - p1[1];
+  b[2] = p0[2] - p1[2];
+  
+  prodve(a, b, c);
+  return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
+}
+
+char float2char(float f)
+{
+  // float normalized in [-1, 1], char in [-127, 127]
+  f *= 127.;
+  if(f > 127.) return 127;
+  else if(f < -127.) return -127;
+  else return (char)f;
+}
+
+float char2float(char c)
+{
+  float f = c;
+  f /= 127.;
+  if(f > 1.) return 1.;
+  else if(f < -1) return -1.;
+  else return f;
+}
+
+double InterpolateIso(double *X, double *Y, double *Z,
+                      double *Val, double V, int I1, int I2,
+                      double *XI, double *YI, double *ZI)
+{
+  if(Val[I1] == Val[I2]) {
+    *XI = X[I1];
+    *YI = Y[I1];
+    *ZI = Z[I1];
+    return 0;
+  }
+  else {
+    double coef = (V - Val[I1]) / (Val[I2] - Val[I1]);
+    *XI = coef * (X[I2] - X[I1]) + X[I1];
+    *YI = coef * (Y[I2] - Y[I1]) + Y[I1];
+    *ZI = coef * (Z[I2] - Z[I1]) + Z[I1];
+    return coef;
+  }
+}
+
+void gradSimplex(double *x, double *y, double *z, double *v, double *grad)
+{
+  // p = p1 * (1-u-v-w) + p2 u + p3 v + p4 w
+
+  double mat[3][3];
+  double det, b[3];
+  mat[0][0] = x[1] - x[0];
+  mat[1][0] = x[2] - x[0];
+  mat[2][0] = x[3] - x[0];
+  mat[0][1] = y[1] - y[0];
+  mat[1][1] = y[2] - y[0];
+  mat[2][1] = y[3] - y[0];
+  mat[0][2] = z[1] - z[0];
+  mat[1][2] = z[2] - z[0];
+  mat[2][2] = z[3] - z[0];
+  b[0] = v[1] - v[0];
+  b[1] = v[2] - v[0];
+  b[2] = v[3] - v[0];
+  sys3x3(mat, b, grad, &det);
+}
+
+double ComputeVonMises(double *V)
+{
+  double tr = (V[0] + V[4] + V[8]) / 3.;
+  double v11 = V[0] - tr, v12 = V[1],      v13 = V[2];
+  double v21 = V[3],      v22 = V[4] - tr, v23 = V[5];
+  double v31 = V[6],      v32 = V[7],      v33 = V[8] - tr;
+  return sqrt(1.5 * (v11 * v11 + v12 * v12 + v13 * v13 +
+                     v21 * v21 + v22 * v22 + v23 * v23 +
+                     v31 * v31 + v32 * v32 + v33 * v33));
+}
+
+double ComputeScalarRep(int numComp, double *val)
+{
+  if(numComp == 1)
+    return val[0];
+  else if(numComp == 3)
+    return sqrt(val[0] * val[0] + val[1] * val[1] + val[2] * val[2]);
+  else if(numComp == 9)
+    return ComputeVonMises(val);
+  return 0.;
+}
+
+void eigenvalue(double mat[3][3], double v[3])
+{            
+  // characteristic polynomial of T : find v root of
+  // v^3 - I1 v^2 + I2 T - I3 = 0
+  // I1 : first invariant , trace(T)
+  // I2 : second invariant , 1/2 (I1^2 -trace(T^2))
+  // I3 : third invariant , det T
+
+  double c[4];
+  c[3] = 1.0;
+  c[2] = - trace3x3(mat);
+  c[1] = 0.5 * (c[2] * c[2] - trace2(mat));
+  c[0] = - det3x3(mat);
+  
+  //printf("%g %g %g\n", mat[0][0], mat[0][1], mat[0][2]);
+  //printf("%g %g %g\n", mat[1][0], mat[1][1], mat[1][2]);
+  //printf("%g %g %g\n", mat[2][0], mat[2][1], mat[2][2]);
+  //printf("%g x^3 + %g x^2 + %g x + %g = 0\n", c[3], c[2], c[1], c[0]);
+  
+  double imag[3];
+  FindCubicRoots(c, v, imag);
+  eigsort(v);
+}
+
+void FindCubicRoots(const double coef[4], double real[3], double imag[3])
+{
+  double a = coef[3];
+  double b = coef[2];
+  double c = coef[1];
+  double d = coef[0];
+
+  if(!a || !d){
+    Msg::Error("Degenerate cubic: use a second degree solver!");
+    return;
   }
-  while (status == GSL_CONTINUE && iter < niter);
+
+  b /= a;
+  c /= a;
+  d /= a;
   
-  for (int i = 0; i < N; i++)
-    U[i] = gsl_vector_get(s->x, i);
-  res = s->f;
-  gsl_multimin_fdfminimizer_free(s);
-  gsl_vector_free(x);
-}                                           
+  double q = (3.0*c - (b*b))/9.0;
+  double r = -(27.0*d) + b*(9.0*c - 2.0*(b*b));
+  r /= 54.0;
+
+  double discrim = q*q*q + r*r;
+  imag[0] = 0.0; // The first root is always real.
+  double term1 = (b/3.0);
 
-#else
+  if (discrim > 0) { // one root is real, two are complex
+    double s = r + sqrt(discrim);
+    s = ((s < 0) ? -pow(-s, (1.0/3.0)) : pow(s, (1.0/3.0)));
+    double t = r - sqrt(discrim);
+    t = ((t < 0) ? -pow(-t, (1.0/3.0)) : pow(t, (1.0/3.0)));
+    real[0] = -term1 + s + t;
+    term1 += (s + t)/2.0;
+    real[1] = real[2] = -term1;
+    term1 = sqrt(3.0)*(-t + s)/2;
+    imag[1] = term1;
+    imag[2] = -term1;
+    return;
+  }
 
-int Init_Numeric()
+  // The remaining options are all real
+  imag[1] = imag[2] = 0.0;
+  
+  double r13;
+  if (discrim == 0){ // All roots real, at least two are equal.
+    r13 = ((r < 0) ? -pow(-r,(1.0/3.0)) : pow(r,(1.0/3.0)));
+    real[0] = -term1 + 2.0*r13;
+    real[1] = real[2] = -(r13 + term1);
+    return;
+  }
+
+  // Only option left is that all roots are real and unequal (to get
+  // here, q < 0)
+  q = -q;
+  double dum1 = q*q*q;
+  dum1 = acos(r/sqrt(dum1));
+  r13 = 2.0*sqrt(q);
+  real[0] = -term1 + r13*cos(dum1/3.0);
+  real[1] = -term1 + r13*cos((dum1 + 2.0*M_PI)/3.0);
+  real[2] = -term1 + r13*cos((dum1 + 4.0*M_PI)/3.0);
+}
+
+void eigsort(double d[3])
 {
-  gmsh::exactinit() ;
-  return 1;
+  int k, j, i;
+  double p;
+  
+  for (i=0; i<3; i++) {
+    p=d[k=i];
+    for (j=i+1;j<3;j++)
+      if (d[j] >= p) p=d[k=j];
+    if (k != i) {
+      d[k]=d[i];
+      d[i]=p;
+    }
+  }
 }
 
-void minimize_N(int N, double (*f) (double*, void *data), 
-                void (*df) (double*, double*, double &, void *data) ,
-                void *data,int niter,
-                double *, double &res)
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3])
 {
-  Msg::Error("Gmsh must be compiled with GSL support for minimize_N");
+  int i, j, k, n = 3;
+  double TT[3][3];
+
+  for(i = 1; i <= n; i++) {
+    for(j = 1; j <= n; j++) {
+      II[i - 1][j - 1] = 0.0;
+      TT[i - 1][j - 1] = 0.0;
+    }
+  }
+
+  gmshMatrix<double> M(3, 3), V(3, 3);
+  gmshVector<double> W(3);
+  for(i = 1; i <= n; i++){
+    for(j = 1; j <= n; j++){
+      M(i - 1, j - 1) = MM[i - 1][j - 1];
+    }
+  }
+  M.svd(V, W);
+  for(i = 1; i <= n; i++) {
+    for(j = 1; j <= n; j++) {
+      double ww = W(i - 1);
+      if(fabs(ww) > 1.e-16) { // singular value precision
+        TT[i - 1][j - 1] += M(j - 1, i - 1) / ww;
+      }
+    }
+  }
+  for(i = 1; i <= n; i++) {
+    for(j = 1; j <= n; j++) {
+      for(k = 1; k <= n; k++) {
+        II[i - 1][j - 1] += V(i - 1, k - 1) * TT[k - 1][j - 1];
+      }
+    }
+  }
 }
 
-#endif
+bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
+               gmshVector<double> &x, void *data, double relax, double tolx)
+{
+  const int MAXIT = 50;
+  const double EPS = 1.e-4;
+  const int N = x.size();
+  
+  gmshMatrix<double> J(N, N);
+  gmshVector<double> f(N), feps(N), dx(N);
+  
+  for (int iter = 0; iter < MAXIT; iter++){
+    func(x, f, data);
+
+    for (int j = 0; j < N; j++){
+      double h = EPS * fabs(x(j));
+      if(h == 0.) h = EPS;
+      x(j) += h;
+      func(x, feps, data);
+      for (int i = 0; i < N; i++)
+        J(i, j) = (feps(i) - f(i)) / h;
+      x(j) -= h;
+    }
+    
+    if (N == 1)
+      dx(0) = f(0) / J(0, 0);
+    else
+      J.lu_solve(f, dx);
+    
+    for (int i = 0; i < N; i++)
+      x(i) -= relax * dx(i);
+
+    if(dx.norm() < tolx) return true; 
+  }
+  return false;
+}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 59573e793d605ea5cf53883eb9e065a2a3235d48..aeeeb0b39dd2ac737e793d92558aca1d32ecb6d9 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -6,21 +6,71 @@
 #ifndef _NUMERIC_H_
 #define _NUMERIC_H_
 
-#include "NumericEmbedded.h"
+#include <math.h>
+#include "GmshMatrix.h"
 
-int Init_Numeric();
-void minimize_N(int N, double (*f) (double*, void *data), 
-                void (*df) (double*, double*, double &, void *data),
-                void *data,int niter,
-                double *, double &res);
+#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
+#define sign(x)      (((x)>=0)?1:-1)
 
-// Robust geometrical predicates
-namespace gmsh{
-double exactinit();
-double incircle(double *pa, double *pb, double *pc, double *pd);
-double insphere(double *pa, double *pb, double *pc, double *pd, double *pe);
-double orient2d(double *pa, double *pb, double *pc);
-double orient3d(double *pa, double *pb, double *pc, double *pd);
+double myatan2(double a, double b);
+double myasin(double a);
+double myacos(double a);
+inline void prodve(double a[3], double b[3], double c[3])
+{
+  c[2] = a[0] * b[1] - a[1] * b[0];
+  c[1] = -a[0] * b[2] + a[2] * b[0];
+  c[0] = a[1] * b[2] - a[2] * b[1];
 }
+inline void prosca(double a[3], double b[3], double *c)
+{
+  *c = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+void matvec(double mat[3][3], double vec[3], double res[3]);
+inline double norm3(double a[3])
+{
+  return sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
+}
+inline double norme(double a[3])
+{
+  const double mod = norm3(a);
+  if(mod != 0.0){
+    const double one_over_mod = 1./mod;
+    a[0] *= one_over_mod;
+    a[1] *= one_over_mod;
+    a[2] *= one_over_mod;
+  }
+  return mod;
+}
+void normal3points(double x0, double y0, double z0,
+                   double x1, double y1, double z1,
+                   double x2, double y2, double z2,
+                   double n[3]);
+int sys2x2(double mat[2][2], double b[2], double res[2]);
+int sys3x3(double mat[3][3], double b[3], double res[3], double *det);
+int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det);
+double det2x2(double mat[2][2]);
+double det2x3(double mat[2][3]);
+double det3x3(double mat[3][3]);
+double trace3x3(double mat[3][3]);
+double trace2 (double mat[3][3]);
+double inv3x3(double mat[3][3], double inv[3][3]);
+double inv2x2(double mat[2][2], double inv[2][2]);
+double angle_02pi(double A3);
+double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
+double triangle_area(double p0[3], double p1[3], double p2[3]);
+char float2char(float f);
+float char2float(char c);
+void eigenvalue(double mat[3][3], double re[3]);
+void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
+void eigsort(double d[3]);
+double InterpolateIso(double *X, double *Y, double *Z, 
+                      double *Val, double V, int I1, int I2, 
+                      double *XI, double *YI ,double *ZI);
+void gradSimplex(double *x, double *y, double *z, double *v, double *grad);
+double ComputeVonMises(double *val);
+double ComputeScalarRep(int numComp, double *val);
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
+bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
+               gmshVector<double> &x, void *data, double relax=1., double tolx=1.e-6);
 
 #endif
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index fa1a54e7183ecd89c5e6cc2c15d6542a300dbc2a..79f95453fe9255b403ae1a61c0d119af2f637579 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -6,7 +6,7 @@
 #include "gmshElasticity.h"
 #include "Numeric.h"
 
-void gmshElasticityTerm::elementMatrix(MElement *e, Double_Matrix & m) const
+void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
 {
   int nbNodes = e->getNumVertices();
   int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
@@ -31,10 +31,10 @@ void gmshElasticityTerm::elementMatrix(MElement *e, Double_Matrix & m) const
       {  0,   0,   0,    0, C44,   0}, 
       {  0,   0,   0,    0,   0, C44} };
   
-  Double_Matrix H(6, 6);
-  Double_Matrix B(6, 3 * nbNodes);
-  Double_Matrix BTH(3 * nbNodes, 6);
-  Double_Matrix BT(3 * nbNodes, 6);
+  gmshMatrix<double> H(6, 6);
+  gmshMatrix<double> B(6, 3 * nbNodes);
+  gmshMatrix<double> BTH(3 * nbNodes, 6);
+  gmshMatrix<double> BT(3 * nbNodes, 6);
   for (int i = 0; i < 6; i++)
     for (int j = 0; j < 6; j++)
       H(i, j) = C[i][j];
@@ -70,12 +70,7 @@ void gmshElasticityTerm::elementMatrix(MElement *e, Double_Matrix & m) const
       BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
     }
     BTH.set_all(0.);
-    BTH.blas_dgemm(BT, H); 
-    m.blas_dgemm(BTH, B, weight * detJ, 1.);
+    BTH.gemm(BT, H); 
+    m.gemm(BTH, B, weight * detJ, 1.);
   } 
 }
-
-
-
-
-
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index a3953991a91fc19ec7c78e5f79968a9fb80fc6be..a7dc483e049e1e28aee065b523638f9c8ac1b80c 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -28,7 +28,7 @@ class gmshElasticityTerm : public gmshNodalFemTerm {
  public:
   gmshElasticityTerm(GModel *gm, double E, double nu, int iField = 1) : 
     gmshNodalFemTerm(gm), _E(E), _nu(nu), _iField(iField){}
-  void elementMatrix(MElement *e, Double_Matrix & m) const;
+  void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
 };
 
 #endif
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index d47a79ed27b3e0541fdda74e48b78458c5e3557c..0c9f54af5412e889360fddec534e8b44f83b9a86 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -6,7 +6,7 @@
 #include "gmshLaplace.h"
 #include "Numeric.h"
 
-void gmshLaplaceTerm::elementMatrix(MElement *e, Double_Matrix & m) const
+void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
 {
   int nbNodes = e->getNumVertices();
   int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index 7b56095159fc44e1ada705612fd82fe910a22077..d1ab8d90e148908ed61450ef05006c71aad94aaa 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -29,7 +29,7 @@ class gmshLaplaceTerm : public gmshNodalFemTerm {
  public:
   gmshLaplaceTerm(GModel *gm, gmshFunction *diffusivity, int iField = 0) : 
     gmshNodalFemTerm(gm), _diffusivity(diffusivity), _iField(iField){}
-  virtual void elementMatrix(MElement *e, Double_Matrix &m) const;
+  virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
 };
 
 #endif
diff --git a/Numeric/gmshLinearSystemFull.h b/Numeric/gmshLinearSystemFull.h
index 9302e4e8e047a18fa69ec841194dd793cc00becc..054632dcd6d304ac6dcb685f84f5338f9d101691 100644
--- a/Numeric/gmshLinearSystemFull.h
+++ b/Numeric/gmshLinearSystemFull.h
@@ -14,8 +14,8 @@
 #include "GmshMatrix.h"
 
 class gmshLinearSystemFull : public gmshLinearSystem {
-  Double_Matrix *_a;
-  Double_Vector *_b, *_x;
+  gmshMatrix<double> *_a;
+  gmshVector<double> *_b, *_x;
 public :
   gmshLinearSystemFull () : _a(0), _b(0), _x(0){}
   virtual bool isAllocated () const {return _a != 0;}
@@ -24,9 +24,9 @@ public :
     if (_a) delete _a;
     if (_x) delete _x;
     if (_b) delete _b;
-    _a = new  Double_Matrix(_nbRows,_nbRows);
-    _b = new  Double_Vector(_nbRows);
-    _x = new  Double_Vector(_nbRows);    
+    _a = new  gmshMatrix<double>(_nbRows,_nbRows);
+    _b = new  gmshVector<double>(_nbRows);
+    _x = new  gmshVector<double>(_nbRows);    
   }
   virtual ~gmshLinearSystemFull ()
   {
diff --git a/Numeric/gmshTermOfFormulation.cpp b/Numeric/gmshTermOfFormulation.cpp
index 088433923f1e1c4af6a6aa6dd0ea8efecd0be15e..554b6b1b48157f68ba51c89e9a5d86aa487434b5 100644
--- a/Numeric/gmshTermOfFormulation.cpp
+++ b/Numeric/gmshTermOfFormulation.cpp
@@ -55,7 +55,7 @@ void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys,const std::vector<MElemen
 }
 
 void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, 
-                                   Double_Matrix &localMatrix, 
+                                   gmshMatrix<double> &localMatrix, 
                                    MElement *e) const
 {
   const int nbR = sizeOfR(e);
@@ -78,7 +78,7 @@ void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, MElement *e) const
 {
   const int nbR = sizeOfR(e);
   const int nbC = sizeOfC(e);
-  Double_Matrix localMatrix (nbR, nbC);
+  gmshMatrix<double> localMatrix (nbR, nbC);
   elementMatrix(e, localMatrix);
   addToMatrix(lsys, localMatrix, e);
 }
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index be89cd1ec3f9ce4ac25d436e03c0139ebc81cf09..ffd3cee46af07b836e0e147d90a9cf3e04f8ef20 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -51,13 +51,13 @@ class gmshNodalFemTerm : public gmshTermOfFormulation {
   gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation(gm) {}
   virtual ~gmshNodalFemTerm ();
   // compute the element matrix
-  virtual void elementMatrix(MElement *e, Double_Matrix &m) const = 0;
+  virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const = 0;
 
   void addToMatrix(gmshAssembler &J, MElement *e) const;
   void addToMatrix(gmshAssembler &J, GEntity *ge) const;
   void addToMatrix(gmshAssembler &J) const;
   void addToMatrix(gmshAssembler &J,const std::vector<MElement*> &) const;
-  void addToMatrix(gmshAssembler &Jac, Double_Matrix &localMatrix, MElement *e) const;
+  void addToMatrix(gmshAssembler &Jac, gmshMatrix<double> &localMatrix, MElement *e) const;
 
   void addDirichlet(int physical, int dim, int comp, int field, const gmshFunction &e, 
                     gmshAssembler &);
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 7099ca209f78f02048768bd70213db3ce8e73823..8077401af6406a6a47e148044f838ab2b23b10d3 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -386,7 +386,7 @@ void yyerror(char *s);
 void yymsg(int level, const char *fmt, ...);
 void skip_until(const char *skip, const char *until);
 int PrintListOfDouble(char *format, List_T *list, char *buffer);
-Double_Matrix ListOfListOfDouble2Matrix(List_T *list);
+gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list);
 void FixRelativePath(const char *in, char *out);
 
 
@@ -8097,7 +8097,7 @@ int PrintListOfDouble(char *format, List_T *list, char *buffer)
   return 0;
 }
 
-Double_Matrix ListOfListOfDouble2Matrix(List_T *list)
+gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list)
 {
   int M = List_Nbr(list);
   int N = 0;
@@ -8105,7 +8105,7 @@ Double_Matrix ListOfListOfDouble2Matrix(List_T *list)
     List_T *line = *(List_T**)List_Pointer_Fast(list, i);
     N = std::max(N, List_Nbr(line));
   }
-  Double_Matrix mat(M, N);
+  gmshMatrix<double> mat(M, N);
   for(int i = 0; i < M; i++){
     List_T *line = *(List_T**)List_Pointer_Fast(list, i);
     for(int j = 0; j < List_Nbr(line); j++){
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 8764084d9c517b1f0ded91093c5373ea4385820d..7f5fdeb80ea5097d6d48ea829447f7808516dbdc 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -65,7 +65,7 @@ void yyerror(char *s);
 void yymsg(int level, const char *fmt, ...);
 void skip_until(const char *skip, const char *until);
 int PrintListOfDouble(char *format, List_T *list, char *buffer);
-Double_Matrix ListOfListOfDouble2Matrix(List_T *list);
+gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list);
 void FixRelativePath(const char *in, char *out);
 %}
 
@@ -3503,7 +3503,7 @@ int PrintListOfDouble(char *format, List_T *list, char *buffer)
   return 0;
 }
 
-Double_Matrix ListOfListOfDouble2Matrix(List_T *list)
+gmshMatrix<double> ListOfListOfDouble2Matrix(List_T *list)
 {
   int M = List_Nbr(list);
   int N = 0;
@@ -3511,7 +3511,7 @@ Double_Matrix ListOfListOfDouble2Matrix(List_T *list)
     List_T *line = *(List_T**)List_Pointer_Fast(list, i);
     N = std::max(N, List_Nbr(line));
   }
-  Double_Matrix mat(M, N);
+  gmshMatrix<double> mat(M, N);
   for(int i = 0; i < M; i++){
     List_T *line = *(List_T**)List_Pointer_Fast(list, i);
     for(int j = 0; j < List_Nbr(line); j++){
diff --git a/Parser/Makefile b/Parser/Makefile
index 1d449f731f33901ff4ab402ed1b2fdc0afbc4234..446dfe2d549a2739e868449166e843f6f069ceb1 100644
--- a/Parser/Makefile
+++ b/Parser/Makefile
@@ -56,10 +56,9 @@ depend:
 Gmsh.tab${OBJEXT}: Gmsh.tab.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   ../Numeric/GmshMatrix.h ../Common/MallocUtils.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h ../Geo/GModel.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h ../Geo/GModel.h \
+  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
   ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
@@ -81,9 +80,8 @@ Gmsh.yy${OBJEXT}: Gmsh.yy.cpp ../Common/GmshMessage.h ../Geo/Geo.h \
   ../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
   ../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
-  ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
-  ../Geo/SPoint2.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
-  Gmsh.tab.hpp
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
+  ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
+  ../Common/ListUtils.h ../Geo/SPoint2.h ../Geo/ExtrudeParams.h \
+  ../Common/SmoothData.h Gmsh.tab.hpp
 FunctionManager${OBJEXT}: FunctionManager.cpp FunctionManager.h
diff --git a/Plugin/CutPlane.cpp b/Plugin/CutPlane.cpp
index 4aac591555ccbfe661424bd45a51fb140d62e2d6..b1cca992e1ca8198ac7cd66ae536941709d6fb94 100644
--- a/Plugin/CutPlane.cpp
+++ b/Plugin/CutPlane.cpp
@@ -160,7 +160,7 @@ double GMSH_CutPlanePlugin::levelset(double x, double y, double z, double val) c
     CutPlaneOptions_Number[2].def * z + CutPlaneOptions_Number[3].def;
 }
 
-bool GMSH_CutPlanePlugin::geometricalFilter(Double_Matrix *node_positions) const
+bool GMSH_CutPlanePlugin::geometricalFilter(gmshMatrix<double> *node_positions) const
 {
   const double l0 = levelset((*node_positions)(0, 0),
                              (*node_positions)(0, 1),
diff --git a/Plugin/CutPlane.h b/Plugin/CutPlane.h
index 0df33a8d9690a44b07bfa949ebbcbe8c95c64125..afe0f7d19ba2a346f5593716de295db9fa004499 100644
--- a/Plugin/CutPlane.h
+++ b/Plugin/CutPlane.h
@@ -27,7 +27,7 @@ public:
   int getNbOptions() const;
   StringXNumber *getOption(int iopt);  
   PView *execute(PView *);
-  virtual bool geometricalFilter(Double_Matrix *) const;
+  virtual bool geometricalFilter(gmshMatrix<double> *) const;
 
   static double callbackA(int, int, double);
   static double callbackB(int, int, double);
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 736f01abce6f9473e559838570d90eafefd8b6c6..0cd3036a2d39d4adbe2fca9bff6f5c14052c6ff6 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -84,8 +84,7 @@ Levelset${OBJEXT}: Levelset.cpp Levelset.h Plugin.h ../Common/Options.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/ListUtils.h MakeSimplex.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Post/adaptiveData.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Post/adaptiveData.h
 CutPlane${OBJEXT}: CutPlane.cpp ../Common/GmshConfig.h CutPlane.h Levelset.h \
   Plugin.h ../Common/Options.h ../Post/ColorTable.h \
   ../Common/GmshMessage.h ../Post/PView.h ../Geo/SPoint3.h \
@@ -124,20 +123,19 @@ Lambda2${OBJEXT}: Lambda2.cpp Lambda2.h Plugin.h ../Common/Options.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/ListUtils.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 Eigenvectors${OBJEXT}: Eigenvectors.cpp Eigenvectors.h Plugin.h \
   ../Common/Options.h ../Post/ColorTable.h ../Common/GmshMessage.h \
   ../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
   ../Post/PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Numeric/EigSolve.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Numeric/EigSolve.h
 Eigenvalues${OBJEXT}: Eigenvalues.cpp Eigenvalues.h Plugin.h ../Common/Options.h \
   ../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/ListUtils.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 StreamLines${OBJEXT}: StreamLines.cpp ../Common/GmshConfig.h StreamLines.h \
   Plugin.h ../Common/Options.h ../Post/ColorTable.h \
   ../Common/GmshMessage.h ../Post/PView.h ../Geo/SPoint3.h \
@@ -186,14 +184,13 @@ Warp${OBJEXT}: Warp.cpp Warp.h Plugin.h ../Common/Options.h ../Post/ColorTable.h
   ../Post/PViewDataList.h ../Post/PViewData.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
   ../Common/ListUtils.h ../Common/SmoothData.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 SphericalRaise${OBJEXT}: SphericalRaise.cpp SphericalRaise.h Plugin.h \
   ../Common/Options.h ../Post/ColorTable.h ../Common/GmshMessage.h \
   ../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
   ../Post/PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
 Skin${OBJEXT}: Skin.cpp Skin.h Plugin.h ../Common/Options.h ../Post/ColorTable.h \
   ../Common/GmshMessage.h ../Post/PView.h ../Geo/SPoint3.h \
   ../Post/PViewDataList.h ../Post/PViewData.h ../Geo/SBoundingBox3d.h \
@@ -223,8 +220,7 @@ ExtractElements${OBJEXT}: ExtractElements.cpp ExtractElements.h Plugin.h \
   ../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
   ../Post/PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
   ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/ListUtils.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
 ExtractEdges${OBJEXT}: ExtractEdges.cpp ExtractEdges.h Plugin.h \
   ../Common/Options.h ../Post/ColorTable.h ../Common/GmshMessage.h \
   ../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
@@ -260,28 +256,26 @@ Integrate${OBJEXT}: Integrate.cpp Integrate.h Plugin.h ../Common/Options.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/ListUtils.h ../Post/shapeFunctions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Post/PViewOptions.h ../Post/ColorTable.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Post/PViewOptions.h \
+  ../Post/ColorTable.h
 Gradient${OBJEXT}: Gradient.cpp Gradient.h Plugin.h ../Common/Options.h \
   ../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/ListUtils.h ../Post/shapeFunctions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
 Curl${OBJEXT}: Curl.cpp Curl.h Plugin.h ../Common/Options.h ../Post/ColorTable.h \
   ../Common/GmshMessage.h ../Post/PView.h ../Geo/SPoint3.h \
   ../Post/PViewDataList.h ../Post/PViewData.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
   ../Common/ListUtils.h ../Post/shapeFunctions.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Numeric/GmshMatrix.h
 Divergence${OBJEXT}: Divergence.cpp Divergence.h Plugin.h ../Common/Options.h \
   ../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
   ../Common/GmshConfig.h ../Common/ListUtils.h ../Post/shapeFunctions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
 Annotate${OBJEXT}: Annotate.cpp ../Common/GmshConfig.h Annotate.h Plugin.h \
   ../Common/Options.h ../Post/ColorTable.h ../Common/GmshMessage.h \
   ../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
diff --git a/Plugin/Plugin.h b/Plugin/Plugin.h
index 81ab29cea50eaf15767cd7e0aaed0d77349e37a0..fc68f4d462287da37e0e6781a5fe37e908819f9b 100644
--- a/Plugin/Plugin.h
+++ b/Plugin/Plugin.h
@@ -83,7 +83,7 @@ class GMSH_PostPlugin : public GMSH_Plugin
   // get the data in list format
   virtual PViewDataList *getDataList(PView *view);
   virtual void assignSpecificVisibility() const {}
-  virtual bool geometricalFilter(Double_Matrix *) const { return true; }
+  virtual bool geometricalFilter(gmshMatrix<double> *) const { return true; }
 };
 
 // The base class for solver plugins. The idea is to be able to
diff --git a/Post/Makefile b/Post/Makefile
index bfa6bec1693d67b3f2e26af99caafd1f155f021d..921dcfaa0865ae6b2a8ac224138d07092f8da610 100644
--- a/Post/Makefile
+++ b/Post/Makefile
@@ -51,16 +51,17 @@ depend:
 # DO NOT DELETE THIS LINE
 PView${OBJEXT}: PView.cpp PView.h ../Geo/SPoint3.h PViewDataList.h PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/ListUtils.h \
-  PViewDataGModel.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h \
-  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/GPoint.h \
+  ../Common/GmshConfig.h ../Common/ListUtils.h PViewDataGModel.h \
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/GPoint.h \
   ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
   ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
   ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
   ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h PViewOptions.h ColorTable.h \
-  ../Common/VertexArray.h ../Common/SmoothData.h adaptiveData.h
+  ../Common/VertexArray.h ../Common/SmoothData.h adaptiveData.h \
+  ../Common/GmshMessage.h
 PViewIO${OBJEXT}: PViewIO.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   PView.h ../Geo/SPoint3.h PViewDataList.h PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
@@ -74,29 +75,25 @@ PViewIO${OBJEXT}: PViewIO.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
   ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Common/StringUtils.h
 PViewData${OBJEXT}: PViewData.cpp PViewData.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
-  ../Common/GmshMessage.h ../Common/ListUtils.h adaptiveData.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h
+  ../Common/ListUtils.h adaptiveData.h ../Numeric/Numeric.h \
+  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h
 PViewDataIO${OBJEXT}: PViewDataIO.cpp ../Common/GmshMessage.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h PViewData.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
+  PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h
 PViewDataList${OBJEXT}: PViewDataList.cpp PViewDataList.h PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/ListUtils.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/SmoothData.h ../Common/Context.h \
-  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+  ../Common/GmshConfig.h ../Common/ListUtils.h ../Numeric/Numeric.h \
+  ../Numeric/GmshMatrix.h ../Common/SmoothData.h ../Common/GmshMessage.h \
+  ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
 PViewDataListIO${OBJEXT}: PViewDataListIO.cpp PViewDataList.h PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Common/ListUtils.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
-  ../Mesh/meshPartitionOptions.h adaptiveData.h
+  ../Common/GmshConfig.h ../Common/ListUtils.h ../Numeric/Numeric.h \
+  ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h adaptiveData.h
 PViewDataGModel${OBJEXT}: PViewDataGModel.cpp PViewDataGModel.h PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Common/GmshConfig.h ../Geo/GModel.h ../Geo/GVertex.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
   ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
   ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
@@ -106,8 +103,8 @@ PViewDataGModel${OBJEXT}: PViewDataGModel.cpp PViewDataGModel.h PViewData.h \
   ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
   ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
+  ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+  ../Numeric/GmshMatrix.h ../Numeric/Gauss.h ../Numeric/Numeric.h \
   ../Numeric/GmshMatrix.h
 PViewDataGModelIO${OBJEXT}: PViewDataGModelIO.cpp ../Common/GmshConfig.h \
   ../Common/GmshMessage.h PViewDataGModel.h PViewData.h \
@@ -124,38 +121,36 @@ PViewDataGModelIO${OBJEXT}: PViewDataGModelIO.cpp ../Common/GmshConfig.h \
   ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/StringUtils.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/StringUtils.h
 PViewOptions${OBJEXT}: PViewOptions.cpp ../Common/GmshConfig.h \
   ../Common/GmshMessage.h PViewOptions.h ColorTable.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h
 adaptiveData${OBJEXT}: adaptiveData.cpp adaptiveData.h ../Numeric/GmshMatrix.h \
-  ../Common/GmshConfig.h ../Common/GmshMessage.h ../Plugin/Plugin.h \
-  ../Common/Options.h ../Post/ColorTable.h ../Post/PView.h \
+  ../Common/GmshConfig.h ../Plugin/Plugin.h ../Common/Options.h \
+  ../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
   ../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/ListUtils.h \
   ../Common/OS.h
 shapeFunctions${OBJEXT}: shapeFunctions.cpp shapeFunctions.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
-  ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h
+  ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
+  ../Common/GmshMessage.h
 OctreePost${OBJEXT}: OctreePost.cpp ../Common/Octree.h \
   ../Common/OctreeInternals.h OctreePost.h ../Common/ListUtils.h PView.h \
   ../Geo/SPoint3.h PViewDataList.h PViewData.h ../Geo/SBoundingBox3d.h \
   ../Geo/SPoint3.h ../Numeric/GmshMatrix.h ../Common/GmshConfig.h \
-  ../Common/GmshMessage.h PViewDataGModel.h ../Geo/GModel.h \
-  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
-  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
-  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
-  ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h shapeFunctions.h \
-  ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
-  ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  PViewDataGModel.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h \
+  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/GPoint.h \
+  ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
+  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+  ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
+  ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+  ../Common/GmshMessage.h shapeFunctions.h ../Geo/MElement.h \
+  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+  ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h
 ColorTable${OBJEXT}: ColorTable.cpp ../Common/GmshMessage.h ColorTable.h \
   ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
-  ../Common/GmshConfig.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Numeric/GmshMatrix.h
+  ../Common/GmshConfig.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 53ef81ec40d375ff1ab043803698103e84427557..af43e33958ec369c669e9c302a66eda807d56eb5 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -17,7 +17,7 @@ PViewData::PViewData()
 PViewData::~PViewData()
 {
   if(_adaptive) delete _adaptive;
-  for(std::map<int, std::vector<Double_Matrix*> >::iterator it = _interpolation.begin();
+  for(std::map<int, std::vector<gmshMatrix<double>*> >::iterator it = _interpolation.begin();
       it != _interpolation.end(); it++)
     for(unsigned int i = 0; i < it->second.size(); i++)
       delete it->second[i];
@@ -74,28 +74,28 @@ void PViewData::setValue(int step, int ent, int ele, int nod, int comp, double v
 }
 
 void PViewData::setInterpolationMatrices(int type, 
-                                         const Double_Matrix &coefVal,
-                                         const Double_Matrix &expVal)
+                                         const gmshMatrix<double> &coefVal,
+                                         const gmshMatrix<double> &expVal)
 {
   if(!type || _interpolation[type].size()) return;
-  _interpolation[type].push_back(new Double_Matrix(coefVal));
-  _interpolation[type].push_back(new Double_Matrix(expVal));
+  _interpolation[type].push_back(new gmshMatrix<double>(coefVal));
+  _interpolation[type].push_back(new gmshMatrix<double>(expVal));
 }
 
 void PViewData::setInterpolationMatrices(int type, 
-                                         const Double_Matrix &coefVal,
-                                         const Double_Matrix &expVal, 
-                                         const Double_Matrix &coefGeo,
-                                         const Double_Matrix &expGeo)
+                                         const gmshMatrix<double> &coefVal,
+                                         const gmshMatrix<double> &expVal, 
+                                         const gmshMatrix<double> &coefGeo,
+                                         const gmshMatrix<double> &expGeo)
 {
   if(!type || _interpolation[type].size()) return;
-  _interpolation[type].push_back(new Double_Matrix(coefVal));
-  _interpolation[type].push_back(new Double_Matrix(expVal));
-  _interpolation[type].push_back(new Double_Matrix(coefGeo));
-  _interpolation[type].push_back(new Double_Matrix(expGeo));
+  _interpolation[type].push_back(new gmshMatrix<double>(coefVal));
+  _interpolation[type].push_back(new gmshMatrix<double>(expVal));
+  _interpolation[type].push_back(new gmshMatrix<double>(coefGeo));
+  _interpolation[type].push_back(new gmshMatrix<double>(expGeo));
 }
 
-int PViewData::getInterpolationMatrices(int type, std::vector<Double_Matrix*> &p)
+int PViewData::getInterpolationMatrices(int type, std::vector<gmshMatrix<double>*> &p)
 {
   if(_interpolation.count(type)){
     p = _interpolation[type];
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 66c61094fa7df77cbffd90aef3d66eac580cd9b6..7d72d8c0f803e3a6531475b978641781805ca73d 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -35,7 +35,7 @@ class PViewData {
   adaptiveData *_adaptive;
   // interpolation matrices, indexed by the number of edges per
   // element (1 for lines, 3 for triangles, etc.)
-  std::map<int, std::vector<Double_Matrix*> > _interpolation;
+  std::map<int, std::vector<gmshMatrix<double>*> > _interpolation;
 
  public:
   PViewData();
@@ -179,14 +179,14 @@ class PViewData {
   // set/get the interpolation matrices for elements with "type"
   // number of edges
   void setInterpolationMatrices(int type, 
-                                const Double_Matrix &coefVal,
-                                const Double_Matrix &expVal);
+                                const gmshMatrix<double> &coefVal,
+                                const gmshMatrix<double> &expVal);
   void setInterpolationMatrices(int type, 
-                                const Double_Matrix &coefVal,
-                                const Double_Matrix &expVal,
-                                const Double_Matrix &coefGeo, 
-                                const Double_Matrix &expGeo);
-  int getInterpolationMatrices(int type, std::vector<Double_Matrix*> &p);
+                                const gmshMatrix<double> &coefVal,
+                                const gmshMatrix<double> &expVal,
+                                const gmshMatrix<double> &coefGeo, 
+                                const gmshMatrix<double> &expGeo);
+  int getInterpolationMatrices(int type, std::vector<gmshMatrix<double>*> &p);
   inline bool haveInterpolationMatrices(){ return !_interpolation.empty(); }
 
   // smooth the data in the view (makes it C0)
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 92e5e9bc824dc8a0db0630dfcd3782a2895b8119..0907110f0e7bf5ffaa69cc6e34ba56344a9ca330 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -197,7 +197,7 @@ void PViewDataList::_stat(List_T *list, int nbcomp, int nbelm, int nbnod, int nb
   int nbval = nbcomp * nbnod;
 
   if(haveInterpolationMatrices()){
-    std::vector<Double_Matrix*> im;
+    std::vector<gmshMatrix<double>*> im;
     if(getInterpolationMatrices(nbedg, im) == 4)
       nbnod = im[2]->size1();
     nbval = nbcomp * im[0]->size1();
@@ -248,7 +248,7 @@ void PViewDataList::_setLast(int ele, int dim, int nbnod, int nbcomp, int nbedg,
                              List_T *list, int nblist)
 {
   if(haveInterpolationMatrices()){
-    std::vector<Double_Matrix*> im;
+    std::vector<gmshMatrix<double>*> im;
     if(getInterpolationMatrices(nbedg, im) == 4)
       nbnod = im[2]->size1();
   }
@@ -665,10 +665,10 @@ bool PViewDataList::combineSpace(nameData &nd)
 
     // copy interpolation from first merged dataset, if any
     if(!i){
-      for(std::map<int, std::vector<Double_Matrix*> >::iterator it = 
+      for(std::map<int, std::vector<gmshMatrix<double>*> >::iterator it = 
             l->_interpolation.begin(); it != l->_interpolation.end(); it++)
         for(unsigned int i = 0; i < it->second.size(); i++)
-          _interpolation[it->first].push_back(new Double_Matrix(*it->second[i]));
+          _interpolation[it->first].push_back(new gmshMatrix<double>(*it->second[i]));
     }
     
     // merge elememts
@@ -793,10 +793,10 @@ bool PViewDataList::combineTime(nameData &nd)
   }
   NbT2 = data[0]->NbT2;
   NbT3 = data[0]->NbT3;
-  for(std::map<int, std::vector<Double_Matrix*> >::iterator it = 
+  for(std::map<int, std::vector<gmshMatrix<double>*> >::iterator it = 
         data[0]->_interpolation.begin(); it != data[0]->_interpolation.end(); it++)
     for(unsigned int i = 0; i < it->second.size(); i++)
-      _interpolation[it->first].push_back(new Double_Matrix(*it->second[i]));
+      _interpolation[it->first].push_back(new gmshMatrix<double>(*it->second[i]));
   
   // merge values for all element types
   for(int i = 0; i < 24; i++){
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index e642d5764faacba31c8da54b984baf192461044a..d25b469c631b055169597cc04af0212564085ebb 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -50,9 +50,9 @@ static void cleanElement()
   T::allPoints.clear();
 }
 
-static void computeShapeFunctions(Double_Matrix *coeffs, Double_Matrix *eexps,
-				  double u, double v, double w, Double_Vector *sf,
-                                  Double_Vector *tmp)
+static void computeShapeFunctions(gmshMatrix<double> *coeffs, gmshMatrix<double> *eexps,
+				  double u, double v, double w, gmshVector<double> *sf,
+                                  gmshVector<double> *tmp)
 {
   for(int i = 0; i < eexps->size1(); i++) {
     (*tmp)(i) = pow(u, (*eexps)(i, 0));
@@ -830,7 +830,7 @@ void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
 }
 
 template <class T>
-adaptiveElements<T>::adaptiveElements(std::vector<Double_Matrix*> &p)
+adaptiveElements<T>::adaptiveElements(std::vector<gmshMatrix<double>*> &p)
   : _coeffsVal(0), _eexpsVal(0), _interpolVal(0),
     _coeffsGeom(0), _eexpsGeom(0), _interpolGeom(0)
 {
@@ -864,15 +864,15 @@ void adaptiveElements<T>::init(int level)
   int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
   
   if(_interpolVal) delete _interpolVal;
-  _interpolVal = new Double_Matrix(T::allPoints.size(), numVals);
+  _interpolVal = new gmshMatrix<double>(T::allPoints.size(), numVals);
   
   if(_interpolGeom) delete _interpolGeom;
-  _interpolGeom = new Double_Matrix(T::allPoints.size(), numNodes);
+  _interpolGeom = new gmshMatrix<double>(T::allPoints.size(), numNodes);
   
-  Double_Vector sfv(numVals), *tmpv = 0;
-  Double_Vector sfg(numNodes), *tmpg = 0;
-  if(_eexpsVal) tmpv = new Double_Vector(_eexpsVal->size1());
-  if(_eexpsGeom) tmpg = new Double_Vector(_eexpsGeom->size1());
+  gmshVector<double> sfv(numVals), *tmpv = 0;
+  gmshVector<double> sfg(numNodes), *tmpg = 0;
+  if(_eexpsVal) tmpv = new gmshVector<double>(_eexpsVal->size1());
+  if(_eexpsGeom) tmpg = new gmshVector<double>(_eexpsGeom->size1());
 
   int i = 0;
   for(std::set<adaptivePoint>::iterator it = T::allPoints.begin(); 
@@ -937,7 +937,7 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
   double t1 = GetTimeInSeconds();
 #endif
   
-  Double_Vector val(numVals), res(numPoints);
+  gmshVector<double> val(numVals), res(numPoints);
   if(numComp == 1){
     for(int i = 0; i < numVals; i++)
       val(i) = values[i].v[0];
@@ -957,10 +957,10 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
   }
   if(onlyComputeMinMax) return;
   
-  Double_Matrix *resxyz = 0;
+  gmshMatrix<double> *resxyz = 0;
   if(numComp == 3){
-    Double_Matrix valxyz(numVals, 3);
-    resxyz = new Double_Matrix(numPoints, 3);
+    gmshMatrix<double> valxyz(numVals, 3);
+    resxyz = new gmshMatrix<double>(numPoints, 3);
     for(int i = 0; i < numVals; i++){
       valxyz(i, 0) = values[i].v[0];
       valxyz(i, 1) = values[i].v[1];
@@ -976,7 +976,7 @@ void adaptiveElements<T>::adapt(double tol, int numComp,
     return;
   }
   
-  Double_Matrix xyz(numNodes, 3), XYZ(numPoints, 3);
+  gmshMatrix<double> xyz(numNodes, 3), XYZ(numPoints, 3);
   for(int i = 0; i < numNodes; i++){
     xyz(i, 0) = coords[i].c[0];
     xyz(i, 1) = coords[i].c[1];
@@ -1134,7 +1134,7 @@ adaptiveData::adaptiveData(PViewData *data)
     _tetrahedra(0), _hexahedra(0), _prisms(0)
 {
   _outData = new PViewDataList(true);
-  std::vector<Double_Matrix*> p;
+  std::vector<gmshMatrix<double>*> p;
   if(_inData->getNumLines()){
     _inData->getInterpolationMatrices(1, p);
     _lines = new adaptiveElements<adaptiveLine>(p);
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index 2bcd2a6604d2f73196636cadadd0274413dc8a04..415193c543491029875abd81837d81c94faa43c2 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -53,7 +53,7 @@ class adaptiveLine {
   {
     return (p[0]->val + p[1]->val) / 2.;
   }
-  inline static void GSF(double u, double v, double w, Double_Vector &sf)
+  inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
   {
     sf(0) = (1 - u) / 2.;
     sf(1) = (1 + u) / 2.;
@@ -85,7 +85,7 @@ class adaptiveTriangle {
   {
     return (p[0]->val + p[1]->val + p[2]->val) / 3.;
   }
-  inline static void GSF(double u, double v, double w, Double_Vector &sf)
+  inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
   {
     sf(0) = 1. - u - v;
     sf(1) = u;
@@ -120,7 +120,7 @@ class adaptiveQuadrangle {
   {
     return (p[0]->val + p[1]->val + p[2]->val + p[3]->val) / 4.;
   }
-  inline static void GSF(double u, double v, double w, Double_Vector &sf)
+  inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
   {
     sf(0) = 0.25 * (1. - u) * (1. - v);
     sf(1) = 0.25 * (1. + u) * (1. - v);
@@ -161,7 +161,7 @@ class adaptivePrism {
   {
     return (p[0]->val + p[1]->val + p[2]->val + p[3]->val + p[4]->val + p[5]->val) / 6.;
   }
-  inline static void GSF(double u, double v, double w, Double_Vector &sf)
+  inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
   {
     sf(0) = (1. - u - v) * (1 - w) / 2;
     sf(1) = u * (1-w)/2;
@@ -200,7 +200,7 @@ class adaptiveTetrahedron {
   {
     return (p[0]->val + p[1]->val + p[2]->val + p[3]->val) / 4.;
   }
-  inline static void GSF(double u, double v, double w, Double_Vector &sf)
+  inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
   {
     sf(0) = 1. - u - v - w;
     sf(1) = u;
@@ -243,7 +243,7 @@ class adaptiveHexahedron {
     return (p[0]->val + p[1]->val + p[2]->val+ p[3]->val +
             p[4]->val + p[5]->val + p[6]->val+ p[7]->val) / 8.;
   }
-  inline static void GSF(double u, double v, double w, Double_Vector &sf)
+  inline static void GSF(double u, double v, double w, gmshVector<double> &sf)
   {
     sf(0) = 0.125 * (1 - u) * (1 - v) * (1 - w);
     sf(1) = 0.125 * (1 + u) * (1 - v) * (1 - w);
@@ -285,10 +285,10 @@ class PValues{
 template <class T>
 class adaptiveElements {
  private:
-  Double_Matrix *_coeffsVal, *_eexpsVal, *_interpolVal;
-  Double_Matrix *_coeffsGeom, *_eexpsGeom, *_interpolGeom;
+  gmshMatrix<double> *_coeffsVal, *_eexpsVal, *_interpolVal;
+  gmshMatrix<double> *_coeffsGeom, *_eexpsGeom, *_interpolGeom;
  public:
-  adaptiveElements(std::vector<Double_Matrix*> &interpolationMatrices);
+  adaptiveElements(std::vector<gmshMatrix<double>*> &interpolationMatrices);
   ~adaptiveElements();
   // create the _interpolVal and _interpolGeom matrices at the given
   // refinement level
diff --git a/README b/README
index e02ae898d6abf94ce2f1c9f7ef6870ba3cdcb0d5..46afa6528ebf525387ec2b5ea987f829056959d2 100644
--- a/README
+++ b/README
@@ -7,12 +7,11 @@ To install Gmsh, type
 make
 make install
 
-This requires GSL 1.2 or above (freely available from
-http://sources.redhat.com/gsl/) and FLTK 1.1.7 or above (configured
-with OpenGL support; freely available from http://www.fltk.org). You
-can use the --with-fltk-prefix and --with-gsl-prefix configure options
-(or define the FLTK_PREFIX and GSL_PREFIX environment variables) if
-the libraries are not installed in their default locations.
+This FLTK 1.1.7 or above (configured with OpenGL support; freely
+available from http://www.fltk.org). You can use the
+--with-fltk-prefix configure options (or define the FLTK_PREFIX
+environment variable) if the libraries are not installed in their
+default locations.
 
 To install a non-graphical version of Gmsh (that does not require FLTK
 nor OpenGL), type
diff --git a/configure b/configure
index 12b2a7f6bbbf3bd8afb9df9f127715b057900ffd..413c6ed27e31aee0f164dc10522d139a25200d62 100755
--- a/configure
+++ b/configure
@@ -1263,7 +1263,6 @@ if test -n "$ac_init_help"; then
 Optional Features:
   --disable-FEATURE       do not include FEATURE (same as --enable-FEATURE=no)
   --enable-FEATURE[=ARG]  include FEATURE [ARG=yes]
-  --enable-gsl            use GSL as numerical toolkit (default=yes)
   --enable-gui            build the graphical user interface (default=yes)
   --enable-parser         build the parser (default=yes)
   --enable-post           build the post-processing module (default=yes)
@@ -1298,7 +1297,6 @@ Optional Packages:
   --with-PACKAGE[=ARG]    use PACKAGE [ARG=yes]
   --without-PACKAGE       do not use PACKAGE (same as --with-PACKAGE=no)
   --with-fltk-prefix=PFX  prefix where FLTK is installed
-  --with-gsl-prefix=PFX   prefix where GSL is installed
   --with-jpeg-prefix=PFX  prefix where the JPEG library and includes are
                           installed
   --with-png-prefix=PFX   prefix where the PNG library and includes are
@@ -1765,12 +1763,6 @@ if test "${with_fltk_prefix+set}" = set; then
 fi
 
 
-# Check whether --with-gsl-prefix was given.
-if test "${with_gsl_prefix+set}" = set; then
-  withval=$with_gsl_prefix; GSL_PREFIX=$withval
-fi
-
-
 # Check whether --with-jpeg-prefix was given.
 if test "${with_jpeg_prefix+set}" = set; then
   withval=$with_jpeg_prefix; JPEG_PREFIX=$withval
@@ -1849,11 +1841,6 @@ if test "${with_blas_lapack_prefix+set}" = set; then
 fi
 
 
-# Check whether --enable-gsl was given.
-if test "${enable_gsl+set}" = set; then
-  enableval=$enable_gsl;
-fi
-
 # Check whether --enable-gui was given.
 if test "${enable_gui+set}" = set; then
   enableval=$enable_gui;
@@ -1987,7 +1974,6 @@ fi
 
 if test "x$enable_minimal" = "xyes"; then
   enable_gui=no;
-  enable_gsl=no;
   enable_fm=no;
   enable_netgen=no;
   enable_tetgen=no;
@@ -5185,84 +5171,6 @@ _ACEOF
   fi
 fi
 
-if test "x$enable_gsl" != "xno"; then
-  if test "x${GSL_PREFIX}" != "x"; then
-    LDFLAGS="-L${GSL_PREFIX} -L${GSL_PREFIX}/lib ${LDFLAGS}"
-  fi
-  { echo "$as_me:$LINENO: checking for main in -lgsl" >&5
-echo $ECHO_N "checking for main in -lgsl... $ECHO_C" >&6; }
-if test "${ac_cv_lib_gsl_main+set}" = set; then
-  echo $ECHO_N "(cached) $ECHO_C" >&6
-else
-  ac_check_lib_save_LIBS=$LIBS
-LIBS="-lgsl -lgslcblas $LIBS"
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-
-
-int
-main ()
-{
-return main ();
-  ;
-  return 0;
-}
-_ACEOF
-rm -f conftest.$ac_objext conftest$ac_exeext
-if { (ac_try="$ac_link"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5
-  (eval "$ac_link") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } && {
-	 test -z "$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       } && test -s conftest$ac_exeext &&
-       $as_test_x conftest$ac_exeext; then
-  ac_cv_lib_gsl_main=yes
-else
-  echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-	ac_cv_lib_gsl_main=no
-fi
-
-rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \
-      conftest$ac_exeext conftest.$ac_ext
-LIBS=$ac_check_lib_save_LIBS
-fi
-{ echo "$as_me:$LINENO: result: $ac_cv_lib_gsl_main" >&5
-echo "${ECHO_T}$ac_cv_lib_gsl_main" >&6; }
-if test $ac_cv_lib_gsl_main = yes; then
-  GSL="yes"
-fi
-
-  if test "x${GSL}" = "xyes"; then
-    cat >>confdefs.h <<\_ACEOF
-#define HAVE_GSL 1
-_ACEOF
-
-    BO="${BO} Gsl"
-    if test "x${GSL_PREFIX}" = "x"; then
-      GMSH_LIBS="${GMSH_LIBS} -lgsl"
-    else
-      GMSH_LIBS="${GMSH_LIBS} -L${GSL_PREFIX} -L${GSL_PREFIX}/lib -lgsl"
-      FLAGS="${FLAGS} -I${GSL_PREFIX} -I${GSL_PREFIX}/include"
-    fi
-  fi
-fi
-
 if test "x$enable_fm" != "xno"; then
   if test "x${FM_PREFIX}" != "x"; then
     LDFLAGS="-L${FM_PREFIX}/lib ${LDFLAGS}"
@@ -5415,160 +5323,7 @@ _ACEOF
   fi
 fi
 
-if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
-  LDFLAGS="-L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib ${LDFLAGS}"
-fi
-{ echo "$as_me:$LINENO: checking for cblas_dgemm in -lcblas" >&5
-echo $ECHO_N "checking for cblas_dgemm in -lcblas... $ECHO_C" >&6; }
-if test "${ac_cv_lib_cblas_cblas_dgemm+set}" = set; then
-  echo $ECHO_N "(cached) $ECHO_C" >&6
-else
-  ac_check_lib_save_LIBS=$LIBS
-LIBS="-lcblas  $LIBS"
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-
-/* Override any GCC internal prototype to avoid an error.
-   Use char because int might match the return type of a GCC
-   builtin and then its argument prototype would still apply.  */
-#ifdef __cplusplus
-extern "C"
-#endif
-char cblas_dgemm ();
-int
-main ()
-{
-return cblas_dgemm ();
-  ;
-  return 0;
-}
-_ACEOF
-rm -f conftest.$ac_objext conftest$ac_exeext
-if { (ac_try="$ac_link"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5
-  (eval "$ac_link") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } && {
-	 test -z "$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       } && test -s conftest$ac_exeext &&
-       $as_test_x conftest$ac_exeext; then
-  ac_cv_lib_cblas_cblas_dgemm=yes
-else
-  echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-	ac_cv_lib_cblas_cblas_dgemm=no
-fi
-
-rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \
-      conftest$ac_exeext conftest.$ac_ext
-LIBS=$ac_check_lib_save_LIBS
-fi
-{ echo "$as_me:$LINENO: result: $ac_cv_lib_cblas_cblas_dgemm" >&5
-echo "${ECHO_T}$ac_cv_lib_cblas_cblas_dgemm" >&6; }
-if test $ac_cv_lib_cblas_cblas_dgemm = yes; then
-  CBLAS="yes" BLAS_LIBS="-lcblas"
-fi
-
-if test "x${CBLAS}" != "xyes"; then
-  { echo "$as_me:$LINENO: checking for cblas_dgemm in -lcblas" >&5
-echo $ECHO_N "checking for cblas_dgemm in -lcblas... $ECHO_C" >&6; }
-if test "${ac_cv_lib_cblas_cblas_dgemm+set}" = set; then
-  echo $ECHO_N "(cached) $ECHO_C" >&6
-else
-  ac_check_lib_save_LIBS=$LIBS
-LIBS="-lcblas -latlas $LIBS"
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-
-/* Override any GCC internal prototype to avoid an error.
-   Use char because int might match the return type of a GCC
-   builtin and then its argument prototype would still apply.  */
-#ifdef __cplusplus
-extern "C"
-#endif
-char cblas_dgemm ();
-int
-main ()
-{
-return cblas_dgemm ();
-  ;
-  return 0;
-}
-_ACEOF
-rm -f conftest.$ac_objext conftest$ac_exeext
-if { (ac_try="$ac_link"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5
-  (eval "$ac_link") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } && {
-	 test -z "$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       } && test -s conftest$ac_exeext &&
-       $as_test_x conftest$ac_exeext; then
-  ac_cv_lib_cblas_cblas_dgemm=yes
-else
-  echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-	ac_cv_lib_cblas_cblas_dgemm=no
-fi
-
-rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \
-      conftest$ac_exeext conftest.$ac_ext
-LIBS=$ac_check_lib_save_LIBS
-fi
-{ echo "$as_me:$LINENO: result: $ac_cv_lib_cblas_cblas_dgemm" >&5
-echo "${ECHO_T}$ac_cv_lib_cblas_cblas_dgemm" >&6; }
-if test $ac_cv_lib_cblas_cblas_dgemm = yes; then
-  CBLAS="yes" BLAS_LIBS="-lcblas -latlas"
-fi
-
-fi
-if test "x${CBLAS}" = "xyes"; then
-  cat >>confdefs.h <<\_ACEOF
-#define HAVE_CBLAS 1
-_ACEOF
-
-  BO="${BO} Cblas"
-else
-  if test "x${GSL}" = "xyes"; then
-        BLAS_LIBS="-lgslcblas"
-    cat >>confdefs.h <<\_ACEOF
-#define HAVE_CBLAS 1
-_ACEOF
-
-    BO="${BO} Cblas"
-  fi
-fi
-
-if test "x${FM}" = "xyes" -o "x${GSL}" != "xyes"; then
-  ac_ext=f
+ac_ext=f
 ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5'
 ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5'
 ac_compiler_gnu=$ac_cv_f77_compiler_gnu
@@ -5828,19 +5583,22 @@ ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5'
 ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5'
 ac_compiler_gnu=$ac_cv_cxx_compiler_gnu
 
-  case "${F77}" in
-    *gfortran*)
-      F77LIB="-lgfortran"
-      ;;
-    *g77*)
-      F77LIB="-lg2c"
-      ;;
-    *)
-      F77LIB=""
-      ;;
-  esac
-  LDFLAGS="${LDFLAGS} ${F77LIB}"
-  { echo "$as_me:$LINENO: checking for ATL_xerbla in -latlas" >&5
+case "${F77}" in
+  *gfortran*)
+    F77LIB="-lgfortran"
+    ;;
+  *g77*)
+    F77LIB="-lg2c"
+    ;;
+  *)
+    F77LIB=""
+    ;;
+esac
+LDFLAGS="${LDFLAGS} ${F77LIB}"
+if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
+  LDFLAGS="${LDFLAGS} -L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib"
+fi
+{ echo "$as_me:$LINENO: checking for ATL_xerbla in -latlas" >&5
 echo $ECHO_N "checking for ATL_xerbla in -latlas... $ECHO_C" >&6; }
 if test "${ac_cv_lib_atlas_ATL_xerbla+set}" = set; then
   echo $ECHO_N "(cached) $ECHO_C" >&6
@@ -5969,8 +5727,8 @@ fi
 
 fi
 
-  if test "x${BLAS}" != "xyes"; then
-    { echo "$as_me:$LINENO: checking for dgemm_ in -lblas" >&5
+if test "x${BLAS}" != "xyes"; then
+  { echo "$as_me:$LINENO: checking for dgemm_ in -lblas" >&5
 echo $ECHO_N "checking for dgemm_ in -lblas... $ECHO_C" >&6; }
 if test "${ac_cv_lib_blas_dgemm_+set}" = set; then
   echo $ECHO_N "(cached) $ECHO_C" >&6
@@ -6035,14 +5793,14 @@ if test $ac_cv_lib_blas_dgemm_ = yes; then
   BLAS="yes" BLAS_LIBS="${BLAS_LIBS} -lblas"
 fi
 
-  fi
-  if test "x${BLAS}" = "xyes"; then
-    cat >>confdefs.h <<\_ACEOF
+fi
+if test "x${BLAS}" = "xyes"; then
+  cat >>confdefs.h <<\_ACEOF
 #define HAVE_BLAS 1
 _ACEOF
 
-    BO="${BO} Blas"
-    { echo "$as_me:$LINENO: checking for dbdsqr_ in -llapack" >&5
+  BO="${BO} Blas"
+  { echo "$as_me:$LINENO: checking for dbdsqr_ in -llapack" >&5
 echo $ECHO_N "checking for dbdsqr_ in -llapack... $ECHO_C" >&6; }
 if test "${ac_cv_lib_lapack_dbdsqr_+set}" = set; then
   echo $ECHO_N "(cached) $ECHO_C" >&6
@@ -6107,16 +5865,14 @@ if test $ac_cv_lib_lapack_dbdsqr_ = yes; then
   LAPACK="yes" BLAS_LIBS="-llapack ${BLAS_LIBS}"
 fi
 
-    if test "x${LAPACK}" = "xyes"; then
-      cat >>confdefs.h <<\_ACEOF
+  if test "x${LAPACK}" = "xyes"; then
+    cat >>confdefs.h <<\_ACEOF
 #define HAVE_LAPACK 1
 _ACEOF
 
-      BO="${BO} Lapack"
-    fi
+    BO="${BO} Lapack"
   fi
 fi
-
 if test "x${BLAS_LIBS}" != "x"; then
   if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
     GMSH_LIBS="${GMSH_LIBS} -L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib ${BLAS_LIBS}"
diff --git a/configure.in b/configure.in
index b2207fe87c5bf6f01d5114cdd7854fbcd17c8156..d92e38c2c0622a5a382ba358d7842347bafb4228 100644
--- a/configure.in
+++ b/configure.in
@@ -13,10 +13,6 @@ AC_ARG_WITH(fltk-prefix,
             AC_HELP_STRING([--with-fltk-prefix=PFX],
                            [prefix where FLTK is installed]),
             [FLTK_PREFIX=$withval])
-AC_ARG_WITH(gsl-prefix,
-            AC_HELP_STRING([--with-gsl-prefix=PFX],
-                           [prefix where GSL is installed]),
-            [GSL_PREFIX=$withval])
 AC_ARG_WITH(jpeg-prefix,
             AC_HELP_STRING([--with-jpeg-prefix=PFX],
                            [prefix where the JPEG library and includes are installed]),
@@ -71,9 +67,6 @@ AC_ARG_WITH(blas-lapack-prefix,
             [BLAS_LAPACK_PREFIX=$withval])
 
 dnl Parse '--enable' command line options
-AC_ARG_ENABLE(gsl,
-              AC_HELP_STRING([--enable-gsl],
-                             [use GSL as numerical toolkit (default=yes)]))
 AC_ARG_ENABLE(gui,
               AC_HELP_STRING([--enable-gui],
                              [build the graphical user interface (default=yes)]))
@@ -156,7 +149,6 @@ AC_ARG_ENABLE(minimal,
 dnl "minimal" build shortcut
 if test "x$enable_minimal" = "xyes"; then
   enable_gui=no;
-  enable_gsl=no;
   enable_fm=no;
   enable_netgen=no;
   enable_tetgen=no;
@@ -681,24 +673,6 @@ if test "x${ZLIB}" = "xyes"; then
   fi
 fi 
 
-dnl Check for GSL
-if test "x$enable_gsl" != "xno"; then
-  if test "x${GSL_PREFIX}" != "x"; then
-    LDFLAGS="-L${GSL_PREFIX} -L${GSL_PREFIX}/lib ${LDFLAGS}"
-  fi
-  AC_CHECK_LIB(gsl,main,GSL="yes",[],-lgslcblas)
-  if test "x${GSL}" = "xyes"; then
-    AC_DEFINE(HAVE_GSL)
-    BO="${BO} Gsl"
-    if test "x${GSL_PREFIX}" = "x"; then
-      GMSH_LIBS="${GMSH_LIBS} -lgsl"
-    else
-      GMSH_LIBS="${GMSH_LIBS} -L${GSL_PREFIX} -L${GSL_PREFIX}/lib -lgsl"
-      FLAGS="${FLAGS} -I${GSL_PREFIX} -I${GSL_PREFIX}/include"
-    fi
-  fi
-fi
-
 dnl Check for FourierModel
 if test "x$enable_fm" != "xno"; then
   if test "x${FM_PREFIX}" != "x"; then
@@ -733,60 +707,39 @@ if test "x$enable_fm" != "xno"; then
   fi
 fi
 
-dnl Check for C version of BLAS
+dnl Check for blas and lapack
+AC_PROG_F77
+case "${F77}" in
+  *gfortran*)
+    F77LIB="-lgfortran"
+    ;;
+  *g77*)
+    F77LIB="-lg2c"
+    ;;
+  *)
+    F77LIB=""
+    ;;
+esac
+LDFLAGS="${LDFLAGS} ${F77LIB}"
 if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
-  LDFLAGS="-L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib ${LDFLAGS}"
+  LDFLAGS="${LDFLAGS} -L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib"
 fi
-AC_CHECK_LIB(cblas,cblas_dgemm,CBLAS="yes" BLAS_LIBS="-lcblas")
-if test "x${CBLAS}" != "xyes"; then
-  AC_CHECK_LIB(cblas,cblas_dgemm,CBLAS="yes" BLAS_LIBS="-lcblas -latlas",[],-latlas)
+AC_CHECK_LIB(atlas,ATL_xerbla,
+  AC_CHECK_LIB(f77blas,dgemm_,
+   [BLAS="yes" BLAS_LIBS="${BLAS_LIBS} -lf77blas -latlas"],[],-latlas))
+if test "x${BLAS}" != "xyes"; then
+  AC_CHECK_LIB(blas,dgemm_,[BLAS="yes" BLAS_LIBS="${BLAS_LIBS} -lblas"])
 fi
-if test "x${CBLAS}" = "xyes"; then
-  AC_DEFINE(HAVE_CBLAS)
-  BO="${BO} Cblas"
-else 
-  if test "x${GSL}" = "xyes"; then
-    dnl use unoptimized gsl version
-    BLAS_LIBS="-lgslcblas"
-    AC_DEFINE(HAVE_CBLAS)
-    BO="${BO} Cblas"
+if test "x${BLAS}" = "xyes"; then
+  AC_DEFINE(HAVE_BLAS)
+  BO="${BO} Blas"
+  AC_CHECK_LIB(lapack,dbdsqr_,
+    [LAPACK="yes" BLAS_LIBS="-llapack ${BLAS_LIBS}"],[],${BLAS_LIBS})
+  if test "x${LAPACK}" = "xyes"; then
+    AC_DEFINE(HAVE_LAPACK)
+    BO="${BO} Lapack"
   fi
 fi
-
-dnl Check for Fortran version of blas and lapack (only used when not
-dnl using GSL, or of FourierModel is linked in)
-if test "x${FM}" = "xyes" -o "x${GSL}" != "xyes"; then
-  AC_PROG_F77
-  case "${F77}" in
-    *gfortran*)
-      F77LIB="-lgfortran"
-      ;;
-    *g77*)
-      F77LIB="-lg2c"
-      ;;
-    *)
-      F77LIB=""
-      ;;
-  esac
-  LDFLAGS="${LDFLAGS} ${F77LIB}"
-  AC_CHECK_LIB(atlas,ATL_xerbla,
-    AC_CHECK_LIB(f77blas,dgemm_,
-     [BLAS="yes" BLAS_LIBS="${BLAS_LIBS} -lf77blas -latlas"],[],-latlas))
-  if test "x${BLAS}" != "xyes"; then
-    AC_CHECK_LIB(blas,dgemm_,[BLAS="yes" BLAS_LIBS="${BLAS_LIBS} -lblas"])
-  fi
-  if test "x${BLAS}" = "xyes"; then
-    AC_DEFINE(HAVE_BLAS)
-    BO="${BO} Blas"
-    AC_CHECK_LIB(lapack,dbdsqr_,
-      [LAPACK="yes" BLAS_LIBS="-llapack ${BLAS_LIBS}"],[],${BLAS_LIBS})
-    if test "x${LAPACK}" = "xyes"; then
-      AC_DEFINE(HAVE_LAPACK)
-      BO="${BO} Lapack"
-    fi
-  fi
-fi
-
 if test "x${BLAS_LIBS}" != "x"; then
   if test "x${BLAS_LAPACK_PREFIX}" != "x"; then
     GMSH_LIBS="${GMSH_LIBS} -L${BLAS_LAPACK_PREFIX} -L${BLAS_LAPACK_PREFIX}/lib ${BLAS_LIBS}"
diff --git a/doc/CREDITS.txt b/doc/CREDITS.txt
index f11d58cf3139f2d76b72c7e8c252cff0e72b55b4..c0ebb70c481b2ca1e6ad7eedf4204f7f13ae2b4f 100644
--- a/doc/CREDITS.txt
+++ b/doc/CREDITS.txt
@@ -9,20 +9,20 @@
           <jean-francois.remacle at uclouvain.be>
 
 Major code contributions to Gmsh have been provided by Nicolas Tardieu
-(help with the GSL and Netgen integration), Stephen Guzik (CGNS and
-partitioner integration), Pascale Noyret (MED integration) and
-Jonathan Lambrechts (Fields).
+(help with Netgen integration), Stephen Guzik (CGNS and partitioner
+integration), Pascale Noyret (MED integration) and Jonathan Lambrechts
+(Fields).
 
 Other code contributors include: David Colignon for new colormaps and
 lots of testing and bug reports; Patrick Dular for transfinite mesh
 bug fixes; Laurent Stainier for the eigenvalue solvers and for help
 with the MacOS port and the tensor display code; Pierre Badel for help
-with the GSL integration; Marc Ume for the original list code; Matt
-Gundry for the Plot3d mesh format; Jozef Vesely for help with the
-Tetgen integration; Koen Hillewaert for high order element mappings
-and other improvements; Jacques Lechelle for the DIFFPACK mesh format;
-Ruth Sabariego for pyramids; Gaetan Bricteux for Gauss integration
-routines.
+with root finding and minimization routines; Marc Ume for the original
+list code; Matt Gundry for the Plot3d mesh format; Jozef Vesely for
+help with the Tetgen integration; Koen Hillewaert for high order
+element mappings and other improvements; Jacques Lechelle for the
+DIFFPACK mesh format; Ruth Sabariego for pyramids; Gaetan Bricteux for
+Gauss integration routines.
 
 The AVL tree code (Common/avl.*) and the YUV image code
 (Graphics/gl2yuv.*) are copyright (C) 1988-1993, 1995 The Regents of
diff --git a/doc/FAQ.txt b/doc/FAQ.txt
index 632e7586548f1cfb62bca2b04a7aa985d761f329..1862c19e9cf0d035966815d2e26da3f6c2c08312 100644
--- a/doc/FAQ.txt
+++ b/doc/FAQ.txt
@@ -1,4 +1,4 @@
-$Id: FAQ.txt,v 1.2 2009-01-08 23:58:20 geuzaine Exp $
+$Id: FAQ.txt,v 1.3 2009-02-07 17:14:49 geuzaine Exp $
 
 This is the Gmsh FAQ
 
@@ -53,9 +53,8 @@ found at http://www.mesa3d.org.
 * 2.3 What do I need to compile Gmsh from the sources?
 
 You need a C and a C++ compiler (e.g. the GNU compilers gcc and g++)
-as well as the GSL (version 1.2 or higher; freely available from
-http://sources.redhat.com/gsl/) and FLTK (version 1.1.x, configured
-with OpenGL support; freely available from http://www.fltk.org).
+as well as FLTK (version >= 1.1.7, configured with OpenGL support;
+freely available from http://www.fltk.org).
 
 * 2.4 How do I compile Gmsh?
 
diff --git a/doc/VERSIONS.txt b/doc/VERSIONS.txt
index f5350249d12947a58cc68343c180559af987d46c..0d0df4c2574422c485274c639db9b7e1eb8eaed7 100644
--- a/doc/VERSIONS.txt
+++ b/doc/VERSIONS.txt
@@ -1,6 +1,7 @@
-$Id: VERSIONS.txt,v 1.35 2009-01-30 23:33:51 geuzaine Exp $
+$Id: VERSIONS.txt,v 1.36 2009-02-07 17:14:49 geuzaine Exp $
 
-2.3.1 (?): new per-window visibility; fixes for string options in
+(?): removed GSL dependency (Gmsh now simply requires Blas and
+Lapack); new per-window visibility; fixes for string options in
 parser.
 
 2.3.0 (Jan 23, 2009): major graphics and GUI code refactoring; new
diff --git a/doc/gmsh.html b/doc/gmsh.html
index f48bc5bea65d35279d40472874acc44a31fb8fad..c31d91d8b9a6fadf65c82a2030922888f4e89ea9 100644
--- a/doc/gmsh.html
+++ b/doc/gmsh.html
@@ -284,10 +284,9 @@ Linux) or modify the LD_LIBRARY_PATH/SHLIB_PATH/etc. environment
 variable in order for Gmsh to find the libraries.
 <p>
 <a name="build-footnote"></a><a href="#build-footmark"><sup>2</sup></a>You
-need the <a href="http://sources.redhat.com/gsl/">GSL (>= 1.2)</a>
-and <a href="http://www.fltk.org/">FLTK (>= 1.1.7)</a> libraries
-properly installed on your system in order to compile
-Gmsh. Non-graphical versions can be compiled without FLTK.
+need the <a href="http://www.fltk.org/">FLTK (>= 1.1.7)</a> libraries
+properly installed on your system in order to compile the graphical
+version of Gmsh. Non-graphical versions can be compiled without FLTK.
 
 <p>
 Back to <a href="/">geuz.org</a>
diff --git a/utils/misc/driver.cpp b/utils/misc/driver.cpp
index 084c4555c9b98e6f3c84089b7f3e29db1eb02b8a..9a14be051dcb7d94003f8c8b5a7568c610f805ab 100644
--- a/utils/misc/driver.cpp
+++ b/utils/misc/driver.cpp
@@ -3,7 +3,7 @@
 //   ./configure --disable-gui
 //   make install-lib
 //
-// Then compile this driver with "g++ driver.cpp -lGmsh -lgsl -lgslcblas"
+// Then compile this driver with "g++ driver.cpp -lGmsh -llapack -lblas"
 
 #include <stdio.h>
 #include <gmsh/Gmsh.h>
diff --git a/utils/misc/variables.msvc b/utils/misc/variables.msvc
index 37d8a5d7131ae98e4b745d5c219d581f89ba60e6..95bc48a680dfaf77584c4ad2194ed0119db8691d 100644
--- a/utils/misc/variables.msvc
+++ b/utils/misc/variables.msvc
@@ -25,7 +25,6 @@ else
 endif
 
 # Change the following to select which version to build:
-ENABLE_GSL=1
 ENABLE_GUI=0
 ENABLE_PARSER=1
 ENABLE_POSTPRO=1
@@ -35,9 +34,6 @@ ENABLE_METIS=0
 ENABLE_OCC=0
 ENABLE_MED=0
 
-# If you selected ENABLE_GSL, specify where the GSL is installed
-GSL_PREFIX="E:/src/gsl-1.8"
-
 # If you selected ENABLE_GUI, specify where FLTK is installed
 FLTK_PREFIX="E:/src/fltk-1.1.9"
 
@@ -95,10 +91,6 @@ OPTIM=/O2
 GMSH_DIRS=Common Geo Mesh Numeric contrib/ANN contrib/MathEval
 
 # Optional stuff
-ifeq (${ENABLE_GSL},1)
-  FLAGS+=/DHAVE_GSL /I${GSL_PREFIX}
-  GMSH_LIBS+=${GSL_PREFIX}/lib/gsl.lib ${GSL_PREFIX}/lib/gslcblas.lib
-endif
 ifeq (${ENABLE_PARSER},1)
   GMSH_DIRS+=Parser
 else