Eigen: fold remaining OpenNL code into intern/eigen.

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

6 files changed, 85 insertions, 200 deletions

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__ */