diff --git a/SCHNApps/CMakeLists.txt b/SCHNApps/CMakeLists.txt
index 220dd4aff456c28e18524a6f1b5952cc188214ce..203b84a8a1141541c81cefd97e9745f434ba91a8 100644
--- a/SCHNApps/CMakeLists.txt
+++ b/SCHNApps/CMakeLists.txt
@@ -17,6 +17,7 @@ find_package(GLEW REQUIRED)
 find_package(Qt4 REQUIRED)
 find_package(QGLViewer REQUIRED)
 find_package(PythonLibs 2.7 REQUIRED)
+find_package(SuiteSparse REQUIRED)
 
 SET( QT_USE_QTOPENGL TRUE )
 SET( QT_USE_QTXML TRUE )
@@ -66,6 +67,7 @@ SET (EXT_INCLUDES
 	${QT_INCLUDE_DIR}
 	${QGLVIEWER_INCLUDE_DIR}
 	${PYTHON_INCLUDE_DIRS}
+	${SUITESPARSE_INCLUDE_DIRS}
 )
 
 # define libs for external libs
@@ -81,6 +83,7 @@ SET (EXT_LIBS
 	${QT_LIBRARIES}
 	${QGLVIEWER_LIBRARIES}
 	${PYTHON_LIBRARIES}
+	${SUITESPARSE_LIBRARIES}
 )
 
 
diff --git a/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h b/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h
index cba932927d91b7d18c829dee86aead353ca50919..758c3b4146f4228191ef6816310c8e9111ebfa84 100644
--- a/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h
+++ b/SCHNApps/Plugins/surfaceDeformation/include/surfaceDeformation.h
@@ -6,12 +6,17 @@
 
 #include "ui_surfaceDeformation.h"
 
-#include "Utils/drawer.h"
+#include "Container/fakeAttribute.h"
+
+#include "OpenNL/linear_solver.h"
+#include "Algo/LinearSolving/basic.h"
 
 
 using namespace CGoGN;
 using namespace SCHNApps;
 
+// namespace CGoGN { namespace Utils { class Drawer; } }
+
 
 enum SelectionMode
 {
@@ -20,6 +25,8 @@ enum SelectionMode
 };
 
 
+typedef NoNameIOAttribute<Eigen::Matrix3f> Eigen_Matrix3f ;
+
 struct PerMapParameterSet
 {
 	PerMapParameterSet() {}
@@ -32,6 +39,19 @@ struct PerMapParameterSet
 	SelectionMode verticesSelectionMode;
 	std::vector<unsigned int> locked_vertices;
 	std::vector<unsigned int> handle_vertices;
+
+	VertexAttribute<PFP2::VEC3> positionInit;
+	VertexAttribute<PFP2::VEC3> vertexNormal;
+	EdgeAttribute<PFP2::REAL> edgeAngle;
+	EdgeAttribute<PFP2::REAL> edgeWeight;
+	VertexAttribute<PFP2::REAL> vertexArea;
+	VertexAttribute<PFP2::VEC3> diffCoord;
+	VertexAttribute<Eigen_Matrix3f> vertexRotationMatrix;
+	VertexAttribute<PFP2::VEC3> rotatedDiffCoord;
+
+	VertexAttribute<unsigned int> vIndex;
+	unsigned int nb_vertices;
+	LinearSolver<PFP2::REAL>* solver;
 };
 
 struct ParameterSet
@@ -116,11 +136,14 @@ public slots:
 	void cb_selectLockedVertices(bool b);
 	void cb_selectHandleVertices(bool b);
 
+	void matchDiffCoord(View* view, MapHandlerGen* map);
+	void asRigidAsPossible(View* view, MapHandlerGen* map);
+
 private:
 	SurfaceDeformationDockTab* m_dockTab;
 	QHash<View*, ParameterSet*> h_viewParams;
 
-	Utils::Drawer* m_drawer;
+//	Utils::Drawer* m_drawer;
 
 	bool b_refreshingUI;
 
diff --git a/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp b/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp
index cfa97fa978bfda1b4ec9dfe5a8d1edfd29be1e7b..cd8f7d9e89757238e4230eeef9d37b7bf0bf0bdd 100644
--- a/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp
+++ b/SCHNApps/Plugins/surfaceDeformation/src/surfaceDeformation.cpp
@@ -3,6 +3,11 @@
 #include "Algo/Selection/raySelector.h"
 #include "Algo/Selection/collector.h"
 
+#include "Algo/Geometry/normal.h"
+#include "Algo/Geometry/laplacian.h"
+
+//#include "Utils/drawer.h"
+
 #include <QKeyEvent>
 #include <QMouseEvent>
 
@@ -10,10 +15,7 @@
 PerMapParameterSet::PerMapParameterSet(MapHandlerGen* mh) :
 	verticesSelectionMode(LOCKED)
 {
-	GenericMap* map = mh->getGenericMap();
-
-	lockingMarker = new CellMarker<VERTEX>(*map);
-	handleMarker = new CellMarker<VERTEX>(*map);
+	PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
 
 	AttributeContainer& cont = map->getAttributeContainer<VERTEX>();
 
@@ -42,6 +44,34 @@ PerMapParameterSet::PerMapParameterSet(MapHandlerGen* mh) :
 			}
 		}
 	}
+
+	lockingMarker = new CellMarker<VERTEX>(*map);
+	handleMarker = new CellMarker<VERTEX>(*map);
+
+	positionInit = mh->addAttribute<PFP2::VEC3, VERTEX>("positionInit") ;
+	map->copyAttribute(positionInit, positionAttribute) ;
+
+	vertexNormal = mh->addAttribute<PFP2::VEC3, VERTEX>("vertexNormal");
+	edgeAngle = mh->addAttribute<PFP2::REAL, EDGE>("edgeAngle");
+	edgeWeight = mh->addAttribute<PFP2::REAL, EDGE>("edgeWeight");
+	vertexArea = mh->addAttribute<PFP2::REAL, VERTEX>("vertexArea");
+	diffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("diffCoord");
+	vertexRotationMatrix = mh->addAttribute<Eigen_Matrix3f, VERTEX>("vertexRotationMatrix") ;
+	rotatedDiffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("rotatedDiffCoord") ;
+
+	Algo::Surface::Geometry::computeNormalVertices<PFP2>(*map, positionAttribute, vertexNormal) ;
+	Algo::Surface::Geometry::computeAnglesBetweenNormalsOnEdges<PFP2>(*map, positionAttribute, edgeAngle) ;
+	Algo::Surface::Geometry::computeCotanWeightEdges<PFP2>(*map, positionAttribute, edgeWeight) ;
+	Algo::Surface::Geometry::computeVoronoiAreaVertices<PFP2>(*map, positionAttribute, vertexArea) ;
+	Algo::Surface::Geometry::computeLaplacianCotanVertices<PFP2>(*map, edgeWeight, vertexArea, positionAttribute, diffCoord) ;
+
+	for(unsigned int i = vertexRotationMatrix.begin(); i != vertexRotationMatrix.end() ; vertexRotationMatrix.next(i))
+		vertexRotationMatrix[i] = Eigen::Matrix3f::Identity();
+
+	vIndex = mh->addAttribute<unsigned int, VERTEX>("vIndex");
+	nb_vertices = static_cast<AttribMap*>(map)->computeIndexCells<VERTEX>(vIndex);
+
+	solver = new LinearSolver<PFP2::REAL>(nb_vertices);
 }
 
 PerMapParameterSet::~PerMapParameterSet()
@@ -65,25 +95,22 @@ bool SurfaceDeformationPlugin::enable()
 	connect(m_window, SIGNAL(viewAndPluginUnlinked(View*, Plugin*)), this, SLOT(viewUnlinked(View*, Plugin*)));
 	connect(m_window, SIGNAL(currentViewChanged(View*)), this, SLOT(currentViewChanged(View*)));
 
-	Utils::ShaderColorPerVertex* s = new Utils::ShaderColorPerVertex();
-	registerShader(s);
-
-	m_drawer = new Utils::Drawer();
-	registerShader(m_drawer->getShader());
+//	m_drawer = new Utils::Drawer();
+//	registerShader(m_drawer->getShader());
 
 	return true;
 }
 
 void SurfaceDeformationPlugin::disable()
 {
-	delete m_drawer;
+//	delete m_drawer;
 }
 
 void SurfaceDeformationPlugin::redraw(View* view)
 {
 	if(selecting)
 	{
-		glDisable(GL_LIGHTING) ;
+/*		glDisable(GL_LIGHTING) ;
 		m_drawer->newList(GL_COMPILE_AND_EXECUTE) ;
 		m_drawer->lineWidth(2.0f) ;
 		m_drawer->begin(GL_LINES) ;
@@ -102,7 +129,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
 		m_drawer->vertex(selectionCenter + selectionRadius * PFP2::VEC3(0,0,-1)) ;
 		m_drawer->end() ;
 		m_drawer->endList() ;
-	}
+*/	}
 
 	ParameterSet* params = h_viewParams[view];
 	MapHandlerGen* mh = params->selectedMap;
@@ -112,7 +139,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
 
 		if(!perMap->locked_vertices.empty() || !perMap->handle_vertices.empty())
 		{
-			glDisable(GL_LIGHTING) ;
+/*			glDisable(GL_LIGHTING) ;
 			m_drawer->newList(GL_COMPILE_AND_EXECUTE) ;
 			m_drawer->pointSize(4.0f) ;
 			m_drawer->begin(GL_POINTS) ;
@@ -127,7 +154,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
 				m_drawer->vertex(perMap->positionAttribute[perMap->handle_vertices[i]]) ;
 			m_drawer->end() ;
 			m_drawer->endList() ;
-		}
+*/		}
 	}
 }
 
@@ -176,35 +203,26 @@ void SurfaceDeformationPlugin::mousePress(View* view, QMouseEvent* event)
 			neigh.collectAll(d) ;
 			const std::vector<Dart>& insideV = neigh.getInsideVertices() ;
 
-			if(perMap->verticesSelectionMode == LOCKED)
+			for(unsigned int i = 0; i < insideV.size(); ++i)
 			{
-				for(unsigned int i = 0; i < insideV.size(); ++i)
+				unsigned int v = map->getEmbedding<VERTEX>(insideV[i]) ;
+				if (!perMap->lockingMarker->isMarked(v))
 				{
-					unsigned int v = map->getEmbedding<VERTEX>(insideV[i]) ;
-					if (!perMap->lockingMarker->isMarked(v))
-					{
-						perMap->locked_vertices.push_back(v) ;
-						perMap->lockingMarker->mark(v);
-					}
+					perMap->locked_vertices.push_back(v) ;
+					perMap->lockingMarker->mark(v);
 				}
-			}
-			else
-			{
-				for(unsigned int i = 0; i < insideV.size(); ++i)
+				if(perMap->verticesSelectionMode == HANDLE)
 				{
-					unsigned int v = map->getEmbedding<VERTEX>(insideV[i]) ;
 					if(!perMap->handleMarker->isMarked(v))
 					{
 						perMap->handle_vertices.push_back(v) ;
 						perMap->handleMarker->mark(v);
 					}
-					if(!perMap->lockingMarker->isMarked(v))
-					{
-						perMap->locked_vertices.push_back(v) ;
-						perMap->lockingMarker->mark(v) ;
-					}
 				}
 			}
+
+			LinearSolving::resetSolver(perMap->solver, false) ;
+
 			view->updateGL() ;
 		}
 	}
@@ -243,7 +261,8 @@ void SurfaceDeformationPlugin::mouseRelease(View* view, QMouseEvent* event)
 void SurfaceDeformationPlugin::mouseMove(View* view, QMouseEvent* event)
 {
 	ParameterSet* params = h_viewParams[view];
-	PerMapParameterSet* perMap = params->getCurrentMapParameterSet();
+	MapHandlerGen* map = params->selectedMap;
+	PerMapParameterSet* perMap = params->perMap[map->getName()];
 
 	if(dragging)
 	{
@@ -257,9 +276,9 @@ void SurfaceDeformationPlugin::mouseMove(View* view, QMouseEvent* event)
 
 		dragPrevious = q ;
 
-//		matchDiffCoord() ;
+//		matchDiffCoord(view, map);
 //		for(unsigned int i = 0; i < 2; ++i)
-//			asRigidAsPossible();
+			asRigidAsPossible(view, map);
 
 		params->selectedMap->updateVBO(perMap->positionAttribute);
 
@@ -470,6 +489,150 @@ void SurfaceDeformationPlugin::cb_selectHandleVertices(bool b)
 	}
 }
 
+void SurfaceDeformationPlugin::matchDiffCoord(View* view, MapHandlerGen* mh)
+{
+	PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
+	PerMapParameterSet* perMap = h_viewParams[view]->perMap[mh->getName()];
+
+	LinearSolving::initSolver(perMap->solver, perMap->nb_vertices, true, true) ;
+	for(int coord = 0; coord < 3; ++coord)
+	{
+		LinearSolving::setupVariables<PFP2>(*map, perMap->solver, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
+		LinearSolving::startMatrix(perMap->solver) ;
+		LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
+//		LinearSolving::addRowsRHS_Laplacian_Topo<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->diffCoord, coord);
+		LinearSolving::endMatrix(perMap->solver) ;
+		LinearSolving::solve(perMap->solver) ;
+		LinearSolving::getResult<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->positionAttribute, coord) ;
+		LinearSolving::resetSolver(perMap->solver, true) ;
+	}
+}
+
+void SurfaceDeformationPlugin::asRigidAsPossible(View* view, MapHandlerGen* mh)
+{
+	PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
+	PerMapParameterSet* perMap = h_viewParams[view]->perMap[mh->getName()];
+
+	CellMarkerNoUnmark<VERTEX> m(*map) ;
+
+	for(Dart d = map->begin(); d != map->end(); map->next(d))
+	{
+		if(!m.isMarked(d))
+		{
+			m.mark(d) ;
+
+			Eigen::Matrix3f cov = Eigen::Matrix3f::Zero() ;
+			PFP2::VEC3 p = perMap->positionAttribute[d] ;
+			PFP2::VEC3 pInit = perMap->positionInit[d] ;
+			PFP2::REAL area = perMap->vertexArea[d] ;
+			Dart it = d ;
+			do
+			{
+				Dart neigh = map->phi1(it) ;
+				PFP2::VEC3 v = perMap->positionAttribute[neigh] - p ;
+				PFP2::VEC3 vv = perMap->positionInit[neigh] - pInit ;
+				for(unsigned int i = 0; i < 3; ++i)
+					for(unsigned int j = 0; j < 3; ++j)
+						cov(i,j) += v[i] * vv[j] * perMap->edgeWeight[it] / area ;
+				Dart dboundary = map->phi_1(it) ;
+				if(map->phi2(dboundary) == dboundary)
+				{
+					v = perMap->positionAttribute[dboundary] - p ;
+					vv = perMap->positionInit[dboundary] - p ;
+					for(unsigned int i = 0; i < 3; ++i)
+						for(unsigned int j = 0; j < 3; ++j)
+							cov(i,j) += v[i] * vv[j] * perMap->edgeWeight[dboundary] / area ;
+				}
+				it = map->alpha1(it) ;
+			} while(it != d) ;
+
+			Eigen::JacobiSVD<Eigen::Matrix3f> svd(cov, Eigen::ComputeFullU | Eigen::ComputeFullV) ;
+			Eigen::Matrix3f R = svd.matrixU() * svd.matrixV().transpose() ;
+
+			if(R.determinant() < 0)
+			{
+				Eigen::Matrix3f U = svd.matrixU() ;
+				for(unsigned int i = 0; i < 3; ++i)
+					U(i,2) *= -1 ;
+				R = U * svd.matrixV().transpose() ;
+			}
+
+			perMap->vertexRotationMatrix[d] = R ;
+		}
+	}
+
+	for(Dart d = map->begin(); d != map->end(); map->next(d))
+	{
+		if(m.isMarked(d))
+		{
+			m.unmark(d) ;
+
+			unsigned int degree = 0 ;
+			Eigen::Matrix3f r = Eigen::Matrix3f::Zero() ;
+			Dart it = d ;
+			do
+			{
+				r += perMap->vertexRotationMatrix[map->phi1(it)] ;
+				++degree ;
+				Dart dboundary = map->phi_1(it) ;
+				if(map->phi2(dboundary) == dboundary)
+				{
+					r += perMap->vertexRotationMatrix[dboundary] ;
+					++degree ;
+				}
+				it = map->alpha1(it) ;
+			} while(it != d) ;
+			r += perMap->vertexRotationMatrix[d] ;
+			r /= degree + 1 ;
+			PFP2::VEC3& dc = perMap->diffCoord[d] ;
+			Eigen::Vector3f rdc(dc[0], dc[1], dc[2]) ;
+			rdc = r * rdc ;
+			perMap->rotatedDiffCoord[d] = PFP2::VEC3(rdc[0], rdc[1], rdc[2]) ;
+
+//			Eigen::Vector3f rdc = Eigen::Vector3f::Zero() ;
+//			Dart it = d ;
+//			PFP::REAL vArea = perMap->vertexArea[d] ;
+//			Eigen::Matrix3f& rotM = perMap->vertexRotationMatrix[d] ;
+//			PFP::REAL val = 0 ;
+//			do
+//			{
+//				Dart ddn = map->phi1(it) ;
+//				PFP::REAL w = perMap->edgeWeight[it] / vArea ;
+//				PFP::VEC3 vv = (perMap->positionInit[ddn] - perMap->positionInit[it]) * w ;
+//				val += w ;
+//				Eigen::Matrix3f r = 0.5f * (perMap->vertexRotationMatrix[ddn] + rotM) ;
+//				Eigen::Vector3f vvr = r * Eigen::Vector3f(vv[0], vv[1], vv[2]) ;
+//				rdc += vvr ;
+//				Dart dboundary = map->phi_1(it) ;
+//				if(map->phi2(dboundary) == dboundary)
+//				{
+//					w = perMap->edgeWeight[dboundary] / vArea ;
+//					vv = (perMap->positionInit[dboundary] - perMap->positionInit[it]) * w ;
+//					val += w ;
+//					r = 0.5f * (perMap->vertexRotationMatrix[dboundary] + rotM) ;
+//					vvr = r * Eigen::Vector3f(vv[0], vv[1], vv[2]) ;
+//					rdc += vvr ;
+//				}
+//				it = map->alpha1(it) ;
+//			} while(it != d) ;
+//			rdc /= val ;
+//			perMap->rotatedDiffCoord[d] = PFP::VEC3(rdc[0], rdc[1], rdc[2]) ;
+		}
+	}
+
+	LinearSolving::initSolver(perMap->solver, perMap->nb_vertices, true, true) ;
+	for(int coord = 0; coord < 3; ++coord)
+	{
+		LinearSolving::setupVariables<PFP2>(*map, perMap->solver, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
+		LinearSolving::startMatrix(perMap->solver) ;
+		LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->rotatedDiffCoord, coord) ;
+		LinearSolving::endMatrix(perMap->solver) ;
+		LinearSolving::solve(perMap->solver) ;
+		LinearSolving::getResult<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->positionAttribute, coord) ;
+		LinearSolving::resetSolver(perMap->solver, true) ;
+	}
+}
+
 
 
 void SurfaceDeformationDockTab::refreshUI(ParameterSet* params)
diff --git a/SCHNApps/bin/init.py b/SCHNApps/bin/init.py
index a379176f25fe8b5c140e6d8905234f1b4cabd543..ea11acd184bb2b4a31febcf65ea0aebe5e8108b2 100644
--- a/SCHNApps/bin/init.py
+++ b/SCHNApps/bin/init.py
@@ -5,15 +5,15 @@ differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
 subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
 surfaceDeformationPlugin = schnapps.loadPlugin("SurfaceDeformation");
 
-obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/midRes/egea_remesh_9k.off");
+#obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_good_9k.off");
 
-v = schnapps.getView("view_0");
+#v = schnapps.getView("view_0");
 
-schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
-schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
-schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
+#schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
+#schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
+#schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
 
-schnapps.linkViewAndMap(v.getName(), obj.getName());
+#schnapps.linkViewAndMap(v.getName(), obj.getName());
 
 #differentialPropertiesPlugin.computeNormal(obj.getName());
 #differentialPropertiesPlugin.computeCurvature(obj.getName());
diff --git a/SCHNApps/src/main.cpp b/SCHNApps/src/main.cpp
index 18367da8573ced619aac58778820fa2894240a75..b2bb6472c34a1544da1cce78e9caee22e6f775df 100644
--- a/SCHNApps/src/main.cpp
+++ b/SCHNApps/src/main.cpp
@@ -1,11 +1,9 @@
-
 #include <QSplashScreen>
 #include "window.h"
 #include <QFileInfo>
 #include "PythonQt/PythonQt.h"
 #include "PythonQt/gui/PythonQtScriptingConsole.h"
 
-
 int main(int argc, char* argv[])
 {
 	QApplication app(argc, argv);
diff --git a/ThirdParty/include/OpenNL/cholesky_solver.h b/ThirdParty/include/OpenNL/cholesky_solver.h
index 3856caaa1534b2be1b26d87606d8b9fde58f1be5..8c6219347af01708d83bbf92de817bf8e2593e33 100644
--- a/ThirdParty/include/OpenNL/cholesky_solver.h
+++ b/ThirdParty/include/OpenNL/cholesky_solver.h
@@ -32,17 +32,12 @@
 #define __CHOLESKY_SOLVER__
 
 #include <cassert>
-#include "OpenNL/sparse_matrix.h"
-#include "CHOLMOD/cholmod.h"
+#include <cholmod.h>
 
-template< class MATRIX, class VECTOR >
+template <typename CoeffType>
 class Solver_CHOLESKY
 {
 public:
-    typedef MATRIX Matrix ;
-    typedef VECTOR Vector ;
-    typedef typename Vector::CoeffType CoeffType ;
-
     Solver_CHOLESKY() {
     	N = 0 ;
     	factorized_ = false ;
@@ -58,14 +53,14 @@ public:
 
     bool factorized() { return factorized_ ; }
 
-    void factorize(const MATRIX &A_in) {
-    	assert(A_in.n() == A_in.m()) ;
-		assert(A_in.has_symmetric_storage()) ;
+	void factorize(const Eigen::SparseMatrix<CoeffType>& A_in) {
+		assert(A_in.rows() == A_in.cols()) ;
+//		assert(A_in.has_symmetric_storage()) ;
 
-		N = A_in.n() ;
-		int NNZ = A_in.nnz() ;
+		N = A_in.rows() ;
 
 		// Translate sparse matrix into cholmod format
+		int NNZ = A_in.nonZeros() ;
 		cholmod_sparse* A = cholmod_allocate_sparse(N, N, NNZ, false, true, -1, CHOLMOD_REAL, &c);
 
 		int* colptr = static_cast<int*>(A->p) ;
@@ -75,12 +70,12 @@ public:
 		// Convert local Matrix into CHOLMOD Matrix
 		int count = 0 ;
 		for(int j = 0; j < N; j++) {
-			const typename SparseMatrix<CoeffType>::Column& Cj = A_in.column(j) ;
-			colptr[j] = count ;
-			for(unsigned int ii = 0; ii < Cj.nb_coeffs(); ii++) {
-				a[count] = Cj.coeff(ii).a ;
-				rowind[count] = Cj.coeff(ii).index ;
-				count++ ;
+			colptr[j] = count;
+			for(typename Eigen::SparseMatrix<CoeffType>::InnerIterator it(A_in, j); it; ++it)
+			{
+				a[count] = it.value();
+				rowind[count] = it.row();
+				++count;
 			}
 		}
 		colptr[N] = NNZ ;
@@ -94,10 +89,10 @@ public:
 		cholmod_free_sparse(&A, &c) ;
     }
 
-    void solve(const VECTOR& b_in, VECTOR& x_out) {
+	void solve(const Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>& b_in, Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>& x_out) {
     	assert(factorized_) ;
-        assert(L->n == b_in.dimension()) ;
-        assert(L->n == x_out.dimension()) ;
+		assert(L->n == b_in.rows()) ;
+		assert(L->n == x_out.rows()) ;
 
 		// Translate right-hand side into cholmod format
 		cholmod_dense* b = cholmod_allocate_dense(N, 1, N, CHOLMOD_REAL, &c) ;
@@ -114,8 +109,8 @@ public:
 			x_out[i] = cx[i] ;
 
 		// Cleanup
-		cholmod_free_dense(&x, &c) ;
 		cholmod_free_dense(&b, &c) ;
+		cholmod_free_dense(&x, &c) ;
     }
 
     void reset() {
diff --git a/ThirdParty/include/OpenNL/linear_solver.h b/ThirdParty/include/OpenNL/linear_solver.h
index 972dd2c729f6629f9122545a102405261b76e9f5..1f719ab5d5ba13317cd3fab2d103fdb2c5b7e7eb 100644
--- a/ThirdParty/include/OpenNL/linear_solver.h
+++ b/ThirdParty/include/OpenNL/linear_solver.h
@@ -37,6 +37,7 @@
 #include <cstdlib>
 
 #include <Eigen/Sparse>
+#include "cholesky_solver.h"
 
 
 template <class CoeffType>
@@ -80,7 +81,7 @@ public:
 		delete A_ ;
 		delete x_ ;
 		delete b_ ;
-//		delete direct_solver_ ;
+		delete direct_solver_ ;
 	}
 
 	// __________________ Parameters ________________________
@@ -188,9 +189,7 @@ private:
 	Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* x_ ;
 	Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* b_ ;
 
-	Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> >* symmetric_solver_;
-	Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> >* nonsymmetric_solver_;
-	Eigen::SimplicialLDLT<Eigen::SparseMatrix<CoeffType> >* direct_solver_;
+	Solver_CHOLESKY<CoeffType>* direct_solver_ ;
 } ;
 
 #include "OpenNL/linear_solver.hpp"
diff --git a/ThirdParty/include/OpenNL/linear_solver.hpp b/ThirdParty/include/OpenNL/linear_solver.hpp
index bdc792b7accef3ced7acd9003cbed8eb78b4187e..e56f766feb3bbf11acf7448013a93f4cd8a15fb6 100644
--- a/ThirdParty/include/OpenNL/linear_solver.hpp
+++ b/ThirdParty/include/OpenNL/linear_solver.hpp
@@ -33,9 +33,7 @@ LinearSolver<CoeffType>::LinearSolver(unsigned int nb_variables) {
 	A_ = NULL ;
 	x_ = NULL ;
 	b_ = NULL ;
-	symmetric_solver_ = NULL;
-	nonsymmetric_solver_ = NULL;
-	direct_solver_ = NULL;
+	direct_solver_ = new Solver_CHOLESKY<CoeffType>() ;
 }
 
 template <typename CoeffType>
@@ -138,7 +136,7 @@ void LinearSolver<CoeffType>::end_row() {
 		if(!matrix_already_set_) {
 			for(unsigned int i = 0; i < nf; i++) {
 				for(unsigned int j = 0; j < nf; j++) {
-					A_->coeffRef(if_[i], if_[j]) += af_[i] * af_[j];
+					A_->coeffRef(if_[i], if_[j]) += af_[i] * af_[j] ;
 				}
 			}
 		}
@@ -156,11 +154,11 @@ void LinearSolver<CoeffType>::end_row() {
 		// Construct the matrix coefficients of the current row
 		if(!matrix_already_set_) {
 			for(unsigned int i = 0; i < nf; i++) {
-				A_->coeffRef(current_row_, if_[i]) += af_[i];
+				A_->coeffRef(current_row_, if_[i]) += af_[i] ;
 			}
 		}
 		// Construct the right-hand side of the current row
-		(*b_)[current_row_] = bk_ ;
+		(*b_)[current_row_] = -bk_ ;
 		for(unsigned int i = 0; i < nl; i++) {
 			(*b_)[current_row_] -= al_[i] * xl_[i] ;
 		}
@@ -171,10 +169,10 @@ void LinearSolver<CoeffType>::end_row() {
 
 template <typename CoeffType>
 void LinearSolver<CoeffType>::end_system() {
-	if(least_squares_ && direct_ && !direct_solver_)
-		direct_solver_ = new Eigen::SimplicialLDLT<Eigen::SparseMatrix<CoeffType> >(*A_);
+	if(least_squares_ && direct_ && !direct_solver_->factorized())
+		direct_solver_->factorize(*A_) ;
 
-    transition(IN_SYSTEM, CONSTRUCTED) ;
+	transition(IN_SYSTEM, CONSTRUCTED) ;
 }
 
 template <typename CoeffType>
@@ -183,44 +181,48 @@ void LinearSolver<CoeffType>::solve() {
 	if(least_squares_)
 	{
 		if(direct_) {
-			*x_ = direct_solver_->solve(*b_) ;
+			direct_solver_->solve(*b_, *x_) ;
 		} else {
-			symmetric_solver_ = new Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> >(*A_) ;
-			*x_ = symmetric_solver_->solve(*b_) ;
+			Eigen::ConjugateGradient<Eigen::MatrixXf > solver(*A_) ;
+//			int n = x_->rows() ;
+//			Eigen::Matrix<CoeffType, Eigen::Dynamic, 1> guess(n) ;
+//			for(int i = 0; i < n; ++i)
+//				guess[i] = (*x_)[i] ;
+//			*x_ = solver.solveWithGuess(*b_, guess) ;
+			*x_ = solver.solve(*b_) ;
 		}
 	} else {
-		nonsymmetric_solver_ = new Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> >(*A_) ;
-		*x_ = nonsymmetric_solver_->solve(*b_) ;
+		Eigen::BiCGSTAB<Eigen::MatrixXf > solver(*A_) ;
+		*x_ = solver.solve(*b_) ;
 	}
+
 	vector_to_variables() ;
 	transition(CONSTRUCTED, SOLVED) ;
 }
 
 template <typename CoeffType>
 void LinearSolver<CoeffType>::reset(bool keep_matrix) {
-    if(keep_matrix) {
-    	x_->setZero();
-    	b_->setZero();
-    	matrix_already_set_ = true ;
-    } else {
-    	delete A_ ; A_ = NULL ;
-    	delete x_ ; x_ = NULL ;
-    	delete b_ ; b_ = NULL ;
-		if(symmetric_solver_) delete symmetric_solver_ ; symmetric_solver_ = NULL ;
-		if(nonsymmetric_solver_) delete nonsymmetric_solver_ ; nonsymmetric_solver_ = NULL ;
-		if(direct_solver_) delete direct_solver_ ; direct_solver_ = NULL ;
-    	matrix_already_set_ = false ;
-    	for(unsigned int i = 0; i < nb_variables_; ++i) {
-    		variable_[i].unlock() ;
-    	}
-    }
+	if(keep_matrix) {
+		x_->setZero() ;
+		b_->setZero() ;
+		matrix_already_set_ = true ;
+	} else {
+		delete A_ ; A_ = NULL ;
+		delete x_ ; x_ = NULL ;
+		delete b_ ; b_ = NULL ;
+		direct_solver_->reset();
+		matrix_already_set_ = false ;
+		for(unsigned int i = 0; i < nb_variables_; ++i) {
+			variable_[i].unlock() ;
+		}
+	}
 	state_ = INITIAL ;
 }
 
 template <typename CoeffType>
 void LinearSolver<CoeffType>::vector_to_variables() {
 	for(unsigned int i=0; i < nb_variables(); i++) {
-	Variable<CoeffType>& v = variable(i) ;
+		Variable<CoeffType>& v = variable(i) ;
 		if(!v.is_locked()) {
 			v.set_value((*x_)[v.index()]) ;
 		}
@@ -230,7 +232,7 @@ void LinearSolver<CoeffType>::vector_to_variables() {
 template <typename CoeffType>
 void LinearSolver<CoeffType>::variables_to_vector() {
 	for(unsigned int i=0; i < nb_variables(); i++) {
-	Variable<CoeffType>& v = variable(i) ;
+		Variable<CoeffType>& v = variable(i) ;
 		if(!v.is_locked()) {
 			(*x_)[v.index()] = v.value() ;
 		}
diff --git a/cmake_modules/FindSuiteSparse.cmake b/cmake_modules/FindSuiteSparse.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..75a8579fe97f8d953fb5b9e688aa66dfe2b9003d
--- /dev/null
+++ b/cmake_modules/FindSuiteSparse.cmake
@@ -0,0 +1,126 @@
+# - Try to find SUITESPARSE
+# Once done this will define
+#  
+#  SUITESPARSE_FOUND            - system has SUITESPARSE
+#  SUITESPARSE_INCLUDE_DIRS     - the SUITESPARSE include directory
+#  SUITESPARSE_LIBRARIES        - Link these to use SUITESPARSE
+#  SUITESPARSE_SPQR_LIBRARY     - name of spqr library (necessary due to error in debian package)
+#  SUITESPARSE_SPQR_LIBRARY_DIR - name of spqr library (necessary due to error in debian package)
+#  SUITESPARSE_LIBRARY_DIR      - Library main directory containing suitesparse libs
+#  SUITESPARSE_LIBRARY_DIRS     - all Library directories containing suitesparse libs
+#  SUITESPARSE_SPQR_VALID       - automatic identification whether or not spqr package is installed correctly
+
+IF (SUITESPARSE_INCLUDE_DIRS)
+  # Already in cache, be silent
+  SET(SUITESPARSE_FIND_QUIETLY TRUE)
+ENDIF (SUITESPARSE_INCLUDE_DIRS)
+
+if( WIN32 )
+   # Find cholmod part of the suitesparse library collection
+
+   FIND_PATH( CHOLMOD_INCLUDE_DIR cholmod.h
+              PATHS "C:\\libs\\win32\\SuiteSparse\\Include"  )
+
+   # Add cholmod include directory to collection include directories
+   IF ( CHOLMOD_INCLUDE_DIR )
+	list ( APPEND SUITESPARSE_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIR} )
+   ENDIF( CHOLMOD_INCLUDE_DIR )
+
+
+   # find path suitesparse library
+   FIND_PATH( SUITESPARSE_LIBRARY_DIRS 
+	         amd.lib
+               PATHS "C:\\libs\\win32\\SuiteSparse\\libs" )
+
+   # if we found the library, add it to the defined libraries
+   IF ( SUITESPARSE_LIBRARY_DIRS )
+	list ( APPEND SUITESPARSE_LIBRARIES optimized;amd;optimized;camd;optimized;ccolamd;optimized;cholmod;optimized;colamd;optimized;metis;optimized;spqr;optimized;umfpack;debug;amdd;debug;camdd;debug;ccolamdd;debug;cholmodd;debug;spqrd;debug;umfpackd;debug;colamdd;debug;metisd;optimized;blas;optimized;libf2c;optimized;lapack;debug;blasd;debug;libf2cd;debug;lapackd )
+   ENDIF( SUITESPARSE_LIBRARY_DIRS )  
+
+else( WIN32 )
+   IF(APPLE)
+	   FIND_PATH( CHOLMOD_INCLUDE_DIR cholmod.h
+        	      PATHS /opt/local/include/ufsparse
+			    /usr/local/include )
+
+           FIND_PATH( SUITESPARSE_LIBRARY_DIR
+                      NAMES libcholmod.a
+                      PATHS /opt/local/lib
+			    /usr/local/lib )
+   ELSE(APPLE)
+	   FIND_PATH( CHOLMOD_INCLUDE_DIR cholmod.h
+        	      PATHS /usr/local/include 
+        	            /usr/include 
+        	            /usr/include/suitesparse/ 
+        	            ${CMAKE_SOURCE_DIR}/MacOS/Libs/cholmod
+              	      PATH_SUFFIXES cholmod/ CHOLMOD/ )
+   	
+           FIND_PATH( SUITESPARSE_LIBRARY_DIR
+                      NAMES libcholmod.so 
+                      PATHS /usr/lib 
+                            /usr/lib64 
+                            /usr/local/lib )
+   ENDIF(APPLE)
+
+   # Add cholmod include directory to collection include directories
+   IF ( CHOLMOD_INCLUDE_DIR )
+	list ( APPEND SUITESPARSE_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIR} )
+   ENDIF( CHOLMOD_INCLUDE_DIR )
+
+   # if we found the library, add it to the defined libraries
+   IF ( SUITESPARSE_LIBRARY_DIR )
+
+       list ( APPEND SUITESPARSE_LIBRARIES amd)
+       list ( APPEND SUITESPARSE_LIBRARIES btf)
+       list ( APPEND SUITESPARSE_LIBRARIES camd)
+       list ( APPEND SUITESPARSE_LIBRARIES ccolamd)
+       list ( APPEND SUITESPARSE_LIBRARIES cholmod)
+       list ( APPEND SUITESPARSE_LIBRARIES colamd)
+ #     list ( APPEND SUITESPARSE_LIBRARIES csparse)
+       list ( APPEND SUITESPARSE_LIBRARIES cxsparse)
+       list ( APPEND SUITESPARSE_LIBRARIES klu)
+ #     list ( APPEND SUITESPARSE_LIBRARIES spqr)
+       list ( APPEND SUITESPARSE_LIBRARIES umfpack)
+
+       IF (APPLE)
+           list ( APPEND SUITESPARSE_LIBRARIES suitesparseconfig)
+       ENDIF (APPLE)
+
+       # Metis and spqr are optional
+       FIND_LIBRARY( SUITESPARSE_METIS_LIBRARY
+                     NAMES metis
+                     PATHS ${SUITESPARSE_LIBRARY_DIR} )
+       IF (SUITESPARSE_METIS_LIBRARY)			
+	  list ( APPEND SUITESPARSE_LIBRARIES metis)
+       ENDIF(SUITESPARSE_METIS_LIBRARY)
+
+       if(EXISTS  "${CHOLMOD_INCLUDE_DIR}/SuiteSparseQR.hpp")
+	  SET(SUITESPARSE_SPQR_VALID TRUE CACHE BOOL "SuiteSparseSPQR valid")
+       else()
+	  SET(SUITESPARSE_SPQR_VALID false CACHE BOOL "SuiteSparseSPQR valid")
+       endif()
+
+       if(SUITESPARSE_SPQR_VALID)
+	  FIND_LIBRARY( SUITESPARSE_SPQR_LIBRARY
+		      NAMES spqr
+		      PATHS ${SUITESPARSE_LIBRARY_DIR} )
+	  IF (SUITESPARSE_SPQR_LIBRARY)			
+	    list ( APPEND SUITESPARSE_LIBRARIES spqr)
+	  ENDIF (SUITESPARSE_SPQR_LIBRARY)
+       endif()
+       
+    ENDIF( SUITESPARSE_LIBRARY_DIR )  
+   
+endif( WIN32 )
+
+
+IF (SUITESPARSE_INCLUDE_DIRS AND SUITESPARSE_LIBRARIES)
+   IF(WIN32)
+    list (APPEND SUITESPARSE_INCLUDE_DIRS ${CHOLMOD_INCLUDE_DIR}/../../UFconfig )
+   ENDIF(WIN32)
+   SET(SUITESPARSE_FOUND TRUE)
+   MESSAGE(STATUS "Found SuiteSparse")
+ELSE (SUITESPARSE_INCLUDE_DIRS AND SUITESPARSE_LIBRARIES)
+   SET( SUITESPARSE_FOUND FALSE )
+ENDIF (SUITESPARSE_INCLUDE_DIRS AND SUITESPARSE_LIBRARIES)
+