griefith/test
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

OpenNL: significantly simplify code using Eigen / STL.

Browse Source
This commit is contained in:
Brecht Van Lommel
2015-11-22 05:25:32 +01:00
parent e6c58df74e
commit 47ce2d7bef
8 changed files with 231 additions and 795 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
intern/opennl/extern/ONL_opennl.h vendored
View File

@@ -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
}

978
intern/opennl/intern/opennl.cpp
View File

@@ -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 */
typedef struct {
NLuint index;
NLdouble value;
} __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);
}
/************************************************************************************/
/* SparseMatrix data structure */
#define __NL_ROWS 1
#define __NL_COLUMNS 2
#define __NL_SYMMETRIC 4
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];
}
}
}
/************************************************************************************/
/* NLContext data structure */
typedef void(*__NLMatrixFunc)(double* x, double* y);
typedef struct {
NLuint index;
NLdouble value;
} NLCoeff;
typedef struct {
NLdouble value[4];
NLdouble value[4];
NLboolean locked;
NLuint index;
__NLRowColumn *a;
} __NLVariable;
NLuint index;
std::vector<NLCoeff> a;
} NLVariable;
#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_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
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 NLContext {
NLenum state;
struct __NLEigenContext {
EigenSparseSolver *sparse_solver;
} eigen;
} __NLContext;
NLuint n;
NLuint m;
static __NLContext* __nlCurrentContext = NULL;
std::vector<EigenTriplet> Mtriplets;
EigenSparseMatrix M;
EigenSparseMatrix MtM;
std::vector<EigenVectorX> b;
std::vector<EigenVectorX> Mtb;
std::vector<EigenVectorX> x;
static void __nlMatrixVectorProd_default(NLdouble* x, NLdouble* y) {
__nlSparseMatrixMult(&(__nlCurrentContext->M), x, y);
}
EigenSparseSolver *sparse_solver;
NLuint nb_variables;
std::vector<NLVariable> variable;
NLContext nlNewContext(void) {
__NLContext* result = __NL_NEW(__NLContext);
result->state = __NL_STATE_INITIAL;
result->matrix_vector_prod = __nlMatrixVectorProd_default;
NLuint nb_rows;
NLuint nb_rhs;
NLboolean least_squares;
NLboolean solve_again;
};
static NLContext* __nlCurrentContext = NULL;
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.resize(context->nb_rhs);
context->x.resize(context->nb_rhs);
context->b = __NL_NEW_ARRAY(NLdouble, m*context->nb_rhs);
context->alloc_b = NL_TRUE;
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();
context->Mtb =
__NL_NEW_ARRAY(NLdouble, context->n*context->nb_rhs);
context->alloc_Mtb = NL_TRUE;
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(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 */
void nlPrintMatrix(void)
{
NLContext *context = __nlCurrentContext;
/* 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) {
std::cout << "A:" << context->M << std::endl;
/* OpenNL Context */
__NLSparseMatrix* M = (context->least_squares)? &context->MtM: &context->M;
NLuint n = context->n;
for(NLuint rhs=0; rhs<context->nb_rhs; rhs++)
std::cout << "b " << rhs << ":" << context->b[rhs] << std::endl;
/* 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);
/* Performance Sparse LU factorization */
EigenSparseSolver *sparse_solver = new EigenSparseSolver();
context->eigen.sparse_solver = sparse_solver;
sparse_solver->analyzePattern(A);
sparse_solver->factorize(A);
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());
/* 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);
}
return true;
}
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");
}
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]);
memset(value, 0.0, sizeof(*value)*m);
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");
}
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(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);
}

2
source/blender/bmesh/operators/bmo_smooth_laplacian.c
View File

@@ -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);
}

6
source/blender/editors/armature/meshlaplacian.c
View File

@@ -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++)

2
source/blender/editors/armature/reeb.c
View File

@@ -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;

6
source/blender/editors/uvedit/uvedit_parametrizer.c
View File

@@ -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;
}

9
source/blender/modifiers/intern/MOD_laplaciandeform.c
View File

@@ -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;
}

2
source/blender/modifiers/intern/MOD_laplaciansmooth.c
View File

@@ -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);
}
}