diff options
| -rw-r--r-- | intern/opennl/extern/ONL_opennl.h | 21 | ||||
| -rw-r--r-- | intern/opennl/intern/opennl.cpp | 960 | ||||
| -rw-r--r-- | source/blender/bmesh/operators/bmo_smooth_laplacian.c | 2 | ||||
| -rw-r--r-- | source/blender/editors/armature/meshlaplacian.c | 6 | ||||
| -rw-r--r-- | source/blender/editors/armature/reeb.c | 2 | ||||
| -rw-r--r-- | source/blender/editors/uvedit/uvedit_parametrizer.c | 6 | ||||
| -rw-r--r-- | source/blender/modifiers/intern/MOD_laplaciandeform.c | 9 | ||||
| -rw-r--r-- | source/blender/modifiers/intern/MOD_laplaciansmooth.c | 2 |
8 files changed, 222 insertions, 786 deletions
diff --git a/intern/opennl/extern/ONL_opennl.h b/intern/opennl/extern/ONL_opennl.h index d550a638693..61a11fa8993 100644 --- a/intern/opennl/extern/ONL_opennl.h +++ b/intern/opennl/extern/ONL_opennl.h @@ -37,11 +37,6 @@ * the Software into proprietary programs. */ -/* -#define NL_DEBUG -#define NL_PARANOID -*/ - #ifndef nlOPENNL_H #define nlOPENNL_H @@ -49,8 +44,6 @@ extern "C" { #endif -#define NL_VERSION_0_0 1 - /* Datatypes */ typedef unsigned int NLenum; @@ -59,7 +52,7 @@ typedef int NLint; /* 4-byte signed */ typedef unsigned int NLuint; /* 4-byte unsigned */ typedef double NLdouble; /* double precision float */ -typedef void* NLContext; +typedef struct NLContext NLContext; /* Constants */ @@ -76,17 +69,16 @@ typedef void* NLContext; #define NL_SOLVER 0x100 #define NL_NB_VARIABLES 0x101 #define NL_LEAST_SQUARES 0x102 -#define NL_SYMMETRIC 0x106 #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); -void nlMakeCurrent(NLContext context); -NLContext nlGetCurrent(void); +NLContext *nlNewContext(void); +void nlDeleteContext(NLContext *context); +void nlMakeCurrent(NLContext *context); +NLContext *nlGetCurrent(void); /* State get/set */ @@ -113,8 +105,7 @@ void nlRightHandSideSet(NLuint rhsindex, NLuint index, NLdouble value); /* Solve */ void nlPrintMatrix(void); -NLboolean nlSolve(void); -NLboolean nlSolveAdvanced(NLint *permutation, NLboolean solveAgain); +NLboolean nlSolve(NLboolean solveAgain); #ifdef __cplusplus } diff --git a/intern/opennl/intern/opennl.cpp b/intern/opennl/intern/opennl.cpp index 085375b0d8c..de128053342 100644 --- a/intern/opennl/intern/opennl.cpp +++ b/intern/opennl/intern/opennl.cpp @@ -40,688 +40,232 @@ #include "ONL_opennl.h" -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <math.h> - -#ifdef NL_PARANOID -#ifndef NL_DEBUG -#define NL_DEBUG -#endif -#endif - #include <Eigen/Sparse> + +#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; -/************************************************************************************/ -/* Assertions */ - - -static void __nl_assertion_failed(const char* cond, const char* file, int line) { - fprintf( - stderr, - "OpenNL assertion failed: %s, file:%s, line:%d\n", - cond,file,line - ); - abort(); -} - -static void __nl_range_assertion_failed( - double x, double min_val, double max_val, const char* file, int line -) { - fprintf( - stderr, - "OpenNL range assertion failed: %f in [ %f ... %f ], file:%s, line:%d\n", - x, min_val, max_val, file,line - ); - abort(); -} - -static void __nl_should_not_have_reached(const char* file, int line) { - fprintf( - stderr, - "OpenNL should not have reached this point: file:%s, line:%d\n", - file,line - ); - abort(); -} - - -#define __nl_assert(x) { \ - if(!(x)) { \ - __nl_assertion_failed(#x,__FILE__, __LINE__); \ - } \ -} - -#define __nl_range_assert(x,max_val) { \ - if(((x) > (max_val))) { \ - __nl_range_assertion_failed(x, 0.0, max_val, \ - __FILE__, __LINE__ \ - ); \ - } \ -} - -#define __nl_assert_not_reached { \ - __nl_should_not_have_reached(__FILE__, __LINE__); \ -} - -#ifdef NL_DEBUG -#define __nl_debug_assert(x) __nl_assert(x) -#define __nl_debug_range_assert(x,min_val,max_val) __nl_range_assert(x,min_val,max_val) -#else -#define __nl_debug_assert(x) -#define __nl_debug_range_assert(x,min_val,max_val) -#endif - -#ifdef NL_PARANOID -#define __nl_parano_assert(x) __nl_assert(x) -#define __nl_parano_range_assert(x,min_val,max_val) __nl_range_assert(x,min_val,max_val) -#else -#define __nl_parano_assert(x) -#define __nl_parano_range_assert(x,min_val,max_val) -#endif - -/************************************************************************************/ -/* classic macros */ - -#ifndef MIN -#define MIN(x,y) (((x) < (y)) ? (x) : (y)) -#endif - -#ifndef MAX -#define MAX(x,y) (((x) > (y)) ? (x) : (y)) -#endif - -/************************************************************************************/ -/* memory management */ - -#define __NL_NEW(T) (T*)(calloc(1, sizeof(T))) -#define __NL_NEW_ARRAY(T,NB) (T*)(calloc(MAX(NB, 1),sizeof(T))) -#define __NL_RENEW_ARRAY(T,x,NB) (T*)(realloc(x,(NB)*sizeof(T))) -#define __NL_DELETE(x) if(x) free(x); x = NULL -#define __NL_DELETE_ARRAY(x) if(x) free(x); x = NULL - -#define __NL_CLEAR(T, x) memset(x, 0, sizeof(T)) -#define __NL_CLEAR_ARRAY(T,x,NB) if(NB) memset(x, 0, (NB)*sizeof(T)) - -/************************************************************************************/ -/* Dynamic arrays for sparse row/columns */ +/* NLContext data structure */ typedef struct { - NLuint index; + NLuint index; NLdouble value; -} __NLCoeff; +} NLCoeff; typedef struct { - NLuint size; - NLuint capacity; - __NLCoeff* coeff; -} __NLRowColumn; - -static void __nlRowColumnConstruct(__NLRowColumn* c) { - c->size = 0; - c->capacity = 0; - c->coeff = NULL; -} - -static void __nlRowColumnDestroy(__NLRowColumn* c) { - __NL_DELETE_ARRAY(c->coeff); -#ifdef NL_PARANOID - __NL_CLEAR(__NLRowColumn, c); -#endif -} - -static void __nlRowColumnGrow(__NLRowColumn* c) { - if(c->capacity != 0) { - c->capacity = 2 * c->capacity; - c->coeff = __NL_RENEW_ARRAY(__NLCoeff, c->coeff, c->capacity); - } else { - c->capacity = 4; - c->coeff = __NL_NEW_ARRAY(__NLCoeff, c->capacity); - } -} - -static void __nlRowColumnAdd(__NLRowColumn* c, NLint index, NLdouble value) { - NLuint i; - for(i=0; i<c->size; i++) { - if(c->coeff[i].index == (NLuint)index) { - c->coeff[i].value += value; - return; - } - } - if(c->size == c->capacity) { - __nlRowColumnGrow(c); - } - c->coeff[c->size].index = index; - c->coeff[c->size].value = value; - c->size++; -} - -/* Does not check whether the index already exists */ -static void __nlRowColumnAppend(__NLRowColumn* c, NLint index, NLdouble value) { - if(c->size == c->capacity) { - __nlRowColumnGrow(c); - } - c->coeff[c->size].index = index; - c->coeff[c->size].value = value; - c->size++; -} - -static void __nlRowColumnClear(__NLRowColumn* c) { - c->size = 0; - c->capacity = 0; - __NL_DELETE_ARRAY(c->coeff); -} + NLdouble value[4]; + NLboolean locked; + NLuint index; + std::vector<NLCoeff> a; +} NLVariable; -/************************************************************************************/ -/* SparseMatrix data structure */ +#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 -#define __NL_ROWS 1 -#define __NL_COLUMNS 2 -#define __NL_SYMMETRIC 4 +struct NLContext { + NLenum state; -typedef struct { - NLuint m; NLuint n; - NLuint diag_size; - NLenum storage; - __NLRowColumn* row; - __NLRowColumn* column; - NLdouble* diag; -} __NLSparseMatrix; - - -static void __nlSparseMatrixConstruct( - __NLSparseMatrix* M, NLuint m, NLuint n, NLenum storage -) { - NLuint i; - M->m = m; - M->n = n; - M->storage = storage; - if(storage & __NL_ROWS) { - M->row = __NL_NEW_ARRAY(__NLRowColumn, m); - for(i=0; i<m; i++) { - __nlRowColumnConstruct(&(M->row[i])); - } - } else { - M->row = NULL; - } - - if(storage & __NL_COLUMNS) { - M->column = __NL_NEW_ARRAY(__NLRowColumn, n); - for(i=0; i<n; i++) { - __nlRowColumnConstruct(&(M->column[i])); - } - } else { - M->column = NULL; - } - - M->diag_size = MIN(m,n); - M->diag = __NL_NEW_ARRAY(NLdouble, M->diag_size); -} - -static void __nlSparseMatrixDestroy(__NLSparseMatrix* M) { - NLuint i; - __NL_DELETE_ARRAY(M->diag); - if(M->storage & __NL_ROWS) { - for(i=0; i<M->m; i++) { - __nlRowColumnDestroy(&(M->row[i])); - } - __NL_DELETE_ARRAY(M->row); - } - if(M->storage & __NL_COLUMNS) { - for(i=0; i<M->n; i++) { - __nlRowColumnDestroy(&(M->column[i])); - } - __NL_DELETE_ARRAY(M->column); - } -#ifdef NL_PARANOID - __NL_CLEAR(__NLSparseMatrix,M); -#endif -} - -static void __nlSparseMatrixAdd( - __NLSparseMatrix* M, NLuint i, NLuint j, NLdouble value -) { - __nl_parano_range_assert(i, 0, M->m - 1); - __nl_parano_range_assert(j, 0, M->n - 1); - if((M->storage & __NL_SYMMETRIC) && (j > i)) { - return; - } - if(i == j) { - M->diag[i] += value; - } - if(M->storage & __NL_ROWS) { - __nlRowColumnAdd(&(M->row[i]), j, value); - } - if(M->storage & __NL_COLUMNS) { - __nlRowColumnAdd(&(M->column[j]), i, value); - } -} - -static void __nlSparseMatrixClear( __NLSparseMatrix* M) { - NLuint i; - if(M->storage & __NL_ROWS) { - for(i=0; i<M->m; i++) { - __nlRowColumnClear(&(M->row[i])); - } - } - if(M->storage & __NL_COLUMNS) { - for(i=0; i<M->n; i++) { - __nlRowColumnClear(&(M->column[i])); - } - } - __NL_CLEAR_ARRAY(NLdouble, M->diag, M->diag_size); -} - -/************************************************************************************/ -/* SparseMatrix x Vector routines, internal helper routines */ - -static void __nlSparseMatrix_mult_rows_symmetric( - __NLSparseMatrix* A, NLdouble* x, NLdouble* y -) { - NLuint m = A->m; - NLuint i,ij; - __NLRowColumn* Ri = NULL; - __NLCoeff* c = NULL; - for(i=0; i<m; i++) { - y[i] = 0; - Ri = &(A->row[i]); - for(ij=0; ij<Ri->size; ij++) { - c = &(Ri->coeff[ij]); - y[i] += c->value * x[c->index]; - if(i != c->index) { - y[c->index] += c->value * x[i]; - } - } - } -} - -static void __nlSparseMatrix_mult_rows( - __NLSparseMatrix* A, NLdouble* x, NLdouble* y -) { - NLuint m = A->m; - NLuint i,ij; - __NLRowColumn* Ri = NULL; - __NLCoeff* c = NULL; - for(i=0; i<m; i++) { - y[i] = 0; - Ri = &(A->row[i]); - for(ij=0; ij<Ri->size; ij++) { - c = &(Ri->coeff[ij]); - y[i] += c->value * x[c->index]; - } - } -} - -static void __nlSparseMatrix_mult_cols_symmetric( - __NLSparseMatrix* A, NLdouble* x, NLdouble* y -) { - NLuint n = A->n; - NLuint j,ii; - __NLRowColumn* Cj = NULL; - __NLCoeff* c = NULL; - for(j=0; j<n; j++) { - y[j] = 0; - Cj = &(A->column[j]); - for(ii=0; ii<Cj->size; ii++) { - c = &(Cj->coeff[ii]); - y[c->index] += c->value * x[j]; - if(j != c->index) { - y[j] += c->value * x[c->index]; - } - } - } -} - -static void __nlSparseMatrix_mult_cols( - __NLSparseMatrix* A, NLdouble* x, NLdouble* y -) { - NLuint n = A->n; - NLuint j,ii; - __NLRowColumn* Cj = NULL; - __NLCoeff* c = NULL; - __NL_CLEAR_ARRAY(NLdouble, y, A->m); - for(j=0; j<n; j++) { - Cj = &(A->column[j]); - for(ii=0; ii<Cj->size; ii++) { - c = &(Cj->coeff[ii]); - y[c->index] += c->value * x[j]; - } - } -} - -/************************************************************************************/ -/* SparseMatrix x Vector routines, main driver routine */ - -static void __nlSparseMatrixMult(__NLSparseMatrix* A, NLdouble* x, NLdouble* y) { - if(A->storage & __NL_ROWS) { - if(A->storage & __NL_SYMMETRIC) { - __nlSparseMatrix_mult_rows_symmetric(A, x, y); - } else { - __nlSparseMatrix_mult_rows(A, x, y); - } - } else { - if(A->storage & __NL_SYMMETRIC) { - __nlSparseMatrix_mult_cols_symmetric(A, x, y); - } else { - __nlSparseMatrix_mult_cols(A, x, y); - } - } -} - -/* ****************** Routines for least squares ******************* */ - -static void __nlSparseMatrix_square( - __NLSparseMatrix* AtA, __NLSparseMatrix *A -) { - NLuint m = A->m; - NLuint n = A->n; - NLuint i, j0, j1; - __NLRowColumn *Ri = NULL; - __NLCoeff *c0 = NULL, *c1 = NULL; - double value; - - __nlSparseMatrixConstruct(AtA, n, n, A->storage); - - for(i=0; i<m; i++) { - Ri = &(A->row[i]); - - for(j0=0; j0<Ri->size; j0++) { - c0 = &(Ri->coeff[j0]); - for(j1=0; j1<Ri->size; j1++) { - c1 = &(Ri->coeff[j1]); - - value = c0->value*c1->value; - __nlSparseMatrixAdd(AtA, c0->index, c1->index, value); - } - } - } -} - -static void __nlSparseMatrix_transpose_mult_rows( - __NLSparseMatrix* A, NLdouble* x, NLdouble* y -) { - NLuint m = A->m; - NLuint n = A->n; - NLuint i,ij; - __NLRowColumn* Ri = NULL; - __NLCoeff* c = NULL; - - __NL_CLEAR_ARRAY(NLdouble, y, n); - - for(i=0; i<m; i++) { - Ri = &(A->row[i]); - for(ij=0; ij<Ri->size; ij++) { - c = &(Ri->coeff[ij]); - y[c->index] += c->value * x[i]; - } - } -} + NLuint m; -/************************************************************************************/ -/* NLContext data structure */ + std::vector<EigenTriplet> Mtriplets; + EigenSparseMatrix M; + EigenSparseMatrix MtM; + std::vector<EigenVectorX> b; + std::vector<EigenVectorX> Mtb; + std::vector<EigenVectorX> x; -typedef void(*__NLMatrixFunc)(double* x, double* y); + EigenSparseSolver *sparse_solver; -typedef struct { - NLdouble value[4]; - NLboolean locked; - NLuint index; - __NLRowColumn *a; -} __NLVariable; + NLuint nb_variables; + std::vector<NLVariable> variable; -#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 + NLuint nb_rows; + NLuint nb_rhs; -typedef struct { - NLenum state; - NLuint n; - NLuint m; - __NLVariable* variable; - NLdouble* b; - NLdouble* Mtb; - __NLSparseMatrix M; - __NLSparseMatrix MtM; - NLdouble* x; - NLuint nb_variables; - NLuint nb_rows; - NLboolean least_squares; - NLboolean symmetric; - NLuint nb_rhs; - NLboolean solve_again; - NLboolean alloc_M; - NLboolean alloc_MtM; - NLboolean alloc_variable; - NLboolean alloc_x; - NLboolean alloc_b; - NLboolean alloc_Mtb; - NLdouble error; - __NLMatrixFunc matrix_vector_prod; - - struct __NLEigenContext { - EigenSparseSolver *sparse_solver; - } eigen; -} __NLContext; - -static __NLContext* __nlCurrentContext = NULL; - -static void __nlMatrixVectorProd_default(NLdouble* x, NLdouble* y) { - __nlSparseMatrixMult(&(__nlCurrentContext->M), x, y); -} + NLboolean least_squares; + NLboolean solve_again; +}; +static NLContext* __nlCurrentContext = NULL; -NLContext nlNewContext(void) { - __NLContext* result = __NL_NEW(__NLContext); - result->state = __NL_STATE_INITIAL; - result->matrix_vector_prod = __nlMatrixVectorProd_default; +NLContext *nlNewContext(void) +{ + NLContext* result = new NLContext(); + result->state = NL_STATE_INITIAL; result->nb_rhs = 1; nlMakeCurrent(result); return result; } -static void __nlFree_EIGEN(__NLContext *context); - -void nlDeleteContext(NLContext context_in) { - __NLContext* context = (__NLContext*)(context_in); - int i; - - if(__nlCurrentContext == context) { +void nlDeleteContext(NLContext *context) +{ + if(__nlCurrentContext == context) __nlCurrentContext = NULL; - } - if(context->alloc_M) { - __nlSparseMatrixDestroy(&context->M); - } - if(context->alloc_MtM) { - __nlSparseMatrixDestroy(&context->MtM); - } - if(context->alloc_variable) { - for(i=0; i<context->nb_variables; i++) { - if(context->variable[i].a) { - __nlRowColumnDestroy(context->variable[i].a); - __NL_DELETE(context->variable[i].a); - } - } - __NL_DELETE_ARRAY(context->variable); - } - if(context->alloc_b) { - __NL_DELETE_ARRAY(context->b); - } - if(context->alloc_Mtb) { - __NL_DELETE_ARRAY(context->Mtb); - } - if(context->alloc_x) { - __NL_DELETE_ARRAY(context->x); - } - if (context->eigen.sparse_solver) { - __nlFree_EIGEN(context); - } + context->M.resize(0, 0); + context->MtM.resize(0, 0); + context->b.clear(); + context->Mtb.clear(); + context->x.clear(); + context->variable.clear(); -#ifdef NL_PARANOID - __NL_CLEAR(__NLContext, context); -#endif - __NL_DELETE(context); + delete context->sparse_solver; + context->sparse_solver = NULL; + + delete context; } -void nlMakeCurrent(NLContext context) { - __nlCurrentContext = (__NLContext*)(context); +void nlMakeCurrent(NLContext *context) +{ + __nlCurrentContext = context; } -NLContext nlGetCurrent(void) { +NLContext *nlGetCurrent(void) +{ return __nlCurrentContext; } -static void __nlCheckState(NLenum state) { - __nl_assert(__nlCurrentContext->state == state); +static void __nlCheckState(NLenum state) +{ + assert(__nlCurrentContext->state == state); } -static void __nlTransition(NLenum from_state, NLenum to_state) { +static void __nlTransition(NLenum from_state, NLenum to_state) +{ __nlCheckState(from_state); __nlCurrentContext->state = to_state; } -/************************************************************************************/ /* Get/Set parameters */ -void nlSolverParameteri(NLenum pname, NLint param) { - __nlCheckState(__NL_STATE_INITIAL); +void nlSolverParameteri(NLenum pname, NLint param) +{ + __nlCheckState(NL_STATE_INITIAL); switch(pname) { case NL_NB_VARIABLES: { - __nl_assert(param > 0); + assert(param > 0); __nlCurrentContext->nb_variables = (NLuint)param; } break; case NL_NB_ROWS: { - __nl_assert(param > 0); + assert(param > 0); __nlCurrentContext->nb_rows = (NLuint)param; } break; case NL_LEAST_SQUARES: { __nlCurrentContext->least_squares = (NLboolean)param; } break; - case NL_SYMMETRIC: { - __nlCurrentContext->symmetric = (NLboolean)param; - } break; case NL_NB_RIGHT_HAND_SIDES: { __nlCurrentContext->nb_rhs = (NLuint)param; } break; default: { - __nl_assert_not_reached; + assert(0); } break; } } -/************************************************************************************/ /* Get/Set Lock/Unlock variables */ -void nlSetVariable(NLuint rhsindex, NLuint index, NLdouble value) { - __nlCheckState(__NL_STATE_SYSTEM); - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); +void nlSetVariable(NLuint rhsindex, NLuint index, NLdouble value) +{ + __nlCheckState(NL_STATE_SYSTEM); __nlCurrentContext->variable[index].value[rhsindex] = value; } -NLdouble nlGetVariable(NLuint rhsindex, NLuint index) { - __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL); - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); +NLdouble nlGetVariable(NLuint rhsindex, NLuint index) +{ + assert(__nlCurrentContext->state != NL_STATE_INITIAL); return __nlCurrentContext->variable[index].value[rhsindex]; } -void nlLockVariable(NLuint index) { - __nlCheckState(__NL_STATE_SYSTEM); - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); - __nlCurrentContext->variable[index].locked = NL_TRUE; +void nlLockVariable(NLuint index) +{ + __nlCheckState(NL_STATE_SYSTEM); + __nlCurrentContext->variable[index].locked = true; } -void nlUnlockVariable(NLuint index) { - __nlCheckState(__NL_STATE_SYSTEM); - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); - __nlCurrentContext->variable[index].locked = NL_FALSE; +void nlUnlockVariable(NLuint index) +{ + __nlCheckState(NL_STATE_SYSTEM); + __nlCurrentContext->variable[index].locked = false; } -/************************************************************************************/ /* System construction */ -static void __nlVariablesToVector() { - __NLContext *context = __nlCurrentContext; +static void __nlVariablesToVector() +{ + NLContext *context = __nlCurrentContext; NLuint i, j, nb_rhs; - __nl_assert(context->alloc_x); - __nl_assert(context->alloc_variable); - nb_rhs= context->nb_rhs; for(i=0; i<context->nb_variables; i++) { - __NLVariable* v = &(context->variable[i]); + NLVariable* v = &(context->variable[i]); if(!v->locked) { - __nl_assert(v->index < context->n); - for(j=0; j<nb_rhs; j++) - context->x[context->n*j + v->index] = v->value[j]; + context->x[j][v->index] = v->value[j]; } } } -static void __nlVectorToVariables() { - __NLContext *context = __nlCurrentContext; +static void __nlVectorToVariables() +{ + NLContext *context = __nlCurrentContext; NLuint i, j, nb_rhs; - __nl_assert(context->alloc_x); - __nl_assert(context->alloc_variable); - nb_rhs= context->nb_rhs; for(i=0; i<context->nb_variables; i++) { - __NLVariable* v = &(context->variable[i]); + NLVariable* v = &(context->variable[i]); if(!v->locked) { - __nl_assert(v->index < context->n); - for(j=0; j<nb_rhs; j++) - v->value[j] = context->x[context->n*j + v->index]; + v->value[j] = context->x[j][v->index]; } } } -static void __nlBeginSystem() { - __nl_assert(__nlCurrentContext->nb_variables > 0); +static void __nlBeginSystem() +{ + assert(__nlCurrentContext->nb_variables > 0); if (__nlCurrentContext->solve_again) - __nlTransition(__NL_STATE_SYSTEM_SOLVED, __NL_STATE_SYSTEM); + __nlTransition(NL_STATE_SYSTEM_SOLVED, NL_STATE_SYSTEM); else { - __nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM); + __nlTransition(NL_STATE_INITIAL, NL_STATE_SYSTEM); - __nlCurrentContext->variable = __NL_NEW_ARRAY( - __NLVariable, __nlCurrentContext->nb_variables); - - __nlCurrentContext->alloc_variable = NL_TRUE; + __nlCurrentContext->variable.resize(__nlCurrentContext->nb_variables); } } -static void __nlEndSystem() { - __nlTransition(__NL_STATE_MATRIX_CONSTRUCTED, __NL_STATE_SYSTEM_CONSTRUCTED); +static void __nlEndSystem() +{ + __nlTransition(NL_STATE_MATRIX_CONSTRUCTED, NL_STATE_SYSTEM_CONSTRUCTED); } -static void __nlBeginMatrix() { +static void __nlBeginMatrix() +{ NLuint i; NLuint m = 0, n = 0; - NLenum storage = __NL_ROWS; - __NLContext *context = __nlCurrentContext; + NLContext *context = __nlCurrentContext; - __nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX); + __nlTransition(NL_STATE_SYSTEM, NL_STATE_MATRIX); if (!context->solve_again) { for(i=0; i<context->nb_variables; i++) { - if(context->variable[i].locked) { + if(context->variable[i].locked) context->variable[i].index = ~0; - context->variable[i].a = __NL_NEW(__NLRowColumn); - __nlRowColumnConstruct(context->variable[i].a); - } else context->variable[i].index = n++; } @@ -731,75 +275,76 @@ static void __nlBeginMatrix() { context->m = m; context->n = n; - __nlSparseMatrixConstruct(&context->M, m, n, storage); - context->alloc_M = NL_TRUE; - - context->b = __NL_NEW_ARRAY(NLdouble, m*context->nb_rhs); - context->alloc_b = NL_TRUE; + context->b.resize(context->nb_rhs); + context->x.resize(context->nb_rhs); - context->x = __NL_NEW_ARRAY(NLdouble, n*context->nb_rhs); - context->alloc_x = NL_TRUE; + 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 */ - __NL_CLEAR_ARRAY(NLdouble, context->b, context->m*context->nb_rhs); + for (i=0; i<context->nb_rhs; i++) + context->b[i].setZero(context->m); } __nlVariablesToVector(); } -static void __nlEndMatrixRHS(NLuint rhs) { - __NLContext *context = __nlCurrentContext; - __NLVariable *variable; - __NLRowColumn *a; - NLdouble *b, *Mtb; +static void __nlEndMatrixRHS(NLuint rhs) +{ + NLContext *context = __nlCurrentContext; + NLVariable *variable; NLuint i, j; - b = context->b + context->m*rhs; - Mtb = context->Mtb + context->n*rhs; + EigenVectorX& b = context->b[rhs]; for(i=0; i<__nlCurrentContext->nb_variables; i++) { variable = &(context->variable[i]); if(variable->locked) { - a = variable->a; + std::vector<NLCoeff>& a = variable->a; - for(j=0; j<a->size; j++) { - b[a->coeff[j].index] -= a->coeff[j].value*variable->value[rhs]; + for(j=0; j<a.size(); j++) { + b[a[j].index] -= a[j].value*variable->value[rhs]; } } } if(context->least_squares) - __nlSparseMatrix_transpose_mult_rows(&context->M, b, Mtb); + context->Mtb[rhs] = context->M.transpose() * b; } -static void __nlEndMatrix() { - __NLContext *context = __nlCurrentContext; - NLuint i; +static void __nlEndMatrix() +{ + NLContext *context = __nlCurrentContext; - __nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED); + __nlTransition(NL_STATE_MATRIX, NL_STATE_MATRIX_CONSTRUCTED); - if(context->least_squares) { - if(!__nlCurrentContext->solve_again) { - __nlSparseMatrix_square(&context->MtM, &context->M); - context->alloc_MtM = NL_TRUE; + 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 = - __NL_NEW_ARRAY(NLdouble, context->n*context->nb_rhs); - context->alloc_Mtb = NL_TRUE; + context->Mtb.resize(context->nb_rhs); + for (NLuint rhs=0; rhs<context->nb_rhs; rhs++) + context->Mtb[rhs].setZero(context->n); } } - for(i=0; i<context->nb_rhs; i++) - __nlEndMatrixRHS(i); + for (NLuint rhs=0; rhs<context->nb_rhs; rhs++) + __nlEndMatrixRHS(rhs); } void nlMatrixAdd(NLuint row, NLuint col, NLdouble value) { - __NLContext *context = __nlCurrentContext; + NLContext *context = __nlCurrentContext; - __nlCheckState(__NL_STATE_MATRIX); + __nlCheckState(NL_STATE_MATRIX); if(context->solve_again) return; @@ -808,65 +353,57 @@ void nlMatrixAdd(NLuint row, NLuint col, NLdouble value) else if (context->variable[col].locked) { if(!context->least_squares) row = context->variable[row].index; - __nlRowColumnAppend(context->variable[col].a, row, value); + + NLCoeff coeff = {row, value}; + context->variable[col].a.push_back(coeff); } else { - __NLSparseMatrix* M = &context->M; - if(!context->least_squares) row = context->variable[row].index; col = context->variable[col].index; - - __nl_range_assert(row, context->m - 1); - __nl_range_assert(col, context->n - 1); - __nlSparseMatrixAdd(M, row, col, value); + // direct insert into matrix is too slow, so use triplets + EigenTriplet triplet(row, col, value); + context->Mtriplets.push_back(triplet); } } void nlRightHandSideAdd(NLuint rhsindex, NLuint index, NLdouble value) { - __NLContext *context = __nlCurrentContext; - NLdouble* b = context->b; + NLContext *context = __nlCurrentContext; - __nlCheckState(__NL_STATE_MATRIX); + __nlCheckState(NL_STATE_MATRIX); if(context->least_squares) { - __nl_range_assert(index, context->m - 1); - b[rhsindex*context->m + index] += value; + context->b[rhsindex][index] += value; } else { if(!context->variable[index].locked) { index = context->variable[index].index; - __nl_range_assert(index, context->m - 1); - - b[rhsindex*context->m + index] += value; + context->b[rhsindex][index] += value; } } } void nlRightHandSideSet(NLuint rhsindex, NLuint index, NLdouble value) { - __NLContext *context = __nlCurrentContext; - NLdouble* b = context->b; + NLContext *context = __nlCurrentContext; - __nlCheckState(__NL_STATE_MATRIX); + __nlCheckState(NL_STATE_MATRIX); if(context->least_squares) { - __nl_range_assert(index, context->m - 1); - b[rhsindex*context->m + index] = value; + context->b[rhsindex][index] = value; } else { if(!context->variable[index].locked) { index = context->variable[index].index; - __nl_range_assert(index, context->m - 1); - - b[rhsindex*context->m + index] = value; + context->b[rhsindex][index] = value; } } } -void nlBegin(NLenum prim) { +void nlBegin(NLenum prim) +{ switch(prim) { case NL_SYSTEM: { __nlBeginSystem(); @@ -875,12 +412,13 @@ void nlBegin(NLenum prim) { __nlBeginMatrix(); } break; default: { - __nl_assert_not_reached; + assert(0); } } } -void nlEnd(NLenum prim) { +void nlEnd(NLenum prim) +{ switch(prim) { case NL_SYSTEM: { __nlEndSystem(); @@ -889,166 +427,74 @@ void nlEnd(NLenum prim) { __nlEndMatrix(); } break; default: { - __nl_assert_not_reached; + assert(0); } } } -/************************************************************************/ -/* Eigen wrapper */ - -/* Note: Eigen is difficult to call, but it is worth it. */ -/* Here is a driver inspired by A. Sheffer's "cow flattener". */ -static NLboolean __nlFactorize_EIGEN(__NLContext *context, NLint *permutation) { - - /* OpenNL Context */ - __NLSparseMatrix* M = (context->least_squares)? &context->MtM: &context->M; - NLuint n = context->n; - - /* Temporary variables */ - NLuint i, jj, count; - - __nl_assert(!(M->storage & __NL_SYMMETRIC)); - __nl_assert(M->storage & __NL_ROWS); - __nl_assert(M->m == M->n); - - /* Convert M to compressed column format */ - EigenSparseMatrix A(M->m, M->n); - - for(i=0, count=0; i<n; i++) { - __NLRowColumn *Ri = M->row + i; - - for(jj=0; jj<Ri->size; jj++, count++) - A.insert(i, Ri->coeff[jj].index) = Ri->coeff[jj].value; - } - - A.makeCompressed(); - - /* Free M, don't need it anymore at this point */ - __nlSparseMatrixClear(M); +void nlPrintMatrix(void) +{ + NLContext *context = __nlCurrentContext; - /* Performance Sparse LU factorization */ - EigenSparseSolver *sparse_solver = new EigenSparseSolver(); - context->eigen.sparse_solver = sparse_solver; + std::cout << "A:" << context->M << std::endl; - sparse_solver->analyzePattern(A); - sparse_solver->factorize(A); + for(NLuint rhs=0; rhs<context->nb_rhs; rhs++) + std::cout << "b " << rhs << ":" << context->b[rhs] << std::endl; - return (sparse_solver->info() == Eigen::Success); + if (context->MtM.rows() && context->MtM.cols()) + std::cout << "AtA:" << context->MtM << std::endl; } -static NLboolean __nlInvert_EIGEN(__NLContext *context) { +/* Solving */ - /* OpenNL Context */ - NLdouble* b = (context->least_squares)? context->Mtb: context->b; - NLdouble* x = context->x; - NLuint n = context->n, j; +NLboolean nlSolve(NLboolean solveAgain) +{ + NLContext *context = __nlCurrentContext; + NLboolean result = true; - /* Solve each right hand side */ - for(j=0; j<context->nb_rhs; j++, b+=n, x+=n) { - Eigen::Map<Eigen::VectorXd> eigen_b(b, n); + __nlCheckState(NL_STATE_SYSTEM_CONSTRUCTED); - Eigen::VectorXd eigen_x = context->eigen.sparse_solver->solve(eigen_b); - for (NLuint i = 0; i < n; i++) - x[i] = eigen_x[i]; + if (!__nlCurrentContext->solve_again) { + EigenSparseMatrix& M = (context->least_squares)? context->MtM: context->M; - if (context->eigen.sparse_solver->info() != Eigen::Success) - return false; - } + assert(M.rows() == M.cols()); - return true; -} + /* Convert M to compressed column format */ + M.makeCompressed(); -static void __nlFree_EIGEN(__NLContext *context) { - delete context->eigen.sparse_solver; - context->eigen.sparse_solver = NULL; -} - -void nlPrintMatrix(void) { - __NLContext *context = __nlCurrentContext; - __NLSparseMatrix* M = &(context->M); - __NLSparseMatrix* MtM = &(context->MtM); - double *b = context->b; - NLuint i, jj, k; - NLuint m = context->m; - NLuint n = context->n; - __NLRowColumn* Ri = NULL; - double *value = (double*)malloc(sizeof(*value)*(n+m)); - - printf("A:\n"); - for(i=0; i<m; i++) { - Ri = &(M->row[i]); - - memset(value, 0.0, sizeof(*value)*n); - for(jj=0; jj<Ri->size; jj++) - value[Ri->coeff[jj].index] = Ri->coeff[jj].value; - - for (k = 0; k<n; k++) - printf("%.3f ", value[k]); - printf("\n"); - } + /* Perform sparse LU factorization */ + EigenSparseSolver *sparse_solver = new EigenSparseSolver(); + context->sparse_solver = sparse_solver; - for(k=0; k<context->nb_rhs; k++) { - printf("b (%u):\n", k); - for(i=0; i<n; i++) - printf("%f ", b[context->n*k + i]); - printf("\n"); - } - - if(context->alloc_MtM) { - printf("AtA:\n"); - for(i=0; i<n; i++) { - Ri = &(MtM->row[i]); + sparse_solver->analyzePattern(M); + sparse_solver->factorize(M); - memset(value, 0.0, sizeof(*value)*m); - for(jj=0; jj<Ri->size; jj++) - value[Ri->coeff[jj].index] = Ri->coeff[jj].value; + result = (sparse_solver->info() == Eigen::Success); - for (k = 0; k<n; k++) - printf("%.3f ", value[k]); - printf("\n"); - } - - for(k=0; k<context->nb_rhs; k++) { - printf("Mtb (%u):\n", k); - for(i=0; i<n; i++) - printf("%f ", context->Mtb[context->n*k + i]); - printf("\n"); - } - printf("\n"); + /* Free M, don't need it anymore at this point */ + M.resize(0, 0); } - free(value); -} - -/************************************************************************/ -/* nlSolve() driver routine */ - -NLboolean nlSolveAdvanced(NLint *permutation, NLboolean solveAgain) { - NLboolean result = NL_TRUE; - - __nlCheckState(__NL_STATE_SYSTEM_CONSTRUCTED); - - if (!__nlCurrentContext->solve_again) - result = __nlFactorize_EIGEN(__nlCurrentContext, permutation); - if (result) { - result = __nlInvert_EIGEN(__nlCurrentContext); + /* 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(); if (solveAgain) - __nlCurrentContext->solve_again = NL_TRUE; + __nlCurrentContext->solve_again = true; - __nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SYSTEM_SOLVED); + __nlTransition(NL_STATE_SYSTEM_CONSTRUCTED, NL_STATE_SYSTEM_SOLVED); } } return result; } -NLboolean nlSolve() { - return nlSolveAdvanced(NULL, NL_FALSE); -} - diff --git a/source/blender/bmesh/operators/bmo_smooth_laplacian.c b/source/blender/bmesh/operators/bmo_smooth_laplacian.c index a250c716c16..661b1c21892 100644 --- a/source/blender/bmesh/operators/bmo_smooth_laplacian.c +++ b/source/blender/bmesh/operators/bmo_smooth_laplacian.c @@ -551,7 +551,7 @@ void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op) nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (nlSolveAdvanced(NULL, NL_TRUE) ) { + if (nlSolve(NL_TRUE) ) { validate_solution(sys, usex, usey, usez, preserve_volume); } diff --git a/source/blender/editors/armature/meshlaplacian.c b/source/blender/editors/armature/meshlaplacian.c index d2dd0e474c0..956cbe25c3b 100644 --- a/source/blender/editors/armature/meshlaplacian.c +++ b/source/blender/editors/armature/meshlaplacian.c @@ -64,7 +64,7 @@ static void error(const char *str) { printf("error: %s\n", str); } /************************** Laplacian System *****************************/ struct LaplacianSystem { - NLContext context; /* opennl context */ + NLContext *context; /* opennl context */ int totvert, totface; @@ -341,7 +341,7 @@ int laplacian_system_solve(LaplacianSystem *sys) //nlPrintMatrix(); - return nlSolveAdvanced(NULL, NL_TRUE); + return nlSolve(NL_TRUE); } float laplacian_system_get_solution(int v) @@ -1438,7 +1438,7 @@ static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (nlSolveAdvanced(NULL, NL_TRUE)) { + if (nlSolve(NL_TRUE)) { for (z = 0; z < mdb->size; z++) for (y = 0; y < mdb->size; y++) for (x = 0; x < mdb->size; x++) diff --git a/source/blender/editors/armature/reeb.c b/source/blender/editors/armature/reeb.c index f29d15ff416..e10ae4bcf3f 100644 --- a/source/blender/editors/armature/reeb.c +++ b/source/blender/editors/armature/reeb.c @@ -2628,7 +2628,7 @@ int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges) nlEnd(NL_SYSTEM); - success = nlSolveAdvanced(NULL, NL_TRUE); + success = nlSolve(NL_TRUE); if (success) { rval = 1; diff --git a/source/blender/editors/uvedit/uvedit_parametrizer.c b/source/blender/editors/uvedit/uvedit_parametrizer.c index 00615f9bef7..cb387b4af22 100644 --- a/source/blender/editors/uvedit/uvedit_parametrizer.c +++ b/source/blender/editors/uvedit/uvedit_parametrizer.c @@ -193,7 +193,7 @@ typedef struct PChart { union PChartUnion { struct PChartLscm { - NLContext context; + NLContext *context; float *abf_alpha; PVert *pin1, *pin2; } lscm; @@ -2613,7 +2613,7 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart) nlEnd(NL_SYSTEM); - success = nlSolve(); + success = nlSolve(NL_FALSE); if (success) { for (f = chart->faces; f; f = f->nextlink) { @@ -3223,7 +3223,7 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart) nlEnd(NL_SYSTEM); - if (nlSolveAdvanced(NULL, NL_TRUE)) { + if (nlSolve(NL_TRUE)) { p_chart_lscm_load_solution(chart); return P_TRUE; } diff --git a/source/blender/modifiers/intern/MOD_laplaciandeform.c b/source/blender/modifiers/intern/MOD_laplaciandeform.c index 8b0ca7f99fb..916f8b8a36d 100644 --- a/source/blender/modifiers/intern/MOD_laplaciandeform.c +++ b/source/blender/modifiers/intern/MOD_laplaciandeform.c @@ -399,7 +399,6 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3]) sys->context = nlGetCurrent(); nlSolverParameteri(NL_NB_VARIABLES, n); - nlSolverParameteri(NL_SYMMETRIC, NL_FALSE); nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); nlSolverParameteri(NL_NB_ROWS, n + na); nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3); @@ -434,7 +433,7 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3]) } nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (nlSolveAdvanced(NULL, NL_TRUE)) { + if (nlSolve(NL_TRUE)) { sys->has_solution = true; for (j = 1; j <= sys->repeat; j++) { @@ -451,7 +450,7 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3]) nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (!nlSolveAdvanced(NULL, NL_FALSE)) { + if (!nlSolve(NL_FALSE)) { sys->has_solution = false; break; } @@ -495,7 +494,7 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3]) nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (nlSolveAdvanced(NULL, NL_FALSE)) { + if (nlSolve(NL_FALSE)) { sys->has_solution = true; for (j = 1; j <= sys->repeat; j++) { nlBegin(NL_SYSTEM); @@ -510,7 +509,7 @@ static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3]) } nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (!nlSolveAdvanced(NULL, NL_FALSE)) { + if (!nlSolve(NL_FALSE)) { sys->has_solution = false; break; } diff --git a/source/blender/modifiers/intern/MOD_laplaciansmooth.c b/source/blender/modifiers/intern/MOD_laplaciansmooth.c index f04e092b39e..aef28d24a51 100644 --- a/source/blender/modifiers/intern/MOD_laplaciansmooth.c +++ b/source/blender/modifiers/intern/MOD_laplaciansmooth.c @@ -468,7 +468,7 @@ static void laplaciansmoothModifier_do( nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); - if (nlSolveAdvanced(NULL, NL_TRUE)) { + if (nlSolve(NL_TRUE)) { validate_solution(sys, smd->flag, smd->lambda, smd->lambda_border); } } |