Eigen: fold remaining OpenNL code into intern/eigen.

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

14 files changed, 432 insertions, 707 deletions

diff --git a/intern/CMakeLists.txt b/intern/CMakeLists.txt
index 9b0f2de22f2..275524f788f 100644
--- a/intern/CMakeLists.txt
+++ b/intern/CMakeLists.txt
@@ -74,10 +74,6 @@ if(WITH_BULLET)
 	add_subdirectory(rigidbody)
 endif()
 
-if(WITH_OPENNL)
-	add_subdirectory(opennl)
-endif()
-
 if(WITH_OPENSUBDIV)
 	add_subdirectory(opensubdiv)
 endif()
diff --git a/intern/SConscript b/intern/SConscript
index 3b855d60ac8..736242102cc 100644
--- a/intern/SConscript
+++ b/intern/SConscript
@@ -37,7 +37,6 @@ SConscript(['string/SConscript',
             'itasc/SConscript',
             'eigen/SConscript',
             'opencolorio/SConscript',
-            'opennl/SConscript',
             'mikktspace/SConscript',
             'smoke/SConscript',
             'raskter/SConscript'])
diff --git a/intern/eigen/CMakeLists.txt b/intern/eigen/CMakeLists.txt
index 58964e43224..5811b71de94 100644
--- a/intern/eigen/CMakeLists.txt
+++ b/intern/eigen/CMakeLists.txt
@@ -35,9 +35,11 @@ set(SRC
 	eigen_capi.h
 
 	intern/eigenvalues.cc
+	intern/linear_solver.cc
 	intern/svd.cc
 
 	intern/eigenvalues.h
+	intern/linear_solver.h
 	intern/svd.h
 )
 
diff --git a/intern/eigen/eigen_capi.h b/intern/eigen/eigen_capi.h
index 45ee1c015ec..be42e340274 100644
--- a/intern/eigen/eigen_capi.h
+++ b/intern/eigen/eigen_capi.h
@@ -28,6 +28,7 @@
 #define __EIGEN_C_API_H__
 
 #include "intern/eigenvalues.h"
+#include "intern/linear_solver.h"
 #include "intern/svd.h"
 
 #endif  /* __EIGEN_C_API_H__ */
diff --git a/intern/eigen/intern/eigenvalues.cc b/intern/eigen/intern/eigenvalues.cc
index dcaaee8e9c2..57942a4dc55 100644
--- a/intern/eigen/intern/eigenvalues.cc
+++ b/intern/eigen/intern/eigenvalues.cc
@@ -45,7 +45,7 @@ using Eigen::Map;
 
 using Eigen::Success;
 
-bool EG3_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors)
+bool EIG_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors)
 {
 	SelfAdjointEigenSolver<MatrixXf> eigen_solver;
 
diff --git a/intern/eigen/intern/eigenvalues.h b/intern/eigen/intern/eigenvalues.h
index 93fc06c2339..5c08ab5be39 100644
--- a/intern/eigen/intern/eigenvalues.h
+++ b/intern/eigen/intern/eigenvalues.h
@@ -31,7 +31,7 @@
 extern "C" {
 #endif
 
-bool EG3_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors);
+bool EIG_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors);
 
 #ifdef __cplusplus
 }
diff --git a/intern/eigen/intern/linear_solver.cc b/intern/eigen/intern/linear_solver.cc
new file mode 100644
index 00000000000..181b278b9c0
--- /dev/null
+++ b/intern/eigen/intern/linear_solver.cc
@@ -0,0 +1,354 @@
+/*
+ *  Sparse linear solver.
+ *  Copyright (C) 2004 Bruno Levy
+ *  Copyright (C) 2005-2015 Blender Foundation
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ */
+
+#include "linear_solver.h"
+
+#include <Eigen/Sparse>
+
+#include <algorithm>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+/* Eigen data structures */
+
+typedef Eigen::SparseMatrix<double, Eigen::ColMajor> EigenSparseMatrix;
+typedef Eigen::SparseLU<EigenSparseMatrix> EigenSparseLU;
+typedef Eigen::VectorXd EigenVectorX;
+typedef Eigen::Triplet<double> EigenTriplet;
+
+/* Linear Solver data structure */
+
+struct LinearSolver
+{
+	struct Coeff
+	{
+		Coeff()
+		{
+			index = 0;
+			value = 0.0;
+		}
+
+		int index;
+		double value;
+	};
+
+	struct Variable
+	{
+		Variable()
+		{
+			memset(value, 0, sizeof(value));
+			locked = false;
+			index = 0;
+		}
+
+		double value[4];
+		bool locked;
+		int index;
+		std::vector<Coeff> a;
+	};
+
+	enum State
+	{
+		STATE_VARIABLES_CONSTRUCT,
+		STATE_MATRIX_CONSTRUCT,
+		STATE_MATRIX_SOLVED
+	};
+
+	LinearSolver(int num_rows_, int num_variables_, int num_rhs_, bool lsq_)
+	{
+		assert(num_variables_ > 0);
+		assert(num_rhs_ <= 4);
+
+		state = STATE_VARIABLES_CONSTRUCT;
+		m = 0;
+		n = 0;
+		sparseLU = NULL;
+		num_variables = num_variables_;
+		num_rhs = num_rhs_;
+		num_rows = num_rows_;
+		least_squares = lsq_;
+
+		variable.resize(num_variables);
+	}
+
+	~LinearSolver()
+	{
+		delete sparseLU;
+	}
+
+	State state;
+
+	int n;
+	int m;
+
+	std::vector<EigenTriplet> Mtriplets;
+	EigenSparseMatrix M;
+	EigenSparseMatrix MtM;
+	std::vector<EigenVectorX> b;
+	std::vector<EigenVectorX> x;
+
+	EigenSparseLU *sparseLU;
+
+	int num_variables;
+	std::vector<Variable> variable;
+
+	int num_rows;
+	int num_rhs;
+
+	bool least_squares;
+};
+
+LinearSolver *EIG_linear_solver_new(int num_rows, int num_columns, int num_rhs)
+{
+	return new LinearSolver(num_rows, num_columns, num_rhs, false);
+}
+
+LinearSolver *EIG_linear_least_squares_solver_new(int num_rows, int num_columns, int num_rhs)
+{
+	return new LinearSolver(num_rows, num_columns, num_rhs, true);
+}
+
+void EIG_linear_solver_delete(LinearSolver *solver)
+{
+	delete solver;
+}
+
+/* Variables */
+
+void EIG_linear_solver_variable_set(LinearSolver *solver, int rhs, int index, double value)
+{
+	solver->variable[index].value[rhs] = value;
+}
+
+double EIG_linear_solver_variable_get(LinearSolver *solver, int rhs, int index)
+{
+	return solver->variable[index].value[rhs];
+}
+
+void EIG_linear_solver_variable_lock(LinearSolver *solver, int index)
+{
+	if (!solver->variable[index].locked) {
+		assert(solver->state == LinearSolver::STATE_VARIABLES_CONSTRUCT);
+		solver->variable[index].locked = true;
+	}
+}
+
+static void linear_solver_variables_to_vector(LinearSolver *solver)
+{
+	int num_rhs = solver->num_rhs;
+
+	for (int i = 0; i < solver->num_variables; i++) {
+		LinearSolver::Variable* v = &solver->variable[i];
+		if (!v->locked) {
+			for (int j = 0; j < num_rhs; j++)
+				solver->x[j][v->index] = v->value[j];
+		}
+	}
+}
+
+static void linear_solver_vector_to_variables(LinearSolver *solver)
+{
+	int num_rhs = solver->num_rhs;
+
+	for (int i = 0; i < solver->num_variables; i++) {
+		LinearSolver::Variable* v = &solver->variable[i];
+		if (!v->locked) {
+			for (int j = 0; j < num_rhs; j++)
+				v->value[j] = solver->x[j][v->index];
+		}
+	}
+}
+
+/* Matrix */
+
+static void linear_solver_ensure_matrix_construct(LinearSolver *solver)
+{
+	/* transition to matrix construction if necessary */
+	if (solver->state == LinearSolver::STATE_VARIABLES_CONSTRUCT) {
+		int n = 0;
+
+		for (int i = 0; i < solver->num_variables; i++) {
+			if (solver->variable[i].locked)
+				solver->variable[i].index = ~0;
+			else
+				solver->variable[i].index = n++;
+		}
+
+		int m = (solver->num_rows == 0)? n: solver->num_rows;
+
+		solver->m = m;
+		solver->n = n;
+
+		assert(solver->least_squares || m == n);
+
+		/* reserve reasonable estimate */
+		solver->Mtriplets.clear();
+		solver->Mtriplets.reserve(std::max(m, n)*3);
+
+		solver->b.resize(solver->num_rhs);
+		solver->x.resize(solver->num_rhs);
+
+		for (int i = 0; i < solver->num_rhs; i++) {
+			solver->b[i].setZero(m);
+			solver->x[i].setZero(n);
+		}
+
+		linear_solver_variables_to_vector(solver);
+
+		solver->state = LinearSolver::STATE_MATRIX_CONSTRUCT;
+	}
+}
+
+void EIG_linear_solver_matrix_add(LinearSolver *solver, int row, int col, double value)
+{
+	if (solver->state == LinearSolver::STATE_MATRIX_SOLVED)
+		return;
+
+	linear_solver_ensure_matrix_construct(solver);
+
+	if (!solver->least_squares && solver->variable[row].locked);
+	else if (solver->variable[col].locked) {
+		if (!solver->least_squares)
+			row = solver->variable[row].index;
+
+		LinearSolver::Coeff coeff;
+		coeff.index = row;
+		coeff.value = value;
+		solver->variable[col].a.push_back(coeff);
+	}
+	else {
+		if (!solver->least_squares)
+			row = solver->variable[row].index;
+		col = solver->variable[col].index;
+
+		/* direct insert into matrix is too slow, so use triplets */
+		EigenTriplet triplet(row, col, value);
+		solver->Mtriplets.push_back(triplet);
+	}
+}
+
+/* Right hand side */
+
+void EIG_linear_solver_right_hand_side_add(LinearSolver *solver, int rhs, int index, double value)
+{
+	linear_solver_ensure_matrix_construct(solver);
+
+	if (solver->least_squares) {
+		solver->b[rhs][index] += value;
+	}
+	else if (!solver->variable[index].locked) {
+		index = solver->variable[index].index;
+		solver->b[rhs][index] += value;
+	}
+}
+
+/* Solve */
+
+bool EIG_linear_solver_solve(LinearSolver *solver)
+{
+	bool result = true;
+
+	assert(solver->state != LinearSolver::STATE_VARIABLES_CONSTRUCT);
+
+	if (solver->state == LinearSolver::STATE_MATRIX_CONSTRUCT) {
+		/* create matrix from triplets */
+		solver->M.resize(solver->m, solver->n);
+		solver->M.setFromTriplets(solver->Mtriplets.begin(), solver->Mtriplets.end());
+		solver->Mtriplets.clear();
+
+		/* create least squares matrix */
+		if (solver->least_squares)
+			solver->MtM = solver->M.transpose() * solver->M;
+
+		/* convert M to compressed column format */
+		EigenSparseMatrix& M = (solver->least_squares)? solver->MtM: solver->M;
+		M.makeCompressed();
+
+		/* perform sparse LU factorization */
+		EigenSparseLU *sparseLU = new EigenSparseLU();
+		solver->sparseLU = sparseLU;
+
+		sparseLU->compute(M);
+		result = (sparseLU->info() == Eigen::Success);
+
+		solver->state = LinearSolver::STATE_MATRIX_SOLVED;
+	}
+
+	if (result) {
+		/* solve for each right hand side */
+		for (int rhs = 0; rhs < solver->num_rhs; rhs++) {
+			/* modify for locked variables */
+			EigenVectorX& b = solver->b[rhs];
+
+			for (int i = 0; i < solver->num_variables; i++) {
+				LinearSolver::Variable *variable = &solver->variable[i];
+
+				if (variable->locked) {
+					std::vector<LinearSolver::Coeff>& a = variable->a;
+
+					for (int j = 0; j < a.size(); j++)
+						b[a[j].index] -= a[j].value*variable->value[rhs];
+				}
+			}
+
+			/* solve */
+			if (solver->least_squares) {
+				EigenVectorX Mtb = solver->M.transpose() * b;
+				solver->x[rhs] = solver->sparseLU->solve(Mtb);
+			}
+			else {
+				EigenVectorX& b = solver->b[rhs];
+				solver->x[rhs] = solver->sparseLU->solve(b);
+			}
+
+			if (solver->sparseLU->info() != Eigen::Success)
+				result = false;
+		}
+
+		if (result)
+			linear_solver_vector_to_variables(solver);
+	}
+
+	/* clear for next solve */
+	for (int rhs = 0; rhs < solver->num_rhs; rhs++)
+		solver->b[rhs].setZero(solver->m);
+
+	return result;
+}
+
+/* Debugging */
+
+void EIG_linear_solver_print_matrix(LinearSolver *solver)
+{
+	std::cout << "A:" << solver->M << std::endl;
+
+	for (int rhs = 0; rhs < solver->num_rhs; rhs++)
+		std::cout << "b " << rhs << ":" << solver->b[rhs] << std::endl;
+
+	if (solver->MtM.rows() && solver->MtM.cols())
+		std::cout << "AtA:" << solver->MtM << std::endl;
+}
+
diff --git a/intern/eigen/intern/linear_solver.h b/intern/eigen/intern/linear_solver.h
new file mode 100644
index 00000000000..2dbea4d6f68
--- /dev/null
+++ b/intern/eigen/intern/linear_solver.h
@@ -0,0 +1,71 @@
+/*
+ *  Sparse linear solver.
+ *  Copyright (C) 2004 Bruno Levy
+ *  Copyright (C) 2005-2015 Blender Foundation
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ */
+
+#pragma once
+
+#include <stdbool.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Solvers for Ax = b and AtAx = Atb */
+
+typedef struct LinearSolver LinearSolver;
+
+LinearSolver *EIG_linear_solver_new(
+	int num_rows,
+	int num_columns,
+	int num_right_hand_sides);
+
+LinearSolver *EIG_linear_least_squares_solver_new(
+	int num_rows,
+	int num_columns,
+	int num_right_hand_sides);
+
+void EIG_linear_solver_delete(LinearSolver *solver);
+
+/* Variables (x). Any locking must be done before matrix construction. */
+
+void EIG_linear_solver_variable_set(LinearSolver *solver, int rhs, int index, double value);
+double EIG_linear_solver_variable_get(LinearSolver *solver, int rhs, int index);
+void EIG_linear_solver_variable_lock(LinearSolver *solver, int index);
+
+/* Matrix (A) and right hand side (b) */
+
+void EIG_linear_solver_matrix_add(LinearSolver *solver, int row, int col, double value);
+void EIG_linear_solver_right_hand_side_add(LinearSolver *solver, int rhs, int index, double value);
+
+/* Solve. Repeated solves are supported, by changing b between solves. */
+
+bool EIG_linear_solver_solve(LinearSolver *solver);
+
+/* Debugging */
+
+void EIG_linear_solver_print_matrix(LinearSolver *solver);
+
+#ifdef __cplusplus
+}
+#endif
+
diff --git a/intern/eigen/intern/svd.cc b/intern/eigen/intern/svd.cc
index e39a8261edb..4ed65af2a8f 100644
--- a/intern/eigen/intern/svd.cc
+++ b/intern/eigen/intern/svd.cc
@@ -48,7 +48,7 @@ using Eigen::MatrixXf;
 using Eigen::VectorXf;
 using Eigen::Map;
 
-void EG3_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V)
+void EIG_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V)
 {
 	/* Since our matrix is squared, we can use thinU/V. */
 	unsigned int flags = (r_U ? ComputeThinU : 0) | (r_V ? ComputeThinV : 0);
diff --git a/intern/eigen/intern/svd.h b/intern/eigen/intern/svd.h
index 0ac51108977..feadcc3520a 100644
--- a/intern/eigen/intern/svd.h
+++ b/intern/eigen/intern/svd.h
@@ -31,7 +31,7 @@
 extern "C" {
 #endif
 
-void EG3_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V);
+void EIG_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V);
 
 #ifdef __cplusplus
 }
diff --git a/intern/opennl/CMakeLists.txt b/intern/opennl/CMakeLists.txt
deleted file mode 100644
index 9416bc2b9ea..00000000000
--- a/intern/opennl/CMakeLists.txt
+++ /dev/null
@@ -1,58 +0,0 @@
-# ***** BEGIN GPL LICENSE BLOCK *****
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License
-# as published by the Free Software Foundation; either version 2
-# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-#
-# The Original Code is Copyright (C) 2006, Blender Foundation
-# All rights reserved.
-#
-# The Original Code is: all of this file.
-#
-# Contributor(s): Jacques Beaurain.
-#
-# ***** END GPL LICENSE BLOCK *****
-
-# External project, better not fix warnings.
-remove_strict_flags()
-
-# remove debug flag here since this is not a blender maintained library
-# and debug gives a lot of prints on UV unwrapping. developers can enable if they need to.
-if(MSVC)
-	remove_definitions(-DDEBUG)
-else()
-	add_definitions(-UDEBUG)
-endif()
-
-
-# quiet compiler warnings about undefined defines
-add_definitions(
-	-DDEBUGlevel=0
-	-DPRNTlevel=0
-)
-
-set(INC
-	extern
-)
-
-set(INC_SYS
-	../../extern/colamd/Include
-	../../extern/Eigen3
-)
-
-set(SRC
-	intern/opennl.cpp
-	extern/ONL_opennl.h
-)
-
-blender_add_lib(bf_intern_opennl "${SRC}" "${INC}" "${INC_SYS}")
diff --git a/intern/opennl/SConscript b/intern/opennl/SConscript
deleted file mode 100644
index 99df29b780a..00000000000
--- a/intern/opennl/SConscript
+++ /dev/null
@@ -1,35 +0,0 @@
-#!/usr/bin/env python
-#
-# ***** BEGIN GPL LICENSE BLOCK *****
-#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License
-# as published by the Free Software Foundation; either version 2
-# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-#
-# The Original Code is Copyright (C) 2006, Blender Foundation
-# All rights reserved.
-#
-# The Original Code is: all of this file.
-#
-# Contributor(s): Nathan Letwory.
-#
-# ***** END GPL LICENSE BLOCK *****
-
-Import ('env')
-
-sources = env.Glob('intern/*.cpp')
-
-incs = 'extern ../../extern/colamd/Include ../../extern/Eigen3'
-
-env.BlenderLib ('bf_intern_opennl', sources, Split(incs), [], libtype=['intern','player'], priority=[100,90] )
-
diff --git a/intern/opennl/extern/ONL_opennl.h b/intern/opennl/extern/ONL_opennl.h
deleted file mode 100644
index a414f40824b..00000000000
--- a/intern/opennl/extern/ONL_opennl.h
+++ /dev/null
@@ -1,113 +0,0 @@
-/** \file opennl/extern/ONL_opennl.h
- *  \ingroup opennlextern
- */
-/*
- *  OpenNL: Numerical Library
- *  Copyright (C) 2004 Bruno Levy
- *
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation; either version 2 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
- *
- *  If you modify this software, you should include a notice giving the
- *  name of the person performing the modification, the date of modification,
- *  and the reason for such modification.
- *
- *  Contact: Bruno Levy
- *
- *     levy@loria.fr
- *
- *     ISA Project
- *     LORIA, INRIA Lorraine, 
- *     Campus Scientifique, BP 239
- *     54506 VANDOEUVRE LES NANCY CEDEX 
- *     FRANCE
- *
- *  Note that the GNU General Public License does not permit incorporating
- *  the Software into proprietary programs. 
- */
-
-#ifndef nlOPENNL_H
-#define nlOPENNL_H
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/* Datatypes */
-
-typedef unsigned int	NLenum;
-typedef unsigned char	NLboolean;
-typedef int				NLint;		/* 4-byte signed */
-typedef unsigned int	NLuint;		/* 4-byte unsigned */
-typedef double			NLdouble;	/* double precision float */
-
-typedef struct NLContext NLContext;
-
-/* Constants */
-
-#define NL_FALSE   0x0
-#define NL_TRUE    0x1
-
-/* Primitives */
-
-#define NL_SYSTEM  0x0
-#define NL_MATRIX  0x1
-
-/* Solver Parameters */
-
-#define NL_SOLVER              0x100
-#define NL_NB_VARIABLES        0x101
-#define NL_LEAST_SQUARES       0x102
-#define NL_ERROR               0x108
-#define NL_NB_ROWS             0x110
-#define NL_NB_RIGHT_HAND_SIDES 0x112 /* 4 max */
-
-/* Contexts */
-
-NLContext *nlNewContext(void);
-void nlDeleteContext(NLContext *context);
-
-/* State get/set */
-
-void nlSolverParameteri(NLContext *context, NLenum pname, NLint param);
-
-/* Variables */
-
-void nlSetVariable(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value);
-NLdouble nlGetVariable(NLContext *context, NLuint rhsindex, NLuint index);
-void nlLockVariable(NLContext *context, NLuint index);
-void nlUnlockVariable(NLContext *context, NLuint index);
-
-/* Begin/End */
-
-void nlBegin(NLContext *context, NLenum primitive);
-void nlEnd(NLContext *context, NLenum primitive);
-
-/* Setting elements in matrix/vector */
-
-void nlMatrixAdd(NLContext *context, NLuint row, NLuint col, NLdouble value);
-void nlRightHandSideAdd(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value);
-void nlRightHandSideSet(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value);
-
-/* Solve */
-
-void nlPrintMatrix(NLContext *context);
-NLboolean nlSolve(NLContext *context, NLboolean solveAgain);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
-
diff --git a/intern/opennl/intern/opennl.cpp b/intern/opennl/intern/opennl.cpp
deleted file mode 100644
index 331291e790a..00000000000
--- a/intern/opennl/intern/opennl.cpp
+++ /dev/null
@@ -1,492 +0,0 @@
-/** \file opennl/intern/opennl.c
- *  \ingroup opennlintern
- */
-/*
- *
- *  OpenNL: Numerical Library
- *  Copyright (C) 2004 Bruno Levy
- *
- *  This program is free software; you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation; either version 2 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
- *
- *  If you modify this software, you should include a notice giving the
- *  name of the person performing the modification, the date of modification,
- *  and the reason for such modification.
- *
- *  Contact: Bruno Levy
- *
- *	 levy@loria.fr
- *
- *	 ISA Project
- *	 LORIA, INRIA Lorraine,
- *	 Campus Scientifique, BP 239
- *	 54506 VANDOEUVRE LES NANCY CEDEX
- *	 FRANCE
- *
- *  Note that the GNU General Public License does not permit incorporating
- *  the Software into proprietary programs.
- */
-
-#include "ONL_opennl.h"
-
-#include <Eigen/Sparse>
-
-#include <algorithm>
-#include <cassert>
-#include <cstdlib>
-#include <iostream>
-#include <vector>
-
-/* Eigen data structures */
-
-typedef Eigen::SparseMatrix<double, Eigen::ColMajor> EigenSparseMatrix;
-typedef Eigen::SparseLU<EigenSparseMatrix> EigenSparseSolver;
-typedef Eigen::VectorXd EigenVectorX;
-typedef Eigen::Triplet<double> EigenTriplet;
-
-/* NLContext data structure */
-
-struct NLCoeff {
-	NLCoeff()
-	{
-		index = 0;
-		value = 0.0;
-	}
-	NLuint index;
-	NLdouble value;
-};
-
-struct NLVariable {
-	NLVariable()
-	{
-		memset(value, 0, sizeof(value));
-		locked = false;
-		index = 0;
-	}
-	NLdouble value[4];
-	NLboolean locked;
-	NLuint index;
-	std::vector<NLCoeff> a;
-};
-
-#define NL_STATE_INITIAL            0
-#define NL_STATE_SYSTEM             1
-#define NL_STATE_MATRIX             2
-#define NL_STATE_MATRIX_CONSTRUCTED 3
-#define NL_STATE_SYSTEM_CONSTRUCTED 4
-#define NL_STATE_SYSTEM_SOLVED      5
-
-struct NLContext {
-   NLContext()
-   {
-      state = NL_STATE_INITIAL;
-      n = 0;
-      m = 0;
-      sparse_solver = NULL;
-      nb_variables = 0;
-      nb_rhs = 1;
-      nb_rows = 0;
-      least_squares = false;
-      solve_again = false;
-   }
-
-   ~NLContext()
-   {
-      delete sparse_solver;
-   }
-
-	NLenum state;
-
-	NLuint n;
-	NLuint m;
-
-	std::vector<EigenTriplet> Mtriplets;
-	EigenSparseMatrix M;
-	EigenSparseMatrix MtM;
-	std::vector<EigenVectorX> b;
-	std::vector<EigenVectorX> Mtb;
-	std::vector<EigenVectorX> x;
-
-	EigenSparseSolver *sparse_solver;
-
-	NLuint nb_variables;
-	std::vector<NLVariable> variable;
-
-	NLuint nb_rows;
-	NLuint nb_rhs;
-
-	NLboolean least_squares;
-	NLboolean solve_again;
-};
-
-NLContext *nlNewContext(void)
-{
-	return new NLContext();
-}
-
-void nlDeleteContext(NLContext *context)
-{
-	delete context;
-}
-
-static void __nlCheckState(NLContext *context, NLenum state)
-{
-	assert(context->state == state);
-}
-
-static void __nlTransition(NLContext *context, NLenum from_state, NLenum to_state)
-{
-	__nlCheckState(context, from_state);
-	context->state = to_state;
-}
-
-/* Get/Set parameters */
-
-void nlSolverParameteri(NLContext *context, NLenum pname, NLint param)
-{
-	__nlCheckState(context, NL_STATE_INITIAL);
-	switch(pname) {
-	case NL_NB_VARIABLES: {
-		assert(param > 0);
-		context->nb_variables = (NLuint)param;
-	} break;
-	case NL_NB_ROWS: {
-		assert(param > 0);
-		context->nb_rows = (NLuint)param;
-	} break;
-	case NL_LEAST_SQUARES: {
-		context->least_squares = (NLboolean)param;
-	} break;
-	case NL_NB_RIGHT_HAND_SIDES: {
-		context->nb_rhs = (NLuint)param;
-	} break;
-	default: {
-		assert(0);
-	} break;
-	}
-}
-
-/* Get/Set Lock/Unlock variables */
-
-void nlSetVariable(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value)
-{
-	__nlCheckState(context, NL_STATE_SYSTEM);
-	context->variable[index].value[rhsindex] = value;
-}
-
-NLdouble nlGetVariable(NLContext *context, NLuint rhsindex, NLuint index)
-{
-	assert(context->state != NL_STATE_INITIAL);
-	return context->variable[index].value[rhsindex];
-}
-
-void nlLockVariable(NLContext *context, NLuint index)
-{
-	__nlCheckState(context, NL_STATE_SYSTEM);
-	context->variable[index].locked = true;
-}
-
-void nlUnlockVariable(NLContext *context, NLuint index)
-{
-	__nlCheckState(context, NL_STATE_SYSTEM);
-	context->variable[index].locked = false;
-}
-
-/* System construction */
-
-static void __nlVariablesToVector(NLContext *context)
-{
-	NLuint i, j, nb_rhs;
-
-	nb_rhs= context->nb_rhs;
-
-	for(i=0; i<context->nb_variables; i++) {
-		NLVariable* v = &(context->variable[i]);
-		if(!v->locked) {
-			for(j=0; j<nb_rhs; j++)
-				context->x[j][v->index] = v->value[j];
-		}
-	}
-}
-
-static void __nlVectorToVariables(NLContext *context)
-{
-	NLuint i, j, nb_rhs;
-
-	nb_rhs= context->nb_rhs;
-
-	for(i=0; i<context->nb_variables; i++) {
-		NLVariable* v = &(context->variable[i]);
-		if(!v->locked) {
-			for(j=0; j<nb_rhs; j++)
-				v->value[j] = context->x[j][v->index];
-		}
-	}
-}
-
-static void __nlBeginSystem(NLContext *context)
-{
-	assert(context->nb_variables > 0);
-
-	if (context->solve_again)
-		__nlTransition(context, NL_STATE_SYSTEM_SOLVED, NL_STATE_SYSTEM);
-	else {
-		__nlTransition(context, NL_STATE_INITIAL, NL_STATE_SYSTEM);
-
-		context->variable.resize(context->nb_variables);
-	}
-}
-
-static void __nlEndSystem(NLContext *context)
-{
-	__nlTransition(context, NL_STATE_MATRIX_CONSTRUCTED, NL_STATE_SYSTEM_CONSTRUCTED);
-}
-
-static void __nlBeginMatrix(NLContext *context)
-{
-	NLuint i;
-	NLuint m = 0, n = 0;
-
-	__nlTransition(context, NL_STATE_SYSTEM, NL_STATE_MATRIX);
-
-	if (!context->solve_again) {
-		for(i=0; i<context->nb_variables; i++) {
-			if(context->variable[i].locked)
-				context->variable[i].index = ~0;
-			else
-				context->variable[i].index = n++;
-		}
-
-		m = (context->nb_rows == 0)? n: context->nb_rows;
-
-		context->m = m;
-		context->n = n;
-
-		/* reserve reasonable estimate */
-		context->Mtriplets.clear();
-		context->Mtriplets.reserve(std::max(m, n)*3);
-
-		context->b.resize(context->nb_rhs);
-		context->x.resize(context->nb_rhs);
-
-		for (i=0; i<context->nb_rhs; i++) {
-			context->b[i].setZero(m);
-			context->x[i].setZero(n);
-		}
-	}
-	else {
-		/* need to recompute b only, A is not constructed anymore */
-		for (i=0; i<context->nb_rhs; i++)
-			context->b[i].setZero(context->m);
-	}
-
-	__nlVariablesToVector(context);
-}
-
-static void __nlEndMatrixRHS(NLContext *context, NLuint rhs)
-{
-	NLVariable *variable;
-	NLuint i, j;
-
-	EigenVectorX& b = context->b[rhs];
-
-	for(i=0; i<context->nb_variables; i++) {
-		variable = &(context->variable[i]);
-
-		if(variable->locked) {
-			std::vector<NLCoeff>& a = variable->a;
-
-			for(j=0; j<a.size(); j++) {
-				b[a[j].index] -= a[j].value*variable->value[rhs];
-			}
-		}
-	}
-
-	if(context->least_squares)
-		context->Mtb[rhs] = context->M.transpose() * b;
-}
-
-static void __nlEndMatrix(NLContext *context)
-{
-	__nlTransition(context, NL_STATE_MATRIX, NL_STATE_MATRIX_CONSTRUCTED);
-
-	if(!context->solve_again) {
-		context->M.resize(context->m, context->n);
-		context->M.setFromTriplets(context->Mtriplets.begin(), context->Mtriplets.end());
-		context->Mtriplets.clear();
-
-		if(context->least_squares) {
-			context->MtM = context->M.transpose() * context->M;
-
-			context->Mtb.resize(context->nb_rhs);
-			for (NLuint rhs=0; rhs<context->nb_rhs; rhs++)
-				context->Mtb[rhs].setZero(context->n);
-		}
-	}
-
-	for (NLuint rhs=0; rhs<context->nb_rhs; rhs++)
-		__nlEndMatrixRHS(context, rhs);
-}
-
-void nlMatrixAdd(NLContext *context, NLuint row, NLuint col, NLdouble value)
-{
-	__nlCheckState(context, NL_STATE_MATRIX);
-
-	if(context->solve_again)
-		return;
-
-	if (!context->least_squares && context->variable[row].locked);
-	else if (context->variable[col].locked) {
-		if(!context->least_squares)
-			row = context->variable[row].index;
-
-		NLCoeff coeff;
-		coeff.index = row;
-		coeff.value = value;
-		context->variable[col].a.push_back(coeff);
-	}
-	else {
-		if(!context->least_squares)
-			row = context->variable[row].index;
-		col = context->variable[col].index;
-
-		// direct insert into matrix is too slow, so use triplets
-		EigenTriplet triplet(row, col, value);
-		context->Mtriplets.push_back(triplet);
-	}
-}
-
-void nlRightHandSideAdd(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value)
-{
-	__nlCheckState(context, NL_STATE_MATRIX);
-
-	if(context->least_squares) {
-		context->b[rhsindex][index] += value;
-	}
-	else {
-		if(!context->variable[index].locked) {
-			index = context->variable[index].index;
-			context->b[rhsindex][index] += value;
-		}
-	}
-}
-
-void nlRightHandSideSet(NLContext *context, NLuint rhsindex, NLuint index, NLdouble value)
-{
-	__nlCheckState(context, NL_STATE_MATRIX);
-
-	if(context->least_squares) {
-		context->b[rhsindex][index] = value;
-	}
-	else {
-		if(!context->variable[index].locked) {
-			index = context->variable[index].index;
-			context->b[rhsindex][index] = value;
-		}
-	}
-}
-
-void nlBegin(NLContext *context, NLenum prim)
-{
-	switch(prim) {
-	case NL_SYSTEM: {
-		__nlBeginSystem(context);
-	} break;
-	case NL_MATRIX: {
-		__nlBeginMatrix(context);
-	} break;
-	default: {
-		assert(0);
-	}
-	}
-}
-
-void nlEnd(NLContext *context, NLenum prim)
-{
-	switch(prim) {
-	case NL_SYSTEM: {
-		__nlEndSystem(context);
-	} break;
-	case NL_MATRIX: {
-		__nlEndMatrix(context);
-	} break;
-	default: {
-		assert(0);
-	}
-	}
-}
-
-void nlPrintMatrix(NLContext *context)
-{
-	std::cout << "A:" << context->M << std::endl;
-
-	for(NLuint rhs=0; rhs<context->nb_rhs; rhs++)
-		std::cout << "b " << rhs << ":" << context->b[rhs] << std::endl;
-
-	if (context->MtM.rows() && context->MtM.cols())
-		std::cout << "AtA:" << context->MtM << std::endl;
-}
-
-/* Solving */
-
-NLboolean nlSolve(NLContext *context, NLboolean solveAgain)
-{
-	NLboolean result = true;
-
-	__nlCheckState(context, NL_STATE_SYSTEM_CONSTRUCTED);
-
-	if (!context->solve_again) {
-		EigenSparseMatrix& M = (context->least_squares)? context->MtM: context->M;
-
-		assert(M.rows() == M.cols());
-
-		/* Convert M to compressed column format */
-		M.makeCompressed();
-
-		/* Perform sparse LU factorization */
-		EigenSparseSolver *sparse_solver = new EigenSparseSolver();
-		context->sparse_solver = sparse_solver;
-
-		sparse_solver->analyzePattern(M);
-		sparse_solver->factorize(M);
-
-		result = (sparse_solver->info() == Eigen::Success);
-
-		/* Free M, don't need it anymore at this point */
-		M.resize(0, 0);
-	}
-
-	if (result) {
-		/* Solve each right hand side */
-		for(NLuint rhs=0; rhs<context->nb_rhs; rhs++) {
-			EigenVectorX& b = (context->least_squares)? context->Mtb[rhs]: context->b[rhs];
-			context->x[rhs] = context->sparse_solver->solve(b);
-
-			if (context->sparse_solver->info() != Eigen::Success)
-				result = false;
-		}
-
-		if (result) {
-			__nlVectorToVariables(context);
-
-			if (solveAgain)
-				context->solve_again = true;
-
-			__nlTransition(context, NL_STATE_SYSTEM_CONSTRUCTED, NL_STATE_SYSTEM_SOLVED);
-		}
-	}
-
-	return result;
-}
-