Eigen: fold remaining OpenNL code into intern/eigen.

Differential Revision: https://developer.blender.org/D1662

19 files changed, 248 insertions, 508 deletions

diff --git a/source/blender/blenkernel/intern/softbody.c b/source/blender/blenkernel/intern/softbody.c
index 86ffb541bf8..1c3ff3d3e02 100644
--- a/source/blender/blenkernel/intern/softbody.c
+++ b/source/blender/blenkernel/intern/softbody.c
@@ -81,7 +81,6 @@ variables on the UI for now
 #include "BKE_scene.h"
 
 #include  "PIL_time.h"
-// #include  "ONL_opennl.h" remove linking to ONL for now
 
 /* callbacks for errors and interrupts and some goo */
 static int (*SB_localInterruptCallBack)(void) = NULL;
@@ -1811,14 +1810,14 @@ static void dfdx_spring(int ia, int ic, int op, float dir[3], float L, float len
 			for (j=0;j<3;j++) {
 				delta_ij = (i==j ? (1.0f): (0.0f));
 				m=factor*(dir[i]*dir[j] + (1-L/len)*(delta_ij - dir[i]*dir[j]));
-				nlMatrixAdd(ia+i, op+ic+j, m);
+				EIG_linear_solver_matrix_add(ia+i, op+ic+j, m);
 			}
 	}
 	else {
 		for (i=0;i<3;i++)
 			for (j=0;j<3;j++) {
 				m=factor*dir[i]*dir[j];
-				nlMatrixAdd(ia+i, op+ic+j, m);
+				EIG_linear_solver_matrix_add(ia+i, op+ic+j, m);
 			}
 	}
 }
@@ -1827,13 +1826,13 @@ static void dfdx_spring(int ia, int ic, int op, float dir[3], float L, float len
 static void dfdx_goal(int ia, int ic, int op, float factor)
 {
 	int i;
-	for (i=0;i<3;i++) nlMatrixAdd(ia+i, op+ic+i, factor);
+	for (i=0;i<3;i++) EIG_linear_solver_matrix_add(ia+i, op+ic+i, factor);
 }
 
 static void dfdv_goal(int ia, int ic, float factor)
 {
 	int i;
-	for (i=0;i<3;i++) nlMatrixAdd(ia+i, ic+i, factor);
+	for (i=0;i<3;i++) EIG_linear_solver_matrix_add(ia+i, ic+i, factor);
 }
 */
 static void sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float UNUSED(forcetime))
diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c
index c27c6bea160..641e50c9bde 100644
--- a/source/blender/blenlib/intern/math_solvers.c
+++ b/source/blender/blenlib/intern/math_solvers.c
@@ -57,7 +57,7 @@ bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3
 	}
 #endif
 
-	return EG3_self_adjoint_eigen_solve(3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors);
+	return EIG_self_adjoint_eigen_solve(3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors);
 }
 
 /**
@@ -70,5 +70,5 @@ bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3
  */
 void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[3], float r_V[3][3])
 {
-	EG3_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V);
+	EIG_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V);
 }
diff --git a/source/blender/bmesh/CMakeLists.txt b/source/blender/bmesh/CMakeLists.txt
index 686b8507291..d13d3e1cd0f 100644
--- a/source/blender/bmesh/CMakeLists.txt
+++ b/source/blender/bmesh/CMakeLists.txt
@@ -30,6 +30,7 @@ set(INC
 	../blentranslation
 	../makesdna
 	../../../intern/guardedalloc
+	../../../intern/eigen
 	../../../extern/rangetree
 )
 
@@ -176,13 +177,6 @@ if(WITH_INTERNATIONAL)
 	add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-	add_definitions(-DWITH_OPENNL)
-	list(APPEND INC_SYS
-		../../../intern/opennl/extern
-	)
-endif()
-
 if(WITH_FREESTYLE)
 	add_definitions(-DWITH_FREESTYLE)
 endif()
diff --git a/source/blender/bmesh/SConscript b/source/blender/bmesh/SConscript
index c53974be1a7..6bf4fdf7c78 100644
--- a/source/blender/bmesh/SConscript
+++ b/source/blender/bmesh/SConscript
@@ -40,9 +40,9 @@ incs = [
     '../makesdna',
     '../blenkernel',
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     '#/extern/bullet2/src',
     '#/extern/rangetree',
-    '#/intern/opennl/extern'
     ]
 
 defs = []
diff --git a/source/blender/bmesh/operators/bmo_smooth_laplacian.c b/source/blender/bmesh/operators/bmo_smooth_laplacian.c
index 624fbd93818..95ff3eb4858 100644
--- a/source/blender/bmesh/operators/bmo_smooth_laplacian.c
+++ b/source/blender/bmesh/operators/bmo_smooth_laplacian.c
@@ -28,18 +28,15 @@
 
 #include "MEM_guardedalloc.h"
 
-
 #include "BLI_math.h"
 
+#include "eigen_capi.h"
+
 
 #include "bmesh.h"
 
 #include "intern/bmesh_operators_private.h" /* own include */
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
-
 // #define SMOOTH_LAPLACIAN_AREA_FACTOR 4.0f  /* UNUSED */
 // #define SMOOTH_LAPLACIAN_EDGE_FACTOR 2.0f  /* UNUSED */
 #define SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE 1.8f
@@ -59,7 +56,7 @@ struct BLaplacianSystem {
 	/* Pointers to data*/
 	BMesh *bm;
 	BMOperator *op;
-	NLContext *context;
+	LinearSolver *context;
 
 	/*Data*/
 	float min_area;
@@ -92,7 +89,7 @@ static void delete_laplacian_system(LaplacianSystem *sys)
 	delete_void_pointer(sys->vweights);
 	delete_void_pointer(sys->zerola);
 	if (sys->context) {
-		nlDeleteContext(sys->context);
+		EIG_linear_solver_delete(sys->context);
 	}
 	sys->bm = NULL;
 	sys->op = NULL;
@@ -333,9 +330,9 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 					w4 = w4 / 4.0f;
 
 					if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) {
-						nlMatrixAdd(sys->context, idv1, idv2, w2 * sys->vweights[idv1]);
-						nlMatrixAdd(sys->context, idv1, idv3, w3 * sys->vweights[idv1]);
-						nlMatrixAdd(sys->context, idv1, idv4, w4 * sys->vweights[idv1]);
+						EIG_linear_solver_matrix_add(sys->context, idv1, idv2, w2 * sys->vweights[idv1]);
+						EIG_linear_solver_matrix_add(sys->context, idv1, idv3, w3 * sys->vweights[idv1]);
+						EIG_linear_solver_matrix_add(sys->context, idv1, idv4, w4 * sys->vweights[idv1]);
 					}
 				}
 			}
@@ -346,16 +343,16 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 				/* Is ring if number of faces == number of edges around vertice*/
 				i = BM_elem_index_get(f);
 				if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) {
-					nlMatrixAdd(sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
-					nlMatrixAdd(sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
+					EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
+					EIG_linear_solver_matrix_add(sys->context, idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
 				}
 				if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) {
-					nlMatrixAdd(sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
-					nlMatrixAdd(sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
+					EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
+					EIG_linear_solver_matrix_add(sys->context, idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
 				}
 				if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) {
-					nlMatrixAdd(sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
-					nlMatrixAdd(sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
+					EIG_linear_solver_matrix_add(sys->context, idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
+					EIG_linear_solver_matrix_add(sys->context, idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
 				}
 			}
 		}
@@ -368,8 +365,8 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 			idv2 = BM_elem_index_get(e->v2);
 			if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
 				i = BM_elem_index_get(e);
-				nlMatrixAdd(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
-				nlMatrixAdd(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+				EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+				EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
 			}
 		}
 	}
@@ -434,12 +431,12 @@ static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez
 		idv2 = BM_elem_index_get(e->v2);
 		vi1 = e->v1->co;
 		vi2 =  e->v2->co;
-		ve1[0] = nlGetVariable(sys->context, 0, idv1);
-		ve1[1] = nlGetVariable(sys->context, 1, idv1);
-		ve1[2] = nlGetVariable(sys->context, 2, idv1);
-		ve2[0] = nlGetVariable(sys->context, 0, idv2);
-		ve2[1] = nlGetVariable(sys->context, 1, idv2);
-		ve2[2] = nlGetVariable(sys->context, 2, idv2);
+		ve1[0] = EIG_linear_solver_variable_get(sys->context, 0, idv1);
+		ve1[1] = EIG_linear_solver_variable_get(sys->context, 1, idv1);
+		ve1[2] = EIG_linear_solver_variable_get(sys->context, 2, idv1);
+		ve2[0] = EIG_linear_solver_variable_get(sys->context, 0, idv2);
+		ve2[1] = EIG_linear_solver_variable_get(sys->context, 1, idv2);
+		ve2[2] = EIG_linear_solver_variable_get(sys->context, 2, idv2);
 		leni = len_v3v3(vi1, vi2);
 		lene = len_v3v3(ve1, ve2);
 		if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
@@ -455,13 +452,13 @@ static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez
 		m_vertex_id = BM_elem_index_get(v);
 		if (sys->zerola[m_vertex_id] == 0) {
 			if (usex) {
-				v->co[0] =  nlGetVariable(sys->context, 0, m_vertex_id);
+				v->co[0] =  EIG_linear_solver_variable_get(sys->context, 0, m_vertex_id);
 			}
 			if (usey) {
-				v->co[1] =  nlGetVariable(sys->context, 1, m_vertex_id);
+				v->co[1] =  EIG_linear_solver_variable_get(sys->context, 1, m_vertex_id);
 			}
 			if (usez) {
-				v->co[2] =  nlGetVariable(sys->context, 2, m_vertex_id);
+				v->co[2] =  EIG_linear_solver_variable_get(sys->context, 2, m_vertex_id);
 			}
 		}
 	}
@@ -501,32 +498,24 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
 	preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume");
 
 
-	sys->context = nlNewContext();
-
-	nlSolverParameteri(sys->context, NL_NB_VARIABLES, bm->totvert);
-	nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
-	nlSolverParameteri(sys->context, NL_NB_ROWS, bm->totvert);
-	nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
+	sys->context = EIG_linear_least_squares_solver_new(bm->totvert, bm->totvert, 3);
 
-	nlBegin(sys->context, NL_SYSTEM);
 	for (i = 0; i < bm->totvert; i++) {
-		nlLockVariable(sys->context, i);
+		EIG_linear_solver_variable_lock(sys->context, i);
 	}
 	BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
 		m_vertex_id = BM_elem_index_get(v);
-		nlUnlockVariable(sys->context, m_vertex_id);
-		nlSetVariable(sys->context, 0, m_vertex_id, v->co[0]);
-		nlSetVariable(sys->context, 1, m_vertex_id, v->co[1]);
-		nlSetVariable(sys->context, 2, m_vertex_id, v->co[2]);
+		EIG_linear_solver_variable_set(sys->context, 0, m_vertex_id, v->co[0]);
+		EIG_linear_solver_variable_set(sys->context, 1, m_vertex_id, v->co[1]);
+		EIG_linear_solver_variable_set(sys->context, 2, m_vertex_id, v->co[2]);
 	}
 
-	nlBegin(sys->context, NL_MATRIX);
 	init_laplacian_matrix(sys);
 	BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
 		m_vertex_id = BM_elem_index_get(v);
-		nlRightHandSideAdd(sys->context, 0, m_vertex_id, v->co[0]);
-		nlRightHandSideAdd(sys->context, 1, m_vertex_id, v->co[1]);
-		nlRightHandSideAdd(sys->context, 2, m_vertex_id, v->co[2]);
+		EIG_linear_solver_right_hand_side_add(sys->context, 0, m_vertex_id, v->co[0]);
+		EIG_linear_solver_right_hand_side_add(sys->context, 1, m_vertex_id, v->co[1]);
+		EIG_linear_solver_right_hand_side_add(sys->context, 2, m_vertex_id, v->co[2]);
 		i = m_vertex_id;
 		if (sys->zerola[i] == 0) {
 			w = sys->vweights[i] * sys->ring_areas[i];
@@ -535,34 +524,22 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
 			sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border  * 2.0f / w;
 
 			if (!vert_is_boundary(v)) {
-				nlMatrixAdd(sys->context, i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
+				EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
 			}
 			else {
-				nlMatrixAdd(sys->context, i, i,  1.0f + lambda_border * 2.0f);
+				EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + lambda_border * 2.0f);
 			}
 		}
 		else {
-			nlMatrixAdd(sys->context, i, i, 1.0f);
+			EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f);
 		}
 	}
 	fill_laplacian_matrix(sys);
 
-	nlEnd(sys->context, NL_MATRIX);
-	nlEnd(sys->context, NL_SYSTEM);
-
-	if (nlSolve(sys->context, NL_TRUE) ) {
+	if (EIG_linear_solver_solve(sys->context) ) {
 		validate_solution(sys, usex, usey, usez, preserve_volume);
 	}
 
 	delete_laplacian_system(sys);
 }
 
-#else  /* WITH_OPENNL */
-
-#ifdef __GNUC__
-#  pragma GCC diagnostic ignored "-Wunused-parameter"
-#endif
-
-void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op) {}
-
-#endif  /* WITH_OPENNL */
diff --git a/source/blender/editors/armature/CMakeLists.txt b/source/blender/editors/armature/CMakeLists.txt
index 1ed70b3cd98..b213aca478f 100644
--- a/source/blender/editors/armature/CMakeLists.txt
+++ b/source/blender/editors/armature/CMakeLists.txt
@@ -28,6 +28,7 @@ set(INC
 	../../makesrna
 	../../windowmanager
 	../../../../intern/guardedalloc
+	../../../../intern/eigen
 	../../../../intern/glew-mx
 )
 
@@ -68,13 +69,6 @@ if(WITH_INTERNATIONAL)
 	add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-	add_definitions(-DWITH_OPENNL)
-	list(APPEND INC_SYS
-		../../../../intern/opennl/extern
-	)
-endif()
-
 add_definitions(${GL_DEFINITIONS})
 
 blender_add_lib(bf_editor_armature "${SRC}" "${INC}" "${INC_SYS}")
diff --git a/source/blender/editors/armature/SConscript b/source/blender/editors/armature/SConscript
index 9c3959ecb5b..850834d2cd4 100644
--- a/source/blender/editors/armature/SConscript
+++ b/source/blender/editors/armature/SConscript
@@ -31,9 +31,9 @@ sources = env.Glob('*.c')
 
 incs = [
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     env['BF_GLEW_INC'],
     '#/intern/glew-mx',
-    '#/intern/opennl/extern',
     '../include',
     '../../blenkernel',
     '../../blenlib',
diff --git a/source/blender/editors/armature/armature_skinning.c b/source/blender/editors/armature/armature_skinning.c
index ea1a94fbba6..c621d6f99c0 100644
--- a/source/blender/editors/armature/armature_skinning.c
+++ b/source/blender/editors/armature/armature_skinning.c
@@ -51,12 +51,10 @@
 #include "ED_armature.h"
 #include "ED_mesh.h"
 
+#include "eigen_capi.h"
 
 #include "armature_intern.h"
-
-#ifdef WITH_OPENNL
-#  include "meshlaplacian.h"
-#endif
+#include "meshlaplacian.h"
 
 #if 0
 #include "reeb.h"
@@ -401,12 +399,8 @@ static void add_verts_to_dgroups(ReportList *reports, Scene *scene, Object *ob,
 	if (heat) {
 		const char *error = NULL;
 
-#ifdef WITH_OPENNL
 		heat_bone_weighting(ob, mesh, verts, numbones, dgrouplist, dgroupflip,
 		                    root, tip, selected, &error);
-#else
-		error = "Built without OpenNL";
-#endif
 		if (error) {
 			BKE_report(reports, RPT_WARNING, error);
 		}
diff --git a/source/blender/editors/armature/meshlaplacian.c b/source/blender/editors/armature/meshlaplacian.c
index 3f27a58930d..6bf75413e81 100644
--- a/source/blender/editors/armature/meshlaplacian.c
+++ b/source/blender/editors/armature/meshlaplacian.c
@@ -46,11 +46,9 @@
 #include "ED_mesh.h"
 #include "ED_armature.h"
 
-#include "meshlaplacian.h"
-
-#ifdef WITH_OPENNL
+#include "eigen_capi.h"
 
-#include "ONL_opennl.h"
+#include "meshlaplacian.h"
 
 /* ************* XXX *************** */
 static void waitcursor(int UNUSED(val)) {}
@@ -64,7 +62,7 @@ static void error(const char *str) { printf("error: %s\n", str); }
 /************************** Laplacian System *****************************/
 
 struct LaplacianSystem {
-	NLContext *context;  /* opennl context */
+	LinearSolver *context;  /* linear solver */
 
 	int totvert, totface;
 
@@ -76,7 +74,7 @@ struct LaplacianSystem {
 
 	int areaweights;        /* use area in cotangent weights? */
 	int storeweights;       /* store cotangent weights in fweights */
-	int nlbegun;            /* nlBegin(NL_SYSTEM/NL_MATRIX) done */
+	bool variablesdone;     /* variables set in linear system */
 
 	EdgeHash *edgehash;     /* edge hash for construction */
 
@@ -182,18 +180,18 @@ static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int
 	t2 = cotangent_tri_weight_v3(v2, v3, v1) / laplacian_edge_count(sys->edgehash, i3, i1);
 	t3 = cotangent_tri_weight_v3(v3, v1, v2) / laplacian_edge_count(sys->edgehash, i1, i2);
 
-	nlMatrixAdd(sys->context, i1, i1, (t2 + t3) * varea[i1]);
-	nlMatrixAdd(sys->context, i2, i2, (t1 + t3) * varea[i2]);
-	nlMatrixAdd(sys->context, i3, i3, (t1 + t2) * varea[i3]);
+	EIG_linear_solver_matrix_add(sys->context, i1, i1, (t2 + t3) * varea[i1]);
+	EIG_linear_solver_matrix_add(sys->context, i2, i2, (t1 + t3) * varea[i2]);
+	EIG_linear_solver_matrix_add(sys->context, i3, i3, (t1 + t2) * varea[i3]);
 
-	nlMatrixAdd(sys->context, i1, i2, -t3 * varea[i1]);
-	nlMatrixAdd(sys->context, i2, i1, -t3 * varea[i2]);
+	EIG_linear_solver_matrix_add(sys->context, i1, i2, -t3 * varea[i1]);
+	EIG_linear_solver_matrix_add(sys->context, i2, i1, -t3 * varea[i2]);
 
-	nlMatrixAdd(sys->context, i2, i3, -t1 * varea[i2]);
-	nlMatrixAdd(sys->context, i3, i2, -t1 * varea[i3]);
+	EIG_linear_solver_matrix_add(sys->context, i2, i3, -t1 * varea[i2]);
+	EIG_linear_solver_matrix_add(sys->context, i3, i2, -t1 * varea[i3]);
 
-	nlMatrixAdd(sys->context, i3, i1, -t2 * varea[i3]);
-	nlMatrixAdd(sys->context, i1, i3, -t2 * varea[i1]);
+	EIG_linear_solver_matrix_add(sys->context, i3, i1, -t2 * varea[i3]);
+	EIG_linear_solver_matrix_add(sys->context, i1, i3, -t2 * varea[i1]);
 
 	if (sys->storeweights) {
 		sys->fweights[f][0] = t1 * varea[i1];
@@ -218,11 +216,11 @@ static LaplacianSystem *laplacian_system_construct_begin(int totvert, int totfac
 	sys->areaweights = 1;
 	sys->storeweights = 0;
 
-	/* create opennl context */
-	sys->context = nlNewContext();
-	nlSolverParameteri(sys->context, NL_NB_VARIABLES, totvert);
+	/* create linear solver */
 	if (lsq)
-		nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
+		sys->context = EIG_linear_least_squares_solver_new(0, totvert, 1);
+	else
+		sys->context = EIG_linear_solver_new(0, totvert, 1);
 
 	return sys;
 }
@@ -272,7 +270,7 @@ static void laplacian_system_construct_end(LaplacianSystem *sys)
 
 		/* for heat weighting */
 		if (sys->heat.H)
-			nlMatrixAdd(sys->context, a, a, sys->heat.H[a]);
+			EIG_linear_solver_matrix_add(sys->context, a, a, sys->heat.H[a]);
 	}
 
 	if (sys->storeweights)
@@ -301,7 +299,7 @@ static void laplacian_system_delete(LaplacianSystem *sys)
 	if (sys->faces) MEM_freeN(sys->faces);
 	if (sys->fweights) MEM_freeN(sys->fweights);
 
-	nlDeleteContext(sys->context);
+	EIG_linear_solver_delete(sys->context);
 	MEM_freeN(sys);
 }
 
@@ -309,42 +307,37 @@ void laplacian_begin_solve(LaplacianSystem *sys, int index)
 {
 	int a;
 
-	if (!sys->nlbegun) {
-		nlBegin(sys->context, NL_SYSTEM);
-
+	if (!sys->variablesdone) {
 		if (index >= 0) {
 			for (a = 0; a < sys->totvert; a++) {
 				if (sys->vpinned[a]) {
-					nlSetVariable(sys->context, 0, a, sys->verts[a][index]);
-					nlLockVariable(sys->context, a);
+					EIG_linear_solver_variable_set(sys->context, 0, a, sys->verts[a][index]);
+					EIG_linear_solver_variable_lock(sys->context, a);
 				}
 			}
 		}
 
-		nlBegin(sys->context, NL_MATRIX);
-		sys->nlbegun = 1;
+		sys->variablesdone = true;
 	}
 }
 
 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
 {
-	nlRightHandSideAdd(sys->context, 0, v, value);
+	EIG_linear_solver_right_hand_side_add(sys->context, 0, v, value);
 }
 
 int laplacian_system_solve(LaplacianSystem *sys)
 {
-	nlEnd(sys->context, NL_MATRIX);
-	nlEnd(sys->context, NL_SYSTEM);
-	sys->nlbegun = 0;
+	sys->variablesdone = false;
 
-	//nlPrintMatrix(sys->context, );
+	//EIG_linear_solver_print_matrix(sys->context, );
 
-	return nlSolve(sys->context, NL_TRUE);
+	return EIG_linear_solver_solve(sys->context);
 }
 
 float laplacian_system_get_solution(LaplacianSystem *sys, int v)
 {
-	return nlGetVariable(sys->context, 0, v);
+	return EIG_linear_solver_variable_get(sys->context, 0, v);
 }
 
 /************************* Heat Bone Weighting ******************************/
@@ -1284,7 +1277,7 @@ static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y,
 	return totweight;
 }
 
-static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context, int x, int y, int z)
+static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z)
 {
 	MDefBoundIsect *isect;
 	float weight, totweight;
@@ -1294,7 +1287,7 @@ static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context,
 	if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
 		return;
 
-	nlMatrixAdd(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
+	EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
 	
 	totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
 	for (i = 1; i <= 6; i++) {
@@ -1305,12 +1298,12 @@ static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, NLContext *context,
 		isect = mdb->boundisect[acenter][i - 1];
 		if (!isect) {
 			weight = (1.0f / mdb->width[0]) / totweight;
-			nlMatrixAdd(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
+			EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
 		}
 	}
 }
 
-static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, NLContext *context, int x, int y, int z, int cagevert)
+static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z, int cagevert)
 {
 	MDefBoundIsect *isect;
 	float rhs, weight, totweight;
@@ -1331,7 +1324,7 @@ static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, NLContext *context, i
 		if (isect) {
 			weight = (1.0f / isect->len) / totweight;
 			rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
-			nlRightHandSideAdd(context, 0, mdb->varidx[acenter], rhs);
+			EIG_linear_solver_right_hand_side_add(context, 0, mdb->varidx[acenter], rhs);
 		}
 	}
 }
@@ -1386,7 +1379,7 @@ static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y
 
 static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind *mdb)
 {
-	NLContext *context;
+	LinearSolver *context;
 	float vec[3], gridvec[3];
 	int a, b, x, y, z, totvar;
 	char message[256];
@@ -1403,15 +1396,8 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
 	progress_bar(0, "Starting mesh deform solve");
 
-	/* setup opennl solver */
-	context = nlNewContext();
-
-	nlSolverParameteri(context, NL_NB_VARIABLES, totvar);
-	nlSolverParameteri(context, NL_NB_ROWS, totvar);
-	nlSolverParameteri(context, NL_NB_RIGHT_HAND_SIDES, 1);
-
-	nlBegin(context, NL_SYSTEM);
-	nlBegin(context, NL_MATRIX);
+	/* setup linear solver */
+	context = EIG_linear_solver_new(totvar, totvar, 1);
 
 	/* build matrix */
 	for (z = 0; z < mdb->size; z++)
@@ -1421,21 +1407,13 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
 	/* solve for each cage vert */
 	for (a = 0; a < mdb->totcagevert; a++) {
-		if (a != 0) {
-			nlBegin(context, NL_SYSTEM);
-			nlBegin(context, NL_MATRIX);
-		}
-
 		/* fill in right hand side and solve */
 		for (z = 0; z < mdb->size; z++)
 			for (y = 0; y < mdb->size; y++)
 				for (x = 0; x < mdb->size; x++)
 					meshdeform_matrix_add_rhs(mdb, context, x, y, z, a);
 
-		nlEnd(context, NL_MATRIX);
-		nlEnd(context, NL_SYSTEM);
-
-		if (nlSolve(context, NL_TRUE)) {
+		if (EIG_linear_solver_solve(context)) {
 			for (z = 0; z < mdb->size; z++)
 				for (y = 0; y < mdb->size; y++)
 					for (x = 0; x < mdb->size; x++)
@@ -1448,7 +1426,7 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 
 			for (b = 0; b < mdb->size3; b++) {
 				if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
-					mdb->phi[b] = nlGetVariable(context, 0, mdb->varidx[b]);
+					mdb->phi[b] = EIG_linear_solver_variable_get(context, 0, mdb->varidx[b]);
 				mdb->totalphi[b] += mdb->phi[b];
 			}
 
@@ -1502,7 +1480,7 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind
 	/* free */
 	MEM_freeN(mdb->varidx);
 
-	nlDeleteContext(context);
+	EIG_linear_solver_delete(context);
 }
 
 static void harmonic_coordinates_bind(Scene *UNUSED(scene), MeshDeformModifierData *mmd, MeshDeformBind *mdb)
@@ -1705,13 +1683,3 @@ void mesh_deform_bind(Scene *scene, MeshDeformModifierData *mmd, float *vertexco
 	waitcursor(0);
 }
 
-#else  /* WITH_OPENNL */
-
-#ifdef __GNUC__
-#  pragma GCC diagnostic ignored "-Wunused-parameter"
-#endif
-
-void mesh_deform_bind(Scene *scene, MeshDeformModifierData *mmd, float *vertexcos, int totvert, float cagemat[4][4]) {}
-void *modifier_mdef_compact_influences_link_kludge = modifier_mdef_compact_influences;
-
-#endif  /* WITH_OPENNL */
diff --git a/source/blender/editors/armature/reeb.c b/source/blender/editors/armature/reeb.c
index d53e5b857a0..a1c1f4f115e 100644
--- a/source/blender/editors/armature/reeb.c
+++ b/source/blender/editors/armature/reeb.c
@@ -2497,7 +2497,7 @@ int weightFromLoc(EditMesh *em, int axis)
 	return 1;
 }
 
-static void addTriangle(NLContext *context, EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
+static void addTriangle(LinearSolver *context, EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
 {
 	/* Angle opposite e1 */
 	float t1 = cotangent_tri_weight_v3(v1->co, v2->co, v3->co) / e2;
@@ -2512,23 +2512,23 @@ static void addTriangle(NLContext *context, EditVert *v1, EditVert *v2, EditVert
 	int i2 = indexData(v2);
 	int i3 = indexData(v3);
 	
-	nlMatrixAdd(context, i1, i1, t2 + t3);
-	nlMatrixAdd(context, i2, i2, t1 + t3);
-	nlMatrixAdd(context, i3, i3, t1 + t2);
+	EIG_linear_solver_matrix_add(context, i1, i1, t2 + t3);
+	EIG_linear_solver_matrix_add(context, i2, i2, t1 + t3);
+	EIG_linear_solver_matrix_add(context, i3, i3, t1 + t2);
 
-	nlMatrixAdd(context, i1, i2, -t3);
-	nlMatrixAdd(context, i2, i1, -t3);
+	EIG_linear_solver_matrix_add(context, i1, i2, -t3);
+	EIG_linear_solver_matrix_add(context, i2, i1, -t3);
 
-	nlMatrixAdd(context, i2, i3, -t1);
-	nlMatrixAdd(context, i3, i2, -t1);
+	EIG_linear_solver_matrix_add(context, i2, i3, -t1);
+	EIG_linear_solver_matrix_add(context, i3, i2, -t1);
 
-	nlMatrixAdd(context, i3, i1, -t2);
-	nlMatrixAdd(context, i1, i3, -t2);
+	EIG_linear_solver_matrix_add(context, i3, i1, -t2);
+	EIG_linear_solver_matrix_add(context, i1, i3, -t2);
 }
 
 int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 {
-	NLContext *context;
+	LinearSolver *context;
 	NLboolean success;
 	EditVert *eve;
 	EditEdge *eed;
@@ -2542,14 +2542,10 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 		totvert++;
 	}
 
-	/* Solve with openNL */
+	/* Solve */
 	
-	context = nlNewContext();
+	context = EIG_linear_solver_new(, 0, totvert, 1);
 
-	nlSolverParameteri(context, NL_NB_VARIABLES, totvert);
-
-	nlBegin(context, NL_SYSTEM);
-	
 	/* Find local extrema */
 	for (index = 0, eve = em->verts.first; eve; index++, eve = eve->next) {
 		if (eve->h == 0) {
@@ -2583,8 +2579,8 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 			if (maximum || minimum) {
 				float w = weightData(eve);
 				eve->f1 = 0;
-				nlSetVariable(context, 0, index, w);
-				nlLockVariable(context, index);
+				EIG_linear_solver_variable_set(context, 0, index, w);
+				EIG_linear_solver_variable_lock(context, index);
 			}
 			else {
 				eve->f1 = 1;
@@ -2592,8 +2588,6 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 		}
 	}
 	
-	nlBegin(context, NL_MATRIX);
-
 	/* Zero edge weight */
 	for (eed = em->edges.first; eed; eed = eed->next) {
 		eed->tmp.l = 0;
@@ -2625,23 +2619,19 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
 		}
 	}
 	
-	nlEnd(context, NL_MATRIX);
-
-	nlEnd(context, NL_SYSTEM);
-
-	success = nlSolve(context, NL_TRUE);
+	success = EIG_linear_solver_solve(context);
 
 	if (success) {
 		rval = 1;
 		for (index = 0, eve = em->verts.first; eve; index++, eve = eve->next) {
-			weightSetData(eve, nlGetVariable(context, 0, index));
+			weightSetData(eve, EIG_linear_solver_variable_get(context, 0, index));
 		}
 	}
 	else {
 		rval = 0;
 	}
 
-	nlDeleteContext(context);
+	EIG_linear_solver_delete(context);
 
 	return rval;
 }
diff --git a/source/blender/editors/uvedit/CMakeLists.txt b/source/blender/editors/uvedit/CMakeLists.txt
index a90763eed4e..543ef0e0663 100644
--- a/source/blender/editors/uvedit/CMakeLists.txt
+++ b/source/blender/editors/uvedit/CMakeLists.txt
@@ -29,6 +29,7 @@ set(INC
 	../../makesrna
 	../../windowmanager
 	../../../../intern/guardedalloc
+	../../../../intern/eigen
 	../../../../intern/glew-mx
 )
 
@@ -52,13 +53,6 @@ if(WITH_INTERNATIONAL)
 	add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-	add_definitions(-DWITH_OPENNL)
-	list(APPEND INC_SYS
-		../../../../intern/opennl/extern
-	)
-endif()
-
 add_definitions(${GL_DEFINITIONS})
 
 blender_add_lib(bf_editor_uvedit "${SRC}" "${INC}" "${INC_SYS}")
diff --git a/source/blender/editors/uvedit/SConscript b/source/blender/editors/uvedit/SConscript
index b5cccab4002..85ea7f45536 100644
--- a/source/blender/editors/uvedit/SConscript
+++ b/source/blender/editors/uvedit/SConscript
@@ -34,9 +34,9 @@ sources = env.Glob('*.c')
 
 incs = [
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     env['BF_GLEW_INC'],
     '#/intern/glew-mx',
-    '#/intern/opennl/extern',
     '../include',
     '../../blenkernel',
     '../../blenlib',
diff --git a/source/blender/editors/uvedit/uvedit_parametrizer.c b/source/blender/editors/uvedit/uvedit_parametrizer.c
index 311c152a239..e1495b617f8 100644
--- a/source/blender/editors/uvedit/uvedit_parametrizer.c
+++ b/source/blender/editors/uvedit/uvedit_parametrizer.c
@@ -44,9 +44,7 @@
 
 #include "BLI_sys_types.h"  /* for intptr_t support */
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
+#include "eigen_capi.h"
 
 /* Utils */
 
@@ -193,7 +191,7 @@ typedef struct PChart {
 
 	union PChartUnion {
 		struct PChartLscm {
-			NLContext *context;
+			LinearSolver *context;
 			float *abf_alpha;
 			PVert *pin1, *pin2;
 		} lscm;
@@ -2471,17 +2469,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 	PEdge *e;
 	int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
 	PBool success;
-	NLContext *context;
-
-	context = nlNewContext();
-	nlSolverParameteri(context, NL_NB_VARIABLES, nvar);
+	LinearSolver *context;
 
-	nlBegin(context, NL_SYSTEM);
-
-	nlBegin(context, NL_MATRIX);
+	context = EIG_linear_solver_new(0, nvar, 1);
 
 	for (i = 0; i < nvar; i++)
-		nlRightHandSideAdd(context, 0, i, sys->bInterior[i]);
+		EIG_linear_solver_right_hand_side_add(context, 0, i, sys->bInterior[i]);
 
 	for (f = chart->faces; f; f = f->nextlink) {
 		float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
@@ -2527,8 +2520,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
 			sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
 
-			nlRightHandSideAdd(context, 0, v1->u.id, j2[0][0] * beta[0]);
-			nlRightHandSideAdd(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
+			EIG_linear_solver_right_hand_side_add(context, 0, v1->u.id, j2[0][0] * beta[0]);
+			EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
 
 			row1[0] = j2[0][0] * W[0][0];
 			row2[0] = j2[0][0] * W[1][0];
@@ -2547,8 +2540,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
 			sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
 
-			nlRightHandSideAdd(context, 0, v2->u.id, j2[1][1] * beta[1]);
-			nlRightHandSideAdd(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
+			EIG_linear_solver_right_hand_side_add(context, 0, v2->u.id, j2[1][1] * beta[1]);
+			EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
 
 			row1[1] = j2[1][1] * W[0][1];
 			row2[1] = j2[1][1] * W[1][1];
@@ -2567,8 +2560,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
 			sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
 
-			nlRightHandSideAdd(context, 0, v3->u.id, j2[2][2] * beta[2]);
-			nlRightHandSideAdd(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
+			EIG_linear_solver_right_hand_side_add(context, 0, v3->u.id, j2[2][2] * beta[2]);
+			EIG_linear_solver_right_hand_side_add(context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
 
 			row1[2] = j2[2][2] * W[0][2];
 			row2[2] = j2[2][2] * W[1][2];
@@ -2592,29 +2585,25 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 					continue;
 
 				if (i == 0)
-					nlMatrixAdd(context, r, c, j2[0][i] * row1[j]);
+					EIG_linear_solver_matrix_add(context, r, c, j2[0][i] * row1[j]);
 				else
-					nlMatrixAdd(context, r + ninterior, c, j2[0][i] * row1[j]);
+					EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[0][i] * row1[j]);
 
 				if (i == 1)
-					nlMatrixAdd(context, r, c, j2[1][i] * row2[j]);
+					EIG_linear_solver_matrix_add(context, r, c, j2[1][i] * row2[j]);
 				else
-					nlMatrixAdd(context, r + ninterior, c, j2[1][i] * row2[j]);
+					EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[1][i] * row2[j]);
 
 
 				if (i == 2)
-					nlMatrixAdd(context, r, c, j2[2][i] * row3[j]);
+					EIG_linear_solver_matrix_add(context, r, c, j2[2][i] * row3[j]);
 				else
-					nlMatrixAdd(context, r + ninterior, c, j2[2][i] * row3[j]);
+					EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[2][i] * row3[j]);
 			}
 		}
 	}
 
-	nlEnd(context, NL_MATRIX);
-
-	nlEnd(context, NL_SYSTEM);
-
-	success = nlSolve(context, NL_FALSE);
+	success = EIG_linear_solver_solve(context);
 
 	if (success) {
 		for (f = chart->faces; f; f = f->nextlink) {
@@ -2625,24 +2614,24 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 			pre[0] = pre[1] = pre[2] = 0.0;
 
 			if (v1->flag & PVERT_INTERIOR) {
-				float x = nlGetVariable(context, 0, v1->u.id);
-				float x2 = nlGetVariable(context, 0, ninterior + v1->u.id);
+				float x = EIG_linear_solver_variable_get(context, 0, v1->u.id);
+				float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v1->u.id);
 				pre[0] += sys->J2dt[e1->u.id][0] * x;
 				pre[1] += sys->J2dt[e2->u.id][0] * x2;
 				pre[2] += sys->J2dt[e3->u.id][0] * x2;
 			}
 
 			if (v2->flag & PVERT_INTERIOR) {
-				float x = nlGetVariable(context, 0, v2->u.id);
-				float x2 = nlGetVariable(context, 0, ninterior + v2->u.id);
+				float x = EIG_linear_solver_variable_get(context, 0, v2->u.id);
+				float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v2->u.id);
 				pre[0] += sys->J2dt[e1->u.id][1] * x2;
 				pre[1] += sys->J2dt[e2->u.id][1] * x;
 				pre[2] += sys->J2dt[e3->u.id][1] * x2;
 			}
 
 			if (v3->flag & PVERT_INTERIOR) {
-				float x = nlGetVariable(context, 0, v3->u.id);
-				float x2 = nlGetVariable(context, 0, ninterior + v3->u.id);
+				float x = EIG_linear_solver_variable_get(context, 0, v3->u.id);
+				float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v3->u.id);
 				pre[0] += sys->J2dt[e1->u.id][2] * x2;
 				pre[1] += sys->J2dt[e2->u.id][2] * x2;
 				pre[2] += sys->J2dt[e3->u.id][2] * x;
@@ -2673,12 +2662,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
 		}
 
 		for (i = 0; i < ninterior; i++) {
-			sys->lambdaPlanar[i] += (float)nlGetVariable(context, 0, i);
-			sys->lambdaLength[i] += (float)nlGetVariable(context, 0, ninterior + i);
+			sys->lambdaPlanar[i] += (float)EIG_linear_solver_variable_get(context, 0, i);
+			sys->lambdaLength[i] += (float)EIG_linear_solver_variable_get(context, 0, ninterior + i);
 		}
 	}
 
-	nlDeleteContext(context);
+	EIG_linear_solver_delete(context);
 
 	return success;
 }
@@ -3004,12 +2993,12 @@ static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
 
 static void p_chart_lscm_load_solution(PChart *chart)
 {
-	NLContext *context = chart->u.lscm.context;
+	LinearSolver *context = chart->u.lscm.context;
 	PVert *v;
 
 	for (v = chart->verts; v; v = v->nextlink) {
-		v->uv[0] = nlGetVariable(context, 0, 2 * v->u.id);
-		v->uv[1] = nlGetVariable(context, 0, 2 * v->u.id + 1);
+		v->uv[0] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id);
+		v->uv[1] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id + 1);
 	}
 }
 
@@ -3064,16 +3053,13 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
 		for (v = chart->verts; v; v = v->nextlink)
 			v->u.id = id++;
 
-		chart->u.lscm.context = nlNewContext();
-		nlSolverParameteri(chart->u.lscm.context, NL_NB_VARIABLES, 2 * chart->nverts);
-		nlSolverParameteri(chart->u.lscm.context, NL_NB_ROWS, 2 * chart->nfaces);
-		nlSolverParameteri(chart->u.lscm.context, NL_LEAST_SQUARES, NL_TRUE);
+		chart->u.lscm.context = EIG_linear_least_squares_solver_new(2 * chart->nfaces, 2 * chart->nverts, 1);
 	}
 }
 
 static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 {
-	NLContext *context = chart->u.lscm.context;
+	LinearSolver *context = chart->u.lscm.context;
 	PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
 	PFace *f;
 	float *alpha = chart->u.lscm.abf_alpha;
@@ -3081,8 +3067,6 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 	bool flip_faces;
 	int row;
 
-	nlBegin(context, NL_SYSTEM);
-
 #if 0
 	/* TODO: make loading pins work for simplify/complexify. */
 #endif
@@ -3092,25 +3076,25 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 			p_vert_load_pin_select_uvs(handle, v);  /* reload for live */
 
 	if (chart->u.lscm.pin1) {
-		nlLockVariable(context, 2 * pin1->u.id);
-		nlLockVariable(context, 2 * pin1->u.id + 1);
-		nlLockVariable(context, 2 * pin2->u.id);
-		nlLockVariable(context, 2 * pin2->u.id + 1);
+		EIG_linear_solver_variable_lock(context, 2 * pin1->u.id);
+		EIG_linear_solver_variable_lock(context, 2 * pin1->u.id + 1);
+		EIG_linear_solver_variable_lock(context, 2 * pin2->u.id);
+		EIG_linear_solver_variable_lock(context, 2 * pin2->u.id + 1);
 
-		nlSetVariable(context, 0, 2 * pin1->u.id, pin1->uv[0]);
-		nlSetVariable(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
-		nlSetVariable(context, 0, 2 * pin2->u.id, pin2->uv[0]);
-		nlSetVariable(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
+		EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id, pin1->uv[0]);
+		EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
+		EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id, pin2->uv[0]);
+		EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
 	}
 	else {
 		/* set and lock the pins */
 		for (v = chart->verts; v; v = v->nextlink) {
 			if (v->flag & PVERT_PIN) {
-				nlLockVariable(context, 2 * v->u.id);
-				nlLockVariable(context, 2 * v->u.id + 1);
+				EIG_linear_solver_variable_lock(context, 2 * v->u.id);
+				EIG_linear_solver_variable_lock(context, 2 * v->u.id + 1);
 
-				nlSetVariable(context, 0, 2 * v->u.id, v->uv[0]);
-				nlSetVariable(context, 0, 2 * v->u.id + 1, v->uv[1]);
+				EIG_linear_solver_variable_set(context, 0, 2 * v->u.id, v->uv[0]);
+				EIG_linear_solver_variable_set(context, 0, 2 * v->u.id + 1, v->uv[1]);
 			}
 		}
 	}
@@ -3137,8 +3121,6 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 
 	/* construct matrix */
 
-	nlBegin(context, NL_MATRIX);
-
 	row = 0;
 	for (f = chart->faces; f; f = f->nextlink) {
 		PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
@@ -3185,26 +3167,22 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 		cosine = cosf(a1) * ratio;
 		sine = sina1 * ratio;
 
-		nlMatrixAdd(context, row, 2 * v1->u.id,   cosine - 1.0f);
-		nlMatrixAdd(context, row, 2 * v1->u.id + 1, -sine);
-		nlMatrixAdd(context, row, 2 * v2->u.id,   -cosine);
-		nlMatrixAdd(context, row, 2 * v2->u.id + 1, sine);
-		nlMatrixAdd(context, row, 2 * v3->u.id,   1.0);
+		EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id,   cosine - 1.0f);
+		EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, -sine);
+		EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id,   -cosine);
+		EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, sine);
+		EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id,   1.0);
 		row++;
 
-		nlMatrixAdd(context, row, 2 * v1->u.id,   sine);
-		nlMatrixAdd(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
-		nlMatrixAdd(context, row, 2 * v2->u.id,   -sine);
-		nlMatrixAdd(context, row, 2 * v2->u.id + 1, -cosine);
-		nlMatrixAdd(context, row, 2 * v3->u.id + 1, 1.0);
+		EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id,   sine);
+		EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
+		EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id,   -sine);
+		EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, -cosine);
+		EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id + 1, 1.0);
 		row++;
 	}
 
-	nlEnd(context, NL_MATRIX);
-
-	nlEnd(context, NL_SYSTEM);
-
-	if (nlSolve(context, NL_TRUE)) {
+	if (EIG_linear_solver_solve(context)) {
 		p_chart_lscm_load_solution(chart);
 		return P_TRUE;
 	}
@@ -3221,7 +3199,7 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 static void p_chart_lscm_end(PChart *chart)
 {
 	if (chart->u.lscm.context)
-		nlDeleteContext(chart->u.lscm.context);
+		EIG_linear_solver_delete(chart->u.lscm.context);
 	
 	if (chart->u.lscm.abf_alpha) {
 		MEM_freeN(chart->u.lscm.abf_alpha);
@@ -4700,36 +4678,3 @@ void param_flush_restore(ParamHandle *handle)
 	}
 }
 
-#else  /* WITH_OPENNL */
-
-#ifdef __GNUC__
-#  pragma GCC diagnostic ignored "-Wunused-parameter"
-#endif
-
-/* stubs */
-void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
-                    ParamKey *vkeys, float **co, float **uv,
-                    ParamBool *pin, ParamBool *select, float normal[3]) {}
-void param_edge_set_seam(ParamHandle *handle,
-                         ParamKey *vkeys) {}
-void param_aspect_ratio(ParamHandle *handle, float aspx, float aspy) {}
-ParamHandle *param_construct_begin(void) { return NULL; }
-void param_construct_end(ParamHandle *handle, ParamBool fill, ParamBool impl) {}
-void param_delete(ParamHandle *handle) {}
-
-void param_stretch_begin(ParamHandle *handle) {}
-void param_stretch_blend(ParamHandle *handle, float blend) {}
-void param_stretch_iter(ParamHandle *handle) {}
-void param_stretch_end(ParamHandle *handle) {}
-
-void param_pack(ParamHandle *handle, float margin, bool do_rotate) {}
-void param_average(ParamHandle *handle) {}
-
-void param_flush(ParamHandle *handle) {}
-void param_flush_restore(ParamHandle *handle) {}
-
-void param_lscm_begin(ParamHandle *handle, ParamBool live, ParamBool abf) {}
-void param_lscm_solve(ParamHandle *handle) {}
-void param_lscm_end(ParamHandle *handle) {}
-
-#endif  /* WITH_OPENNL */
diff --git a/source/blender/modifiers/CMakeLists.txt b/source/blender/modifiers/CMakeLists.txt
index ad230dede24..0de7676e8f8 100644
--- a/source/blender/modifiers/CMakeLists.txt
+++ b/source/blender/modifiers/CMakeLists.txt
@@ -37,6 +37,7 @@ set(INC
 	../render/extern/include
 	../../../intern/elbeem/extern
 	../../../intern/guardedalloc
+	../../../intern/eigen
 )
 
 set(INC_SYS
@@ -144,13 +145,6 @@ if(WITH_INTERNATIONAL)
 	add_definitions(-DWITH_INTERNATIONAL)
 endif()
 
-if(WITH_OPENNL)
-	add_definitions(-DWITH_OPENNL)
-	list(APPEND INC_SYS
-		../../../intern/opennl/extern
-	)
-endif()
-
 if(WITH_OPENSUBDIV)
 	add_definitions(-DWITH_OPENSUBDIV)
 endif()
diff --git a/source/blender/modifiers/SConscript b/source/blender/modifiers/SConscript
index 7be295aa1a0..30007118562 100644
--- a/source/blender/modifiers/SConscript
+++ b/source/blender/modifiers/SConscript
@@ -33,8 +33,8 @@ incs = [
     '.',
     './intern', 
     '#/intern/guardedalloc',
+    '#/intern/eigen',
     '#/intern/elbeem/extern',
-    '#/intern/opennl/extern',
     '../render/extern/include',
     '../bmesh',
     '../include',
diff --git a/source/blender/modifiers/intern/MOD_laplaciandeform.c b/source/blender/modifiers/intern/MOD_laplaciandeform.c
index fdaacc7cd9e..d4f02d923d3 100644
--- a/source/blender/modifiers/intern/MOD_laplaciandeform.c
+++ b/source/blender/modifiers/intern/MOD_laplaciandeform.c
@@ -42,6 +42,7 @@
 
 #include "MOD_util.h"
 
+#include "eigen_capi.h"
 
 enum {
 	LAPDEFORM_SYSTEM_NOT_CHANGE = 0,
@@ -54,10 +55,6 @@ enum {
 	LAPDEFORM_SYSTEM_CHANGE_NOT_VALID_GROUP,
 };
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
-
 typedef struct LaplacianSystem {
 	bool is_matrix_computed;
 	bool has_solution;
@@ -75,7 +72,7 @@ typedef struct LaplacianSystem {
 	int *unit_verts;			/* Unit vectors of projected edges onto the plane orthogonal to n */
 	int *ringf_indices;			/* Indices of faces per vertex */
 	int *ringv_indices;			/* Indices of neighbors(vertex) per vertex */
-	NLContext *context;			/* System for solve general implicit rotations */
+	LinearSolver *context;			/* System for solve general implicit rotations */
 	MeshElemMap *ringf_map;		/* Map of faces per vertex */
 	MeshElemMap *ringv_map;		/* Map of vertex per vertex */
 } LaplacianSystem;
@@ -134,7 +131,7 @@ static void deleteLaplacianSystem(LaplacianSystem *sys)
 	MEM_SAFE_FREE(sys->ringv_map);
 
 	if (sys->context) {
-		nlDeleteContext(sys->context);
+		EIG_linear_solver_delete(sys->context);
 	}
 	MEM_SAFE_FREE(sys);
 }
@@ -283,9 +280,9 @@ static void initLaplacianMatrix(LaplacianSystem *sys)
 			sys->delta[idv[0]][1] -= v3[1] * w3;
 			sys->delta[idv[0]][2] -= v3[2] * w3;
 
-			nlMatrixAdd(sys->context, idv[0], idv[1], -w2);
-			nlMatrixAdd(sys->context, idv[0], idv[2], -w3);
-			nlMatrixAdd(sys->context, idv[0], idv[0], w2 + w3);
+			EIG_linear_solver_matrix_add(sys->context, idv[0], idv[1], -w2);
+			EIG_linear_solver_matrix_add(sys->context, idv[0], idv[2], -w3);
+			EIG_linear_solver_matrix_add(sys->context, idv[0], idv[0], w2 + w3);
 		}
 	}
 }
@@ -338,9 +335,9 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
 		beta = dot_v3v3(uij, di);
 		gamma = dot_v3v3(e2, di);
 
-		pi[0] = nlGetVariable(sys->context, 0, i);
-		pi[1] = nlGetVariable(sys->context, 1, i);
-		pi[2] = nlGetVariable(sys->context, 2, i);
+		pi[0] = EIG_linear_solver_variable_get(sys->context, 0, i);
+		pi[1] = EIG_linear_solver_variable_get(sys->context, 1, i);
+		pi[2] = EIG_linear_solver_variable_get(sys->context, 2, i);
 		zero_v3(ni);
 		num_fni = 0;
 		num_fni = sys->ringf_map[i].count;
@@ -349,9 +346,9 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
 			fidn = sys->ringf_map[i].indices;
 			vin = sys->tris[fidn[fi]];
 			for (j = 0; j < 3; j++) {
-				vn[j][0] = nlGetVariable(sys->context, 0, vin[j]);
-				vn[j][1] = nlGetVariable(sys->context, 1, vin[j]);
-				vn[j][2] = nlGetVariable(sys->context, 2, vin[j]);
+				vn[j][0] = EIG_linear_solver_variable_get(sys->context, 0, vin[j]);
+				vn[j][1] = EIG_linear_solver_variable_get(sys->context, 1, vin[j]);
+				vn[j][2] = EIG_linear_solver_variable_get(sys->context, 2, vin[j]);
 				if (vin[j] == sys->unit_verts[i]) {
 					copy_v3_v3(pj, vn[j]);
 				}
@@ -372,14 +369,14 @@ static void rotateDifferentialCoordinates(LaplacianSystem *sys)
 		fni[2] = alpha * ni[2] + beta * uij[2] + gamma * e2[2];
 
 		if (len_squared_v3(fni) > FLT_EPSILON) {
-			nlRightHandSideSet(sys->context, 0, i, fni[0]);
-			nlRightHandSideSet(sys->context, 1, i, fni[1]);
-			nlRightHandSideSet(sys->context, 2, i, fni[2]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, i, fni[0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, i, fni[1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, i, fni[2]);
 		}
 		else {
-			nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
-			nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
-			nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, i, sys->delta[i][0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, i, sys->delta[i][1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, i, sys->delta[i][2]);
 		}
 	}
 }
@@ -390,75 +387,59 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
 	n = sys->total_verts;
 	na = sys->total_anchors;
 
-#ifdef OPENNL_THREADING_HACK
-	modifier_opennl_lock();
-#endif
-
 	if (!sys->is_matrix_computed) {
-		sys->context = nlNewContext();
+		sys->context = EIG_linear_least_squares_solver_new(n + na, n, 3);
 
-		nlSolverParameteri(sys->context, NL_NB_VARIABLES, n);
-		nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
-		nlSolverParameteri(sys->context, NL_NB_ROWS, n + na);
-		nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
-		nlBegin(sys->context, NL_SYSTEM);
 		for (i = 0; i < n; i++) {
-			nlSetVariable(sys->context, 0, i, sys->co[i][0]);
-			nlSetVariable(sys->context, 1, i, sys->co[i][1]);
-			nlSetVariable(sys->context, 2, i, sys->co[i][2]);
+			EIG_linear_solver_variable_set(sys->context, 0, i, sys->co[i][0]);
+			EIG_linear_solver_variable_set(sys->context, 1, i, sys->co[i][1]);
+			EIG_linear_solver_variable_set(sys->context, 2, i, sys->co[i][2]);
 		}
 		for (i = 0; i < na; i++) {
 			vid = sys->index_anchors[i];
-			nlSetVariable(sys->context, 0, vid, vertexCos[vid][0]);
-			nlSetVariable(sys->context, 1, vid, vertexCos[vid][1]);
-			nlSetVariable(sys->context, 2, vid, vertexCos[vid][2]);
+			EIG_linear_solver_variable_set(sys->context, 0, vid, vertexCos[vid][0]);
+			EIG_linear_solver_variable_set(sys->context, 1, vid, vertexCos[vid][1]);
+			EIG_linear_solver_variable_set(sys->context, 2, vid, vertexCos[vid][2]);
 		}
-		nlBegin(sys->context, NL_MATRIX);
 
 		initLaplacianMatrix(sys);
 		computeImplictRotations(sys);
 
 		for (i = 0; i < n; i++) {
-			nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
-			nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
-			nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, i, sys->delta[i][0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, i, sys->delta[i][1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, i, sys->delta[i][2]);
 		}
 		for (i = 0; i < na; i++) {
 			vid = sys->index_anchors[i];
-			nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
-			nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
-			nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
-			nlMatrixAdd(sys->context, n + i, vid, 1.0f);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, n + i, vertexCos[vid][0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, n + i, vertexCos[vid][1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, n + i, vertexCos[vid][2]);
+			EIG_linear_solver_matrix_add(sys->context, n + i, vid, 1.0f);
 		}
-		nlEnd(sys->context, NL_MATRIX);
-		nlEnd(sys->context, NL_SYSTEM);
-		if (nlSolve(sys->context, NL_TRUE)) {
+		if (EIG_linear_solver_solve(sys->context)) {
 			sys->has_solution = true;
 
 			for (j = 1; j <= sys->repeat; j++) {
-				nlBegin(sys->context, NL_SYSTEM);
-				nlBegin(sys->context, NL_MATRIX);
 				rotateDifferentialCoordinates(sys);
 
 				for (i = 0; i < na; i++) {
 					vid = sys->index_anchors[i];
-					nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
-					nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
-					nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
+					EIG_linear_solver_right_hand_side_add(sys->context, 0, n + i, vertexCos[vid][0]);
+					EIG_linear_solver_right_hand_side_add(sys->context, 1, n + i, vertexCos[vid][1]);
+					EIG_linear_solver_right_hand_side_add(sys->context, 2, n + i, vertexCos[vid][2]);
 				}
 
-				nlEnd(sys->context, NL_MATRIX);
-				nlEnd(sys->context, NL_SYSTEM);
-				if (!nlSolve(sys->context, NL_FALSE)) {
+				if (!EIG_linear_solver_solve(sys->context)) {
 					sys->has_solution = false;
 					break;
 				}
 			}
 			if (sys->has_solution) {
 				for (vid = 0; vid < sys->total_verts; vid++) {
-					vertexCos[vid][0] = nlGetVariable(sys->context, 0, vid);
-					vertexCos[vid][1] = nlGetVariable(sys->context, 1, vid);
-					vertexCos[vid][2] = nlGetVariable(sys->context, 2, vid);
+					vertexCos[vid][0] = EIG_linear_solver_variable_get(sys->context, 0, vid);
+					vertexCos[vid][1] = EIG_linear_solver_variable_get(sys->context, 1, vid);
+					vertexCos[vid][2] = EIG_linear_solver_variable_get(sys->context, 2, vid);
 				}
 			}
 			else {
@@ -473,49 +454,40 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
 
 	}
 	else if (sys->has_solution) {
-		nlBegin(sys->context, NL_SYSTEM);
-		nlBegin(sys->context, NL_MATRIX);
-
 		for (i = 0; i < n; i++) {
-			nlRightHandSideSet(sys->context, 0, i, sys->delta[i][0]);
-			nlRightHandSideSet(sys->context, 1, i, sys->delta[i][1]);
-			nlRightHandSideSet(sys->context, 2, i, sys->delta[i][2]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, i, sys->delta[i][0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, i, sys->delta[i][1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, i, sys->delta[i][2]);
 		}
 		for (i = 0; i < na; i++) {
 			vid = sys->index_anchors[i];
-			nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
-			nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
-			nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
-			nlMatrixAdd(sys->context, n + i, vid, 1.0f);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, n + i, vertexCos[vid][0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, n + i, vertexCos[vid][1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, n + i, vertexCos[vid][2]);
+			EIG_linear_solver_matrix_add(sys->context, n + i, vid, 1.0f);
 		}
 
-		nlEnd(sys->context, NL_MATRIX);
-		nlEnd(sys->context, NL_SYSTEM);
-		if (nlSolve(sys->context, NL_FALSE)) {
+		if (EIG_linear_solver_solve(sys->context)) {
 			sys->has_solution = true;
 			for (j = 1; j <= sys->repeat; j++) {
-				nlBegin(sys->context, NL_SYSTEM);
-				nlBegin(sys->context, NL_MATRIX);
 				rotateDifferentialCoordinates(sys);
 
 				for (i = 0; i < na; i++) {
 					vid = sys->index_anchors[i];
-					nlRightHandSideSet(sys->context, 0, n + i, vertexCos[vid][0]);
-					nlRightHandSideSet(sys->context, 1, n + i, vertexCos[vid][1]);
-					nlRightHandSideSet(sys->context, 2, n + i, vertexCos[vid][2]);
+					EIG_linear_solver_right_hand_side_add(sys->context, 0, n + i, vertexCos[vid][0]);
+					EIG_linear_solver_right_hand_side_add(sys->context, 1, n + i, vertexCos[vid][1]);
+					EIG_linear_solver_right_hand_side_add(sys->context, 2, n + i, vertexCos[vid][2]);
 				}
-				nlEnd(sys->context, NL_MATRIX);
-				nlEnd(sys->context, NL_SYSTEM);
-				if (!nlSolve(sys->context, NL_FALSE)) {
+				if (!EIG_linear_solver_solve(sys->context)) {
 					sys->has_solution = false;
 					break;
 				}
 			}
 			if (sys->has_solution) {
 				for (vid = 0; vid < sys->total_verts; vid++) {
-					vertexCos[vid][0] = nlGetVariable(sys->context, 0, vid);
-					vertexCos[vid][1] = nlGetVariable(sys->context, 1, vid);
-					vertexCos[vid][2] = nlGetVariable(sys->context, 2, vid);
+					vertexCos[vid][0] = EIG_linear_solver_variable_get(sys->context, 0, vid);
+					vertexCos[vid][1] = EIG_linear_solver_variable_get(sys->context, 1, vid);
+					vertexCos[vid][2] = EIG_linear_solver_variable_get(sys->context, 2, vid);
 				}
 			}
 			else {
@@ -526,10 +498,6 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
 			sys->has_solution = false;
 		}
 	}
-
-#ifdef OPENNL_THREADING_HACK
-	modifier_opennl_unlock();
-#endif
 }
 
 static bool isValidVertexGroup(LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm)
@@ -720,15 +688,6 @@ static void LaplacianDeformModifier_do(
 	}
 }
 
-#else  /* WITH_OPENNL */
-static void LaplacianDeformModifier_do(
-        LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm,
-        float (*vertexCos)[3], int numVerts)
-{
-	UNUSED_VARS(lmd, ob, dm, vertexCos, numVerts);
-}
-#endif  /* WITH_OPENNL */
-
 static void initData(ModifierData *md)
 {
 	LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
@@ -792,12 +751,10 @@ static void deformVertsEM(
 static void freeData(ModifierData *md)
 {
 	LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
-#ifdef WITH_OPENNL
 	LaplacianSystem *sys = (LaplacianSystem *)lmd->cache_system;
 	if (sys) {
 		deleteLaplacianSystem(sys);
 	}
-#endif
 	MEM_SAFE_FREE(lmd->vertexco);
 	lmd->total_verts = 0;
 }
diff --git a/source/blender/modifiers/intern/MOD_laplaciansmooth.c b/source/blender/modifiers/intern/MOD_laplaciansmooth.c
index 189ceb11d08..f1216ff462a 100644
--- a/source/blender/modifiers/intern/MOD_laplaciansmooth.c
+++ b/source/blender/modifiers/intern/MOD_laplaciansmooth.c
@@ -43,9 +43,7 @@
 
 #include "MOD_util.h"
 
-#ifdef WITH_OPENNL
-
-#include "ONL_opennl.h"
+#include "eigen_capi.h"
 
 #if 0
 #define MOD_LAPLACIANSMOOTH_MAX_EDGE_PERCENTAGE 1.8f
@@ -71,7 +69,7 @@ struct BLaplacianSystem {
 	const MPoly *mpoly;
 	const MLoop *mloop;
 	const MEdge *medges;
-	NLContext *context;
+	LinearSolver *context;
 
 	/*Data*/
 	float min_area;
@@ -104,7 +102,7 @@ static void delete_laplacian_system(LaplacianSystem *sys)
 	MEM_SAFE_FREE(sys->zerola);
 
 	if (sys->context) {
-		nlDeleteContext(sys->context);
+		EIG_linear_solver_delete(sys->context);
 	}
 	sys->vertexCos = NULL;
 	sys->mpoly = NULL;
@@ -300,16 +298,16 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 
 			/* Is ring if number of faces == number of edges around vertice*/
 			if (sys->numNeEd[l_curr->v] == sys->numNeFa[l_curr->v] && sys->zerola[l_curr->v] == 0) {
-				nlMatrixAdd(sys->context, l_curr->v, l_next->v, sys->fweights[l_curr_index][2] * sys->vweights[l_curr->v]);
-				nlMatrixAdd(sys->context, l_curr->v, l_prev->v, sys->fweights[l_curr_index][1] * sys->vweights[l_curr->v]);
+				EIG_linear_solver_matrix_add(sys->context, l_curr->v, l_next->v, sys->fweights[l_curr_index][2] * sys->vweights[l_curr->v]);
+				EIG_linear_solver_matrix_add(sys->context, l_curr->v, l_prev->v, sys->fweights[l_curr_index][1] * sys->vweights[l_curr->v]);
 			}
 			if (sys->numNeEd[l_next->v] == sys->numNeFa[l_next->v] && sys->zerola[l_next->v] == 0) {
-				nlMatrixAdd(sys->context, l_next->v, l_curr->v, sys->fweights[l_curr_index][2] * sys->vweights[l_next->v]);
-				nlMatrixAdd(sys->context, l_next->v, l_prev->v, sys->fweights[l_curr_index][0] * sys->vweights[l_next->v]);
+				EIG_linear_solver_matrix_add(sys->context, l_next->v, l_curr->v, sys->fweights[l_curr_index][2] * sys->vweights[l_next->v]);
+				EIG_linear_solver_matrix_add(sys->context, l_next->v, l_prev->v, sys->fweights[l_curr_index][0] * sys->vweights[l_next->v]);
 			}
 			if (sys->numNeEd[l_prev->v] == sys->numNeFa[l_prev->v] && sys->zerola[l_prev->v] == 0) {
-				nlMatrixAdd(sys->context, l_prev->v, l_curr->v, sys->fweights[l_curr_index][1] * sys->vweights[l_prev->v]);
-				nlMatrixAdd(sys->context, l_prev->v, l_next->v, sys->fweights[l_curr_index][0] * sys->vweights[l_prev->v]);
+				EIG_linear_solver_matrix_add(sys->context, l_prev->v, l_curr->v, sys->fweights[l_curr_index][1] * sys->vweights[l_prev->v]);
+				EIG_linear_solver_matrix_add(sys->context, l_prev->v, l_next->v, sys->fweights[l_curr_index][0] * sys->vweights[l_prev->v]);
 			}
 		}
 	}
@@ -323,8 +321,8 @@ static void fill_laplacian_matrix(LaplacianSystem *sys)
 		    sys->zerola[idv1] == 0 &&
 		    sys->zerola[idv2] == 0)
 		{
-			nlMatrixAdd(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
-			nlMatrixAdd(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+			EIG_linear_solver_matrix_add(sys->context, idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+			EIG_linear_solver_matrix_add(sys->context, idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
 		}
 	}
 }
@@ -342,13 +340,13 @@ static void validate_solution(LaplacianSystem *sys, short flag, float lambda, fl
 		if (sys->zerola[i] == 0) {
 			lam = sys->numNeEd[i] == sys->numNeFa[i] ? (lambda >= 0.0f ? 1.0f : -1.0f) : (lambda_border >= 0.0f ? 1.0f : -1.0f);
 			if (flag & MOD_LAPLACIANSMOOTH_X) {
-				sys->vertexCos[i][0] += lam * ((float)nlGetVariable(sys->context, 0, i) - sys->vertexCos[i][0]);
+				sys->vertexCos[i][0] += lam * ((float)EIG_linear_solver_variable_get(sys->context, 0, i) - sys->vertexCos[i][0]);
 			}
 			if (flag & MOD_LAPLACIANSMOOTH_Y) {
-				sys->vertexCos[i][1] += lam * ((float)nlGetVariable(sys->context, 1, i) - sys->vertexCos[i][1]);
+				sys->vertexCos[i][1] += lam * ((float)EIG_linear_solver_variable_get(sys->context, 1, i) - sys->vertexCos[i][1]);
 			}
 			if (flag & MOD_LAPLACIANSMOOTH_Z) {
-				sys->vertexCos[i][2] += lam * ((float)nlGetVariable(sys->context, 2, i) - sys->vertexCos[i][2]);
+				sys->vertexCos[i][2] += lam * ((float)EIG_linear_solver_variable_get(sys->context, 2, i) - sys->vertexCos[i][2]);
 			}
 		}
 	}
@@ -386,24 +384,15 @@ static void laplaciansmoothModifier_do(
 	sys->vert_centroid[2] = 0.0f;
 	memset_laplacian_system(sys, 0);
 
-#ifdef OPENNL_THREADING_HACK
-	modifier_opennl_lock();
-#endif
-
-	sys->context = nlNewContext();
-	nlSolverParameteri(sys->context, NL_NB_VARIABLES, numVerts);
-	nlSolverParameteri(sys->context, NL_LEAST_SQUARES, NL_TRUE);
-	nlSolverParameteri(sys->context, NL_NB_ROWS, numVerts);
-	nlSolverParameteri(sys->context, NL_NB_RIGHT_HAND_SIDES, 3);
+	sys->context = EIG_linear_least_squares_solver_new(numVerts, numVerts, 3);
 
 	init_laplacian_matrix(sys);
 
 	for (iter = 0; iter < smd->repeat; iter++) {
-		nlBegin(sys->context, NL_SYSTEM);
 		for (i = 0; i < numVerts; i++) {
-			nlSetVariable(sys->context, 0, i, vertexCos[i][0]);
-			nlSetVariable(sys->context, 1, i, vertexCos[i][1]);
-			nlSetVariable(sys->context, 2, i, vertexCos[i][2]);
+			EIG_linear_solver_variable_set(sys->context, 0, i, vertexCos[i][0]);
+			EIG_linear_solver_variable_set(sys->context, 1, i, vertexCos[i][1]);
+			EIG_linear_solver_variable_set(sys->context, 2, i, vertexCos[i][2]);
 			if (iter == 0) {
 				add_v3_v3(sys->vert_centroid, vertexCos[i]);
 			}
@@ -412,12 +401,11 @@ static void laplaciansmoothModifier_do(
 			mul_v3_fl(sys->vert_centroid, 1.0f / (float)numVerts);
 		}
 
-		nlBegin(sys->context, NL_MATRIX);
 		dv = dvert;
 		for (i = 0; i < numVerts; i++) {
-			nlRightHandSideSet(sys->context, 0, i, vertexCos[i][0]);
-			nlRightHandSideSet(sys->context, 1, i, vertexCos[i][1]);
-			nlRightHandSideSet(sys->context, 2, i, vertexCos[i][2]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 0, i, vertexCos[i][0]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 1, i, vertexCos[i][1]);
+			EIG_linear_solver_right_hand_side_add(sys->context, 2, i, vertexCos[i][2]);
 			if (iter == 0) {
 				if (dv) {
 					wpaint = defvert_find_weight(dv, defgrp_index);
@@ -434,10 +422,10 @@ static void laplaciansmoothModifier_do(
 						w = sys->vlengths[i];
 						sys->vlengths[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda_border) * wpaint * 2.0f / w;
 						if (sys->numNeEd[i] == sys->numNeFa[i]) {
-							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda) * wpaint);
+							EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + fabsf(smd->lambda) * wpaint);
 						}
 						else {
-							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
+							EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
 						}
 					}
 					else {
@@ -447,15 +435,15 @@ static void laplaciansmoothModifier_do(
 						sys->vlengths[i] = (w == 0.0f) ? 0.0f : -fabsf(smd->lambda_border) * wpaint * 2.0f / w;
 
 						if (sys->numNeEd[i] == sys->numNeFa[i]) {
-							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda) * wpaint / (4.0f * sys->ring_areas[i]));
+							EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + fabsf(smd->lambda) * wpaint / (4.0f * sys->ring_areas[i]));
 						}
 						else {
-							nlMatrixAdd(sys->context, i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
+							EIG_linear_solver_matrix_add(sys->context, i, i,  1.0f + fabsf(smd->lambda_border) * wpaint * 2.0f);
 						}
 					}
 				}
 				else {
-					nlMatrixAdd(sys->context, i, i, 1.0f);
+					EIG_linear_solver_matrix_add(sys->context, i, i, 1.0f);
 				}
 			}
 		}
@@ -464,32 +452,16 @@ static void laplaciansmoothModifier_do(
 			fill_laplacian_matrix(sys);
 		}
 
-		nlEnd(sys->context, NL_MATRIX);
-		nlEnd(sys->context, NL_SYSTEM);
-
-		if (nlSolve(sys->context, NL_TRUE)) {
+		if (EIG_linear_solver_solve(sys->context)) {
 			validate_solution(sys, smd->flag, smd->lambda, smd->lambda_border);
 		}
 	}
-	nlDeleteContext(sys->context);
+	EIG_linear_solver_delete(sys->context);
 	sys->context = NULL;
 
-#ifdef OPENNL_THREADING_HACK
-	modifier_opennl_unlock();
-#endif
-
 	delete_laplacian_system(sys);
 }
 
-#else  /* WITH_OPENNL */
-static void laplaciansmoothModifier_do(
-        LaplacianSmoothModifierData *smd, Object *ob, DerivedMesh *dm,
-        float (*vertexCos)[3], int numVerts)
-{
-	UNUSED_VARS(smd, ob, dm, vertexCos, numVerts);
-}
-#endif  /* WITH_OPENNL */
-
 static void init_data(ModifierData *md)
 {
 	LaplacianSmoothModifierData *smd = (LaplacianSmoothModifierData *) md;
diff --git a/source/blender/modifiers/intern/MOD_util.c b/source/blender/modifiers/intern/MOD_util.c
index be6f7af7791..f9291fb077f 100644
--- a/source/blender/modifiers/intern/MOD_util.c
+++ b/source/blender/modifiers/intern/MOD_util.c
@@ -55,10 +55,6 @@
 
 #include "MEM_guardedalloc.h"
 
-#ifdef OPENNL_THREADING_HACK
-#include "BLI_threads.h"
-#endif
-
 void modifier_init_texture(const Scene *scene, Tex *tex)
 {
 	if (!tex)
@@ -234,23 +230,6 @@ void modifier_get_vgroup(Object *ob, DerivedMesh *dm, const char *name, MDeformV
 }
 
 
-#ifdef OPENNL_THREADING_HACK
-
-static ThreadMutex opennl_context_mutex = BLI_MUTEX_INITIALIZER;
-
-void modifier_opennl_lock(void)
-{
-	BLI_mutex_lock(&opennl_context_mutex);
-}
-
-void modifier_opennl_unlock(void)
-{
-	BLI_mutex_unlock(&opennl_context_mutex);
-}
-
-#endif
-
-
 /* only called by BKE_modifier.h/modifier.c */
 void modifier_type_init(ModifierTypeInfo *types[])
 {
diff --git a/source/blender/modifiers/intern/MOD_util.h b/source/blender/modifiers/intern/MOD_util.h
index b74ff9c2a25..095d7c278df 100644
--- a/source/blender/modifiers/intern/MOD_util.h
+++ b/source/blender/modifiers/intern/MOD_util.h
@@ -52,21 +52,4 @@ struct DerivedMesh *get_dm_for_modifier(struct Object *ob, ModifierApplyFlag fla
 void modifier_get_vgroup(struct Object *ob, struct DerivedMesh *dm,
                          const char *name, struct MDeformVert **dvert, int *defgrp_index);
 
-/* XXX workaround for non-threadsafe context in OpenNL (T38403)
- * OpenNL uses global pointer for "current context", which causes
- * conflict when multiple modifiers get evaluated in threaded depgraph.
- * This is just a stupid hack to prevent assert failure / crash,
- * otherwise we'd have to modify OpenNL on a large scale.
- * OpenNL should be replaced eventually, there are other options (eigen, ceres).
- * - lukas_t
- */
-#ifdef WITH_OPENNL
-#define OPENNL_THREADING_HACK
-#endif
-
-#ifdef OPENNL_THREADING_HACK
-void modifier_opennl_lock(void);
-void modifier_opennl_unlock(void);
-#endif
-
 #endif /* __MOD_UTIL_H__ */