diff options
author | Brecht Van Lommel <brechtvanlommel@gmail.com> | 2015-11-24 22:42:10 +0300 |
---|---|---|
committer | Brecht Van Lommel <brechtvanlommel@gmail.com> | 2015-12-10 03:58:10 +0300 |
commit | f9047c3f8c72f1a15a4c051b507306d308f44646 (patch) | |
tree | 5116007ad245bbbc95695f13afa7134983e42be4 /intern | |
parent | 858b680a50888a071d5d37af261b0c89b47aea8c (diff) |
Eigen: fold remaining OpenNL code into intern/eigen.
Differential Revision: https://developer.blender.org/D1662
Diffstat (limited to 'intern')
-rw-r--r-- | intern/CMakeLists.txt | 4 | ||||
-rw-r--r-- | intern/SConscript | 1 | ||||
-rw-r--r-- | intern/eigen/CMakeLists.txt | 2 | ||||
-rw-r--r-- | intern/eigen/eigen_capi.h | 1 | ||||
-rw-r--r-- | intern/eigen/intern/eigenvalues.cc | 2 | ||||
-rw-r--r-- | intern/eigen/intern/eigenvalues.h | 2 | ||||
-rw-r--r-- | intern/eigen/intern/linear_solver.cc | 354 | ||||
-rw-r--r-- | intern/eigen/intern/linear_solver.h | 71 | ||||
-rw-r--r-- | intern/eigen/intern/svd.cc | 2 | ||||
-rw-r--r-- | intern/eigen/intern/svd.h | 2 | ||||
-rw-r--r-- | intern/opennl/CMakeLists.txt | 58 | ||||
-rw-r--r-- | intern/opennl/SConscript | 35 | ||||
-rw-r--r-- | intern/opennl/extern/ONL_opennl.h | 113 | ||||
-rw-r--r-- | intern/opennl/intern/opennl.cpp | 492 |
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; -} - |