diff options
author | Brecht Van Lommel <brechtvanlommel@pandora.be> | 2005-10-30 21:38:35 +0300 |
---|---|---|
committer | Brecht Van Lommel <brechtvanlommel@pandora.be> | 2005-10-30 21:38:35 +0300 |
commit | 9cce3710ae84a5fba09333abcf2cb1cdcaaa948b (patch) | |
tree | 56880585f008343cd0bfce5878092f2ca5b011da /intern | |
parent | 1f5193225949b6c42f42feac8ab233e73f24671e (diff) |
Support for adding elements in random positions in an opennl matrix.
Also some code formatting.
Diffstat (limited to 'intern')
-rw-r--r-- | intern/opennl/extern/ONL_opennl.h | 101 | ||||
-rw-r--r-- | intern/opennl/intern/opennl.c | 872 |
2 files changed, 494 insertions, 479 deletions
diff --git a/intern/opennl/extern/ONL_opennl.h b/intern/opennl/extern/ONL_opennl.h index a84a153709d..128f1b137bc 100644 --- a/intern/opennl/extern/ONL_opennl.h +++ b/intern/opennl/extern/ONL_opennl.h @@ -52,33 +52,25 @@ extern "C" { #define NL_VERSION_0_0 1 -/* - * - * Datatypes - * - */ +/* Datatypes */ typedef unsigned int NLenum; typedef unsigned char NLboolean; typedef unsigned int NLbitfield; -typedef void NLvoid; -typedef signed char NLbyte; /* 1-byte signed */ -typedef short NLshort; /* 2-byte signed */ -typedef int NLint; /* 4-byte signed */ +typedef void NLvoid; +typedef signed char NLbyte; /* 1-byte signed */ +typedef short NLshort; /* 2-byte signed */ +typedef int NLint; /* 4-byte signed */ typedef unsigned char NLubyte; /* 1-byte unsigned */ typedef unsigned short NLushort; /* 2-byte unsigned */ typedef unsigned int NLuint; /* 4-byte unsigned */ -typedef int NLsizei; /* 4-byte signed */ -typedef float NLfloat; /* single precision float */ -typedef double NLdouble; /* double precision float */ +typedef int NLsizei; /* 4-byte signed */ +typedef float NLfloat; /* single precision float */ +typedef double NLdouble; /* double precision float */ -typedef void* NLContext ; +typedef void* NLContext; -/* - * - * Constants - * - */ +/* Constants */ #define NL_FALSE 0x0 #define NL_TRUE 0x1 @@ -106,54 +98,51 @@ typedef void* NLContext ; #define NL_RIGHT_HAND_SIDE 0x500 #define NL_ROW_SCALING 0x501 -/* - * Contexts - */ - NLContext nlNewContext(void) ; - void nlDeleteContext(NLContext context) ; - void nlMakeCurrent(NLContext context) ; - NLContext nlGetCurrent(void) ; +/* Contexts */ -/* - * State set/get - */ +NLContext nlNewContext(void); +void nlDeleteContext(NLContext context); +void nlMakeCurrent(NLContext context); +NLContext nlGetCurrent(void); - void nlSolverParameterf(NLenum pname, NLfloat param) ; - void nlSolverParameteri(NLenum pname, NLint param) ; +/* State get/set */ - void nlRowParameterf(NLenum pname, NLfloat param) ; - void nlRowParameteri(NLenum pname, NLint param) ; +void nlSolverParameterf(NLenum pname, NLfloat param); +void nlSolverParameteri(NLenum pname, NLint param); - void nlGetBooleanv(NLenum pname, NLboolean* params) ; - void nlGetFloatv(NLenum pname, NLfloat* params) ; - void nlGetIntergerv(NLenum pname, NLint* params) ; +void nlRowParameterf(NLenum pname, NLfloat param); +void nlRowParameteri(NLenum pname, NLint param); - void nlEnable(NLenum pname) ; - void nlDisable(NLenum pname) ; - NLboolean nlIsEnabled(NLenum pname) ; +void nlGetBooleanv(NLenum pname, NLboolean* params); +void nlGetFloatv(NLenum pname, NLfloat* params); +void nlGetIntergerv(NLenum pname, NLint* params); -/* - * Variables - */ - void nlSetVariable(NLuint index, NLfloat value) ; - NLfloat nlGetVariable(NLuint index) ; - void nlLockVariable(NLuint index) ; - void nlUnlockVariable(NLuint index) ; - NLboolean nlVariableIsLocked(NLuint index) ; +void nlEnable(NLenum pname); +void nlDisable(NLenum pname); +NLboolean nlIsEnabled(NLenum pname); -/* - * Begin/End - */ +/* Variables */ - void nlBegin(NLenum primitive) ; - void nlEnd(NLenum primitive) ; - void nlCoefficient(NLuint index, NLfloat value) ; +void nlSetVariable(NLuint index, NLfloat value); +NLfloat nlGetVariable(NLuint index); +void nlLockVariable(NLuint index); +void nlUnlockVariable(NLuint index); +NLboolean nlVariableIsLocked(NLuint index); -/* - * Solve - */ +/* Begin/End */ + +void nlBegin(NLenum primitive); +void nlEnd(NLenum primitive); +void nlCoefficient(NLuint index, NLfloat value); + +/* Setting random elements matrix/vector - not for least squares! */ + +void nlMatrixAdd(NLuint row, NLuint col, NLfloat value); +void nlRightHandSideAdd(NLuint index, NLfloat value); + +/* Solve */ - NLboolean nlSolve(void) ; +NLboolean nlSolve(void); #ifdef __cplusplus } diff --git a/intern/opennl/intern/opennl.c b/intern/opennl/intern/opennl.c index 3bc87c8dc76..7773582e027 100644 --- a/intern/opennl/intern/opennl.c +++ b/intern/opennl/intern/opennl.c @@ -62,8 +62,8 @@ static void __nl_assertion_failed(char* cond, char* file, int line) { stderr, "OpenNL assertion failed: %s, file:%s, line:%d\n", cond,file,line - ) ; - abort() ; + ); + abort(); } static void __nl_range_assertion_failed( @@ -73,8 +73,8 @@ static void __nl_range_assertion_failed( stderr, "OpenNL range assertion failed: %f in [ %f ... %f ], file:%s, line:%d\n", x, min_val, max_val, file,line - ) ; - abort() ; + ); + abort(); } static void __nl_should_not_have_reached(char* file, int line) { @@ -82,14 +82,14 @@ static void __nl_should_not_have_reached(char* file, int line) { stderr, "OpenNL should not have reached this point: file:%s, line:%d\n", file,line - ) ; - abort() ; + ); + abort(); } #define __nl_assert(x) { \ if(!(x)) { \ - __nl_assertion_failed(#x,__FILE__, __LINE__) ; \ + __nl_assertion_failed(#x,__FILE__, __LINE__); \ } \ } @@ -97,12 +97,12 @@ static void __nl_should_not_have_reached(char* file, int line) { if(((x) < (min_val)) || ((x) > (max_val))) { \ __nl_range_assertion_failed(x, min_val, max_val, \ __FILE__, __LINE__ \ - ) ; \ + ); \ } \ } #define __nl_assert_not_reached { \ - __nl_should_not_have_reached(__FILE__, __LINE__) ; \ + __nl_should_not_have_reached(__FILE__, __LINE__); \ } #ifdef NL_DEBUG @@ -148,73 +148,73 @@ static void __nl_should_not_have_reached(char* file, int line) { /* Dynamic arrays for sparse row/columns */ typedef struct { - NLuint index ; - NLfloat value ; -} __NLCoeff ; + NLuint index; + NLfloat value; +} __NLCoeff; typedef struct { - NLuint size ; - NLuint capacity ; - __NLCoeff* coeff ; -} __NLRowColumn ; + NLuint size; + NLuint capacity; + __NLCoeff* coeff; +} __NLRowColumn; static void __nlRowColumnConstruct(__NLRowColumn* c) { - c->size = 0 ; - c->capacity = 0 ; - c->coeff = NULL ; + c->size = 0; + c->capacity = 0; + c->coeff = NULL; } static void __nlRowColumnDestroy(__NLRowColumn* c) { - __NL_DELETE_ARRAY(c->coeff) ; + __NL_DELETE_ARRAY(c->coeff); #ifdef NL_PARANOID - __NL_CLEAR(__NLRowColumn, c) ; + __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) ; + 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) ; + c->capacity = 4; + c->coeff = __NL_NEW_ARRAY(__NLCoeff, c->capacity); } } static void __nlRowColumnAdd(__NLRowColumn* c, NLint index, NLfloat value) { - NLuint i ; + NLuint i; for(i=0; i<c->size; i++) { if(c->coeff[i].index == (NLuint)index) { - c->coeff[i].value += value ; - return ; + c->coeff[i].value += value; + return; } } if(c->size == c->capacity) { - __nlRowColumnGrow(c) ; + __nlRowColumnGrow(c); } - c->coeff[c->size].index = index ; - c->coeff[c->size].value = value ; - c->size++ ; + 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, NLfloat value) { if(c->size == c->capacity) { - __nlRowColumnGrow(c) ; + __nlRowColumnGrow(c); } - c->coeff[c->size].index = index ; - c->coeff[c->size].value = value ; - c->size++ ; + c->coeff[c->size].index = index; + c->coeff[c->size].value = value; + c->size++; } static void __nlRowColumnZero(__NLRowColumn* c) { - c->size = 0 ; + c->size = 0; } static void __nlRowColumnClear(__NLRowColumn* c) { - c->size = 0 ; - c->capacity = 0 ; - __NL_DELETE_ARRAY(c->coeff) ; + c->size = 0; + c->capacity = 0; + __NL_DELETE_ARRAY(c->coeff); } /************************************************************************************/ @@ -225,115 +225,115 @@ static void __nlRowColumnClear(__NLRowColumn* c) { #define __NL_SYMMETRIC 4 typedef struct { - NLuint m ; - NLuint n ; - NLuint diag_size ; - NLenum storage ; - __NLRowColumn* row ; - __NLRowColumn* column ; - NLfloat* diag ; -} __NLSparseMatrix ; + NLuint m; + NLuint n; + NLuint diag_size; + NLenum storage; + __NLRowColumn* row; + __NLRowColumn* column; + NLfloat* diag; +} __NLSparseMatrix; static void __nlSparseMatrixConstruct( __NLSparseMatrix* M, NLuint m, NLuint n, NLenum storage ) { - NLuint i ; - M->m = m ; - M->n = n ; - M->storage = storage ; + NLuint i; + M->m = m; + M->n = n; + M->storage = storage; if(storage & __NL_ROWS) { - M->row = __NL_NEW_ARRAY(__NLRowColumn, m) ; + M->row = __NL_NEW_ARRAY(__NLRowColumn, m); for(i=0; i<n; i++) { - __nlRowColumnConstruct(&(M->row[i])) ; + __nlRowColumnConstruct(&(M->row[i])); } } else { - M->row = NULL ; + M->row = NULL; } if(storage & __NL_COLUMNS) { - M->column = __NL_NEW_ARRAY(__NLRowColumn, n) ; + M->column = __NL_NEW_ARRAY(__NLRowColumn, n); for(i=0; i<n; i++) { - __nlRowColumnConstruct(&(M->column[i])) ; + __nlRowColumnConstruct(&(M->column[i])); } } else { - M->column = NULL ; + M->column = NULL; } - M->diag_size = MIN(m,n) ; - M->diag = __NL_NEW_ARRAY(NLfloat, M->diag_size) ; + M->diag_size = MIN(m,n); + M->diag = __NL_NEW_ARRAY(NLfloat, M->diag_size); } static void __nlSparseMatrixDestroy(__NLSparseMatrix* M) { - NLuint i ; - __NL_DELETE_ARRAY(M->diag) ; + NLuint i; + __NL_DELETE_ARRAY(M->diag); if(M->storage & __NL_ROWS) { for(i=0; i<M->m; i++) { - __nlRowColumnDestroy(&(M->row[i])) ; + __nlRowColumnDestroy(&(M->row[i])); } - __NL_DELETE_ARRAY(M->row) ; + __NL_DELETE_ARRAY(M->row); } if(M->storage & __NL_COLUMNS) { for(i=0; i<M->n; i++) { - __nlRowColumnDestroy(&(M->column[i])) ; + __nlRowColumnDestroy(&(M->column[i])); } - __NL_DELETE_ARRAY(M->column) ; + __NL_DELETE_ARRAY(M->column); } #ifdef NL_PARANOID - __NL_CLEAR(__NLSparseMatrix,M) ; + __NL_CLEAR(__NLSparseMatrix,M); #endif } static void __nlSparseMatrixAdd( __NLSparseMatrix* M, NLuint i, NLuint j, NLfloat value ) { - __nl_parano_range_assert(i, 0, M->m - 1) ; - __nl_parano_range_assert(j, 0, M->n - 1) ; + __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 ; + return; } if(i == j) { - M->diag[i] += value ; + M->diag[i] += value; } if(M->storage & __NL_ROWS) { - __nlRowColumnAdd(&(M->row[i]), j, value) ; + __nlRowColumnAdd(&(M->row[i]), j, value); } if(M->storage & __NL_COLUMNS) { - __nlRowColumnAdd(&(M->column[j]), i, value) ; + __nlRowColumnAdd(&(M->column[j]), i, value); } } static void __nlSparseMatrixClear( __NLSparseMatrix* M) { - NLuint i ; + NLuint i; if(M->storage & __NL_ROWS) { for(i=0; i<M->m; i++) { - __nlRowColumnClear(&(M->row[i])) ; + __nlRowColumnClear(&(M->row[i])); } } if(M->storage & __NL_COLUMNS) { for(i=0; i<M->n; i++) { - __nlRowColumnClear(&(M->column[i])) ; + __nlRowColumnClear(&(M->column[i])); } } - __NL_CLEAR_ARRAY(NLfloat, M->diag, M->diag_size) ; + __NL_CLEAR_ARRAY(NLfloat, M->diag, M->diag_size); } /* Returns the number of non-zero coefficients */ static NLuint __nlSparseMatrixNNZ( __NLSparseMatrix* M) { - NLuint nnz = 0 ; - NLuint i ; + NLuint nnz = 0; + NLuint i; if(M->storage & __NL_ROWS) { for(i = 0; i<M->m; i++) { - nnz += M->row[i].size ; + nnz += M->row[i].size; } } else if (M->storage & __NL_COLUMNS) { for(i = 0; i<M->n; i++) { - nnz += M->column[i].size ; + nnz += M->column[i].size; } } else { - __nl_assert_not_reached ; + __nl_assert_not_reached; } - return nnz ; + return nnz; } /************************************************************************************/ @@ -342,18 +342,18 @@ static NLuint __nlSparseMatrixNNZ( __NLSparseMatrix* M) { static void __nlSparseMatrix_mult_rows_symmetric( __NLSparseMatrix* A, NLfloat* x, NLfloat* y ) { - NLuint m = A->m ; - NLuint i,ij ; - __NLRowColumn* Ri = NULL ; - __NLCoeff* c = NULL ; + 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]) ; + 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] ; + c = &(Ri->coeff[ij]); + y[i] += c->value * x[c->index]; if(i != c->index) { - y[c->index] += c->value * x[i] ; + y[c->index] += c->value * x[i]; } } } @@ -362,16 +362,16 @@ static void __nlSparseMatrix_mult_rows_symmetric( static void __nlSparseMatrix_mult_rows( __NLSparseMatrix* A, NLfloat* x, NLfloat* y ) { - NLuint m = A->m ; - NLuint i,ij ; - __NLRowColumn* Ri = NULL ; - __NLCoeff* c = NULL ; + 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]) ; + 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] ; + c = &(Ri->coeff[ij]); + y[i] += c->value * x[c->index]; } } } @@ -379,18 +379,18 @@ static void __nlSparseMatrix_mult_rows( static void __nlSparseMatrix_mult_cols_symmetric( __NLSparseMatrix* A, NLfloat* x, NLfloat* y ) { - NLuint n = A->n ; - NLuint j,ii ; - __NLRowColumn* Cj = NULL ; - __NLCoeff* c = NULL ; + 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]) ; + 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] ; + c = &(Cj->coeff[ii]); + y[c->index] += c->value * x[j]; if(j != c->index) { - y[j] += c->value * x[c->index] ; + y[j] += c->value * x[c->index]; } } } @@ -399,16 +399,16 @@ static void __nlSparseMatrix_mult_cols_symmetric( static void __nlSparseMatrix_mult_cols( __NLSparseMatrix* A, NLfloat* x, NLfloat* y ) { - NLuint n = A->n ; - NLuint j,ii ; - __NLRowColumn* Cj = NULL ; - __NLCoeff* c = NULL ; - __NL_CLEAR_ARRAY(NLfloat, y, A->m) ; + NLuint n = A->n; + NLuint j,ii; + __NLRowColumn* Cj = NULL; + __NLCoeff* c = NULL; + __NL_CLEAR_ARRAY(NLfloat, y, A->m); for(j=0; j<n; j++) { - Cj = &(A->column[j]) ; + Cj = &(A->column[j]); for(ii=0; ii<Cj->size; ii++) { - c = &(Cj->coeff[ii]) ; - y[c->index] += c->value * x[j] ; + c = &(Cj->coeff[ii]); + y[c->index] += c->value * x[j]; } } } @@ -419,15 +419,15 @@ static void __nlSparseMatrix_mult_cols( static void __nlSparseMatrixMult(__NLSparseMatrix* A, NLfloat* x, NLfloat* y) { if(A->storage & __NL_ROWS) { if(A->storage & __NL_SYMMETRIC) { - __nlSparseMatrix_mult_rows_symmetric(A, x, y) ; + __nlSparseMatrix_mult_rows_symmetric(A, x, y); } else { - __nlSparseMatrix_mult_rows(A, x, y) ; + __nlSparseMatrix_mult_rows(A, x, y); } } else { if(A->storage & __NL_SYMMETRIC) { - __nlSparseMatrix_mult_cols_symmetric(A, x, y) ; + __nlSparseMatrix_mult_cols_symmetric(A, x, y); } else { - __nlSparseMatrix_mult_cols(A, x, y) ; + __nlSparseMatrix_mult_cols(A, x, y); } } } @@ -435,13 +435,13 @@ static void __nlSparseMatrixMult(__NLSparseMatrix* A, NLfloat* x, NLfloat* y) { /************************************************************************************/ /* NLContext data structure */ -typedef void(*__NLMatrixFunc)(float* x, float* y) ; +typedef void(*__NLMatrixFunc)(float* x, float* y); typedef struct { - NLfloat value ; - NLboolean locked ; - NLuint index ; -} __NLVariable ; + NLfloat value; + NLboolean locked; + NLuint index; +} __NLVariable; #define __NL_STATE_INITIAL 0 #define __NL_STATE_SYSTEM 1 @@ -452,213 +452,213 @@ typedef struct { #define __NL_STATE_SOLVED 6 typedef struct { - NLenum state ; - __NLVariable* variable ; - NLuint n ; - __NLSparseMatrix M ; - __NLRowColumn af ; - __NLRowColumn al ; - __NLRowColumn xl ; - NLfloat* x ; - NLfloat* b ; - NLfloat right_hand_side ; - NLfloat row_scaling ; - NLuint nb_variables ; - NLuint current_row ; - NLboolean least_squares ; - NLboolean symmetric ; - NLboolean normalize_rows ; - NLboolean alloc_M ; - NLboolean alloc_af ; - NLboolean alloc_al ; - NLboolean alloc_xl ; - NLboolean alloc_variable ; - NLboolean alloc_x ; - NLboolean alloc_b ; - NLfloat error ; - __NLMatrixFunc matrix_vector_prod ; -} __NLContext ; - -static __NLContext* __nlCurrentContext = NULL ; + NLenum state; + __NLVariable* variable; + NLuint n; + __NLSparseMatrix M; + __NLRowColumn af; + __NLRowColumn al; + __NLRowColumn xl; + NLfloat* x; + NLfloat* b; + NLfloat right_hand_side; + NLfloat row_scaling; + NLuint nb_variables; + NLuint current_row; + NLboolean least_squares; + NLboolean symmetric; + NLboolean normalize_rows; + NLboolean alloc_M; + NLboolean alloc_af; + NLboolean alloc_al; + NLboolean alloc_xl; + NLboolean alloc_variable; + NLboolean alloc_x; + NLboolean alloc_b; + NLfloat error; + __NLMatrixFunc matrix_vector_prod; +} __NLContext; + +static __NLContext* __nlCurrentContext = NULL; static void __nlMatrixVectorProd_default(NLfloat* x, NLfloat* y) { - __nlSparseMatrixMult(&(__nlCurrentContext->M), x, y) ; + __nlSparseMatrixMult(&(__nlCurrentContext->M), x, y); } NLContext nlNewContext(void) { - __NLContext* result = __NL_NEW(__NLContext) ; - result->state = __NL_STATE_INITIAL ; - result->row_scaling = 1.0 ; - result->right_hand_side = 0.0 ; - result->matrix_vector_prod = __nlMatrixVectorProd_default ; - nlMakeCurrent(result) ; - return result ; + __NLContext* result = __NL_NEW(__NLContext); + result->state = __NL_STATE_INITIAL; + result->row_scaling = 1.0; + result->right_hand_side = 0.0; + result->matrix_vector_prod = __nlMatrixVectorProd_default; + nlMakeCurrent(result); + return result; } void nlDeleteContext(NLContext context_in) { - __NLContext* context = (__NLContext*)(context_in) ; + __NLContext* context = (__NLContext*)(context_in); if(__nlCurrentContext == context) { - __nlCurrentContext = NULL ; + __nlCurrentContext = NULL; } if(context->alloc_M) { - __nlSparseMatrixDestroy(&context->M) ; + __nlSparseMatrixDestroy(&context->M); } if(context->alloc_af) { - __nlRowColumnDestroy(&context->af) ; + __nlRowColumnDestroy(&context->af); } if(context->alloc_al) { - __nlRowColumnDestroy(&context->al) ; + __nlRowColumnDestroy(&context->al); } if(context->alloc_xl) { - __nlRowColumnDestroy(&context->xl) ; + __nlRowColumnDestroy(&context->xl); } if(context->alloc_variable) { - __NL_DELETE_ARRAY(context->variable) ; + __NL_DELETE_ARRAY(context->variable); } if(context->alloc_x) { - __NL_DELETE_ARRAY(context->x) ; + __NL_DELETE_ARRAY(context->x); } if(context->alloc_b) { - __NL_DELETE_ARRAY(context->b) ; + __NL_DELETE_ARRAY(context->b); } #ifdef NL_PARANOID - __NL_CLEAR(__NLContext, context) ; + __NL_CLEAR(__NLContext, context); #endif - __NL_DELETE(context) ; + __NL_DELETE(context); } void nlMakeCurrent(NLContext context) { - __nlCurrentContext = (__NLContext*)(context) ; + __nlCurrentContext = (__NLContext*)(context); } NLContext nlGetCurrent(void) { - return __nlCurrentContext ; + return __nlCurrentContext; } static void __nlCheckState(NLenum state) { - __nl_assert(__nlCurrentContext->state == state) ; + __nl_assert(__nlCurrentContext->state == state); } static void __nlTransition(NLenum from_state, NLenum to_state) { - __nlCheckState(from_state) ; - __nlCurrentContext->state = to_state ; + __nlCheckState(from_state); + __nlCurrentContext->state = to_state; } /************************************************************************************/ /* Get/Set parameters */ void nlSolverParameterf(NLenum pname, NLfloat param) { - __nlCheckState(__NL_STATE_INITIAL) ; + __nlCheckState(__NL_STATE_INITIAL); switch(pname) { case NL_NB_VARIABLES: { - __nl_assert(param > 0) ; - __nlCurrentContext->nb_variables = (NLuint)param ; - } break ; + __nl_assert(param > 0); + __nlCurrentContext->nb_variables = (NLuint)param; + } break; case NL_LEAST_SQUARES: { - __nlCurrentContext->least_squares = (NLboolean)param ; - } break ; + __nlCurrentContext->least_squares = (NLboolean)param; + } break; case NL_SYMMETRIC: { - __nlCurrentContext->symmetric = (NLboolean)param ; + __nlCurrentContext->symmetric = (NLboolean)param; } default: { - __nl_assert_not_reached ; - } break ; + __nl_assert_not_reached; + } break; } } void nlSolverParameteri(NLenum pname, NLint param) { - __nlCheckState(__NL_STATE_INITIAL) ; + __nlCheckState(__NL_STATE_INITIAL); switch(pname) { case NL_NB_VARIABLES: { - __nl_assert(param > 0) ; - __nlCurrentContext->nb_variables = (NLuint)param ; - } break ; + __nl_assert(param > 0); + __nlCurrentContext->nb_variables = (NLuint)param; + } break; case NL_LEAST_SQUARES: { - __nlCurrentContext->least_squares = (NLboolean)param ; - } break ; + __nlCurrentContext->least_squares = (NLboolean)param; + } break; case NL_SYMMETRIC: { - __nlCurrentContext->symmetric = (NLboolean)param ; + __nlCurrentContext->symmetric = (NLboolean)param; } default: { - __nl_assert_not_reached ; - } break ; + __nl_assert_not_reached; + } break; } } void nlRowParameterf(NLenum pname, NLfloat param) { - __nlCheckState(__NL_STATE_MATRIX) ; + __nlCheckState(__NL_STATE_MATRIX); switch(pname) { case NL_RIGHT_HAND_SIDE: { - __nlCurrentContext->right_hand_side = param ; - } break ; + __nlCurrentContext->right_hand_side = param; + } break; case NL_ROW_SCALING: { - __nlCurrentContext->row_scaling = param ; - } break ; + __nlCurrentContext->row_scaling = param; + } break; } } void nlRowParameteri(NLenum pname, NLint param) { - __nlCheckState(__NL_STATE_MATRIX) ; + __nlCheckState(__NL_STATE_MATRIX); switch(pname) { case NL_RIGHT_HAND_SIDE: { - __nlCurrentContext->right_hand_side = (NLfloat)param ; - } break ; + __nlCurrentContext->right_hand_side = (NLfloat)param; + } break; case NL_ROW_SCALING: { - __nlCurrentContext->row_scaling = (NLfloat)param ; - } break ; + __nlCurrentContext->row_scaling = (NLfloat)param; + } break; } } void nlGetBooleanv(NLenum pname, NLboolean* params) { switch(pname) { case NL_LEAST_SQUARES: { - *params = __nlCurrentContext->least_squares ; - } break ; + *params = __nlCurrentContext->least_squares; + } break; case NL_SYMMETRIC: { - *params = __nlCurrentContext->symmetric ; - } break ; + *params = __nlCurrentContext->symmetric; + } break; default: { - __nl_assert_not_reached ; - } break ; + __nl_assert_not_reached; + } break; } } void nlGetFloatv(NLenum pname, NLfloat* params) { switch(pname) { case NL_NB_VARIABLES: { - *params = (NLfloat)(__nlCurrentContext->nb_variables) ; - } break ; + *params = (NLfloat)(__nlCurrentContext->nb_variables); + } break; case NL_LEAST_SQUARES: { - *params = (NLfloat)(__nlCurrentContext->least_squares) ; - } break ; + *params = (NLfloat)(__nlCurrentContext->least_squares); + } break; case NL_SYMMETRIC: { - *params = (NLfloat)(__nlCurrentContext->symmetric) ; - } break ; + *params = (NLfloat)(__nlCurrentContext->symmetric); + } break; case NL_ERROR: { - *params = (NLfloat)(__nlCurrentContext->error) ; - } break ; + *params = (NLfloat)(__nlCurrentContext->error); + } break; default: { - __nl_assert_not_reached ; - } break ; + __nl_assert_not_reached; + } break; } } void nlGetIntergerv(NLenum pname, NLint* params) { switch(pname) { case NL_NB_VARIABLES: { - *params = (NLint)(__nlCurrentContext->nb_variables) ; - } break ; + *params = (NLint)(__nlCurrentContext->nb_variables); + } break; case NL_LEAST_SQUARES: { - *params = (NLint)(__nlCurrentContext->least_squares) ; - } break ; + *params = (NLint)(__nlCurrentContext->least_squares); + } break; case NL_SYMMETRIC: { - *params = (NLint)(__nlCurrentContext->symmetric) ; - } break ; + *params = (NLint)(__nlCurrentContext->symmetric); + } break; default: { - __nl_assert_not_reached ; - } break ; + __nl_assert_not_reached; + } break; } } @@ -668,11 +668,11 @@ void nlGetIntergerv(NLenum pname, NLint* params) { void nlEnable(NLenum pname) { switch(pname) { case NL_NORMALIZE_ROWS: { - __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW) ; - __nlCurrentContext->normalize_rows = NL_TRUE ; - } break ; + __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW); + __nlCurrentContext->normalize_rows = NL_TRUE; + } break; default: { - __nl_assert_not_reached ; + __nl_assert_not_reached; } } } @@ -680,11 +680,11 @@ void nlEnable(NLenum pname) { void nlDisable(NLenum pname) { switch(pname) { case NL_NORMALIZE_ROWS: { - __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW) ; - __nlCurrentContext->normalize_rows = NL_FALSE ; - } break ; + __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW); + __nlCurrentContext->normalize_rows = NL_FALSE; + } break; default: { - __nl_assert_not_reached ; + __nl_assert_not_reached; } } } @@ -692,217 +692,219 @@ void nlDisable(NLenum pname) { NLboolean nlIsEnabled(NLenum pname) { switch(pname) { case NL_NORMALIZE_ROWS: { - return __nlCurrentContext->normalize_rows ; - } break ; + return __nlCurrentContext->normalize_rows; + } break; default: { - __nl_assert_not_reached ; + __nl_assert_not_reached; } } - return NL_FALSE ; + return NL_FALSE; } /************************************************************************************/ /* Get/Set Lock/Unlock variables */ void nlSetVariable(NLuint index, NLfloat value) { - __nlCheckState(__NL_STATE_SYSTEM) ; - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1) ; - __nlCurrentContext->variable[index].value = value ; + __nlCheckState(__NL_STATE_SYSTEM); + __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); + __nlCurrentContext->variable[index].value = value; } NLfloat nlGetVariable(NLuint index) { - __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL) ; - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1) ; - return __nlCurrentContext->variable[index].value ; + __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL); + __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); + return __nlCurrentContext->variable[index].value; } void nlLockVariable(NLuint index) { - __nlCheckState(__NL_STATE_SYSTEM) ; - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1) ; - __nlCurrentContext->variable[index].locked = NL_TRUE ; + __nlCheckState(__NL_STATE_SYSTEM); + __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); + __nlCurrentContext->variable[index].locked = NL_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 ; + __nlCheckState(__NL_STATE_SYSTEM); + __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); + __nlCurrentContext->variable[index].locked = NL_FALSE; } NLboolean nlVariableIsLocked(NLuint index) { - __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL) ; - __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1) ; - return __nlCurrentContext->variable[index].locked ; + __nl_assert(__nlCurrentContext->state != __NL_STATE_INITIAL); + __nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables - 1); + return __nlCurrentContext->variable[index].locked; } /************************************************************************************/ /* System construction */ static void __nlVariablesToVector() { - NLuint i ; - __nl_assert(__nlCurrentContext->alloc_x) ; - __nl_assert(__nlCurrentContext->alloc_variable) ; + NLuint i; + __nl_assert(__nlCurrentContext->alloc_x); + __nl_assert(__nlCurrentContext->alloc_variable); for(i=0; i<__nlCurrentContext->nb_variables; i++) { - __NLVariable* v = &(__nlCurrentContext->variable[i]) ; + __NLVariable* v = &(__nlCurrentContext->variable[i]); if(!v->locked) { - __nl_assert(v->index < __nlCurrentContext->n) ; - __nlCurrentContext->x[v->index] = v->value ; + __nl_assert(v->index < __nlCurrentContext->n); + __nlCurrentContext->x[v->index] = v->value; } } } static void __nlVectorToVariables() { - NLuint i ; - __nl_assert(__nlCurrentContext->alloc_x) ; - __nl_assert(__nlCurrentContext->alloc_variable) ; + NLuint i; + __nl_assert(__nlCurrentContext->alloc_x); + __nl_assert(__nlCurrentContext->alloc_variable); for(i=0; i<__nlCurrentContext->nb_variables; i++) { - __NLVariable* v = &(__nlCurrentContext->variable[i]) ; + __NLVariable* v = &(__nlCurrentContext->variable[i]); if(!v->locked) { - __nl_assert(v->index < __nlCurrentContext->n) ; - v->value = __nlCurrentContext->x[v->index] ; + __nl_assert(v->index < __nlCurrentContext->n); + v->value = __nlCurrentContext->x[v->index]; } } } static void __nlBeginSystem() { - __nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM) ; - __nl_assert(__nlCurrentContext->nb_variables > 0) ; + __nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM); + __nl_assert(__nlCurrentContext->nb_variables > 0); __nlCurrentContext->variable = __NL_NEW_ARRAY( __NLVariable, __nlCurrentContext->nb_variables - ) ; - __nlCurrentContext->alloc_variable = NL_TRUE ; + ); + __nlCurrentContext->alloc_variable = NL_TRUE; } static void __nlEndSystem() { - __nlTransition(__NL_STATE_MATRIX_CONSTRUCTED, __NL_STATE_SYSTEM_CONSTRUCTED) ; + __nlTransition(__NL_STATE_MATRIX_CONSTRUCTED, __NL_STATE_SYSTEM_CONSTRUCTED); } static void __nlBeginMatrix() { - NLuint i ; - NLuint n = 0 ; - NLenum storage = __NL_ROWS ; + NLuint i; + NLuint n = 0; + NLenum storage = __NL_ROWS; - __nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX) ; + __nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX); for(i=0; i<__nlCurrentContext->nb_variables; i++) { if(!__nlCurrentContext->variable[i].locked) { - __nlCurrentContext->variable[i].index = n ; - n++ ; + __nlCurrentContext->variable[i].index = n; + n++; } else { - __nlCurrentContext->variable[i].index = ~0 ; + __nlCurrentContext->variable[i].index = ~0; } } - __nlCurrentContext->n = n ; + __nlCurrentContext->n = n; /* a least squares problem results in a symmetric matrix */ if(__nlCurrentContext->least_squares) { - __nlCurrentContext->symmetric = NL_TRUE ; + __nlCurrentContext->symmetric = NL_TRUE; } if(__nlCurrentContext->symmetric) { - storage = (storage | __NL_SYMMETRIC) ; + storage = (storage | __NL_SYMMETRIC); } /* SuperLU storage does not support symmetric storage */ - storage = (storage & ~__NL_SYMMETRIC) ; + storage = (storage & ~__NL_SYMMETRIC); - __nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage) ; - __nlCurrentContext->alloc_M = NL_TRUE ; + __nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage); + __nlCurrentContext->alloc_M = NL_TRUE; - __nlCurrentContext->x = __NL_NEW_ARRAY(NLfloat, n) ; - __nlCurrentContext->alloc_x = NL_TRUE ; + __nlCurrentContext->x = __NL_NEW_ARRAY(NLfloat, n); + __nlCurrentContext->alloc_x = NL_TRUE; - __nlCurrentContext->b = __NL_NEW_ARRAY(NLfloat, n) ; - __nlCurrentContext->alloc_b = NL_TRUE ; + __nlCurrentContext->b = __NL_NEW_ARRAY(NLfloat, n); + __nlCurrentContext->alloc_b = NL_TRUE; - __nlVariablesToVector() ; + __nlVariablesToVector(); - __nlRowColumnConstruct(&__nlCurrentContext->af) ; - __nlCurrentContext->alloc_af = NL_TRUE ; - __nlRowColumnConstruct(&__nlCurrentContext->al) ; - __nlCurrentContext->alloc_al = NL_TRUE ; - __nlRowColumnConstruct(&__nlCurrentContext->xl) ; - __nlCurrentContext->alloc_xl = NL_TRUE ; + __nlRowColumnConstruct(&__nlCurrentContext->af); + __nlCurrentContext->alloc_af = NL_TRUE; + __nlRowColumnConstruct(&__nlCurrentContext->al); + __nlCurrentContext->alloc_al = NL_TRUE; + __nlRowColumnConstruct(&__nlCurrentContext->xl); + __nlCurrentContext->alloc_xl = NL_TRUE; - __nlCurrentContext->current_row = 0 ; + __nlCurrentContext->current_row = 0; } static void __nlEndMatrix() { - __nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED) ; + __nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED); - __nlRowColumnDestroy(&__nlCurrentContext->af) ; - __nlCurrentContext->alloc_af = NL_FALSE ; - __nlRowColumnDestroy(&__nlCurrentContext->al) ; - __nlCurrentContext->alloc_al = NL_FALSE ; - __nlRowColumnDestroy(&__nlCurrentContext->xl) ; - __nlCurrentContext->alloc_al = NL_FALSE ; + __nlRowColumnDestroy(&__nlCurrentContext->af); + __nlCurrentContext->alloc_af = NL_FALSE; + __nlRowColumnDestroy(&__nlCurrentContext->al); + __nlCurrentContext->alloc_al = NL_FALSE; + __nlRowColumnDestroy(&__nlCurrentContext->xl); + __nlCurrentContext->alloc_al = NL_FALSE; +#if 0 if(!__nlCurrentContext->least_squares) { __nl_assert( __nlCurrentContext->current_row == __nlCurrentContext->n - ) ; + ); } +#endif } static void __nlBeginRow() { - __nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW) ; - __nlRowColumnZero(&__nlCurrentContext->af) ; - __nlRowColumnZero(&__nlCurrentContext->al) ; - __nlRowColumnZero(&__nlCurrentContext->xl) ; + __nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW); + __nlRowColumnZero(&__nlCurrentContext->af); + __nlRowColumnZero(&__nlCurrentContext->al); + __nlRowColumnZero(&__nlCurrentContext->xl); } static void __nlScaleRow(NLfloat s) { - __NLRowColumn* af = &__nlCurrentContext->af ; - __NLRowColumn* al = &__nlCurrentContext->al ; - NLuint nf = af->size ; - NLuint nl = al->size ; - NLuint i ; + __NLRowColumn* af = &__nlCurrentContext->af; + __NLRowColumn* al = &__nlCurrentContext->al; + NLuint nf = af->size; + NLuint nl = al->size; + NLuint i; for(i=0; i<nf; i++) { - af->coeff[i].value *= s ; + af->coeff[i].value *= s; } for(i=0; i<nl; i++) { - al->coeff[i].value *= s ; + al->coeff[i].value *= s; } - __nlCurrentContext->right_hand_side *= s ; + __nlCurrentContext->right_hand_side *= s; } static void __nlNormalizeRow(NLfloat weight) { - __NLRowColumn* af = &__nlCurrentContext->af ; - __NLRowColumn* al = &__nlCurrentContext->al ; - NLuint nf = af->size ; - NLuint nl = al->size ; - NLuint i ; - NLfloat norm = 0.0 ; + __NLRowColumn* af = &__nlCurrentContext->af; + __NLRowColumn* al = &__nlCurrentContext->al; + NLuint nf = af->size; + NLuint nl = al->size; + NLuint i; + NLfloat norm = 0.0; for(i=0; i<nf; i++) { - norm += af->coeff[i].value * af->coeff[i].value ; + norm += af->coeff[i].value * af->coeff[i].value; } for(i=0; i<nl; i++) { - norm += al->coeff[i].value * al->coeff[i].value ; + norm += al->coeff[i].value * al->coeff[i].value; } - norm = sqrt(norm) ; - __nlScaleRow(weight / norm) ; + norm = sqrt(norm); + __nlScaleRow(weight / norm); } static void __nlEndRow() { - __NLRowColumn* af = &__nlCurrentContext->af ; - __NLRowColumn* al = &__nlCurrentContext->al ; - __NLRowColumn* xl = &__nlCurrentContext->xl ; - __NLSparseMatrix* M = &__nlCurrentContext->M ; - NLfloat* b = __nlCurrentContext->b ; - NLuint nf = af->size ; - NLuint nl = al->size ; - NLuint current_row = __nlCurrentContext->current_row ; - NLuint i ; - NLuint j ; - NLfloat S ; - __nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX) ; + __NLRowColumn* af = &__nlCurrentContext->af; + __NLRowColumn* al = &__nlCurrentContext->al; + __NLRowColumn* xl = &__nlCurrentContext->xl; + __NLSparseMatrix* M = &__nlCurrentContext->M; + NLfloat* b = __nlCurrentContext->b; + NLuint nf = af->size; + NLuint nl = al->size; + NLuint current_row = __nlCurrentContext->current_row; + NLuint i; + NLuint j; + NLfloat S; + __nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX); if(__nlCurrentContext->normalize_rows) { - __nlNormalizeRow(__nlCurrentContext->row_scaling) ; + __nlNormalizeRow(__nlCurrentContext->row_scaling); } else { - __nlScaleRow(__nlCurrentContext->row_scaling) ; + __nlScaleRow(__nlCurrentContext->row_scaling); } if(__nlCurrentContext->least_squares) { @@ -911,59 +913,81 @@ static void __nlEndRow() { __nlSparseMatrixAdd( M, af->coeff[i].index, af->coeff[j].index, af->coeff[i].value * af->coeff[j].value - ) ; + ); } } - S = -__nlCurrentContext->right_hand_side ; + S = -__nlCurrentContext->right_hand_side; for(j=0; j<nl; j++) { - S += al->coeff[j].value * xl->coeff[j].value ; + S += al->coeff[j].value * xl->coeff[j].value; } for(i=0; i<nf; i++) { - b[ af->coeff[i].index ] -= af->coeff[i].value * S ; + b[ af->coeff[i].index ] -= af->coeff[i].value * S; } } else { for(i=0; i<nf; i++) { __nlSparseMatrixAdd( M, current_row, af->coeff[i].index, af->coeff[i].value - ) ; + ); } - b[current_row] = -__nlCurrentContext->right_hand_side ; + b[current_row] = -__nlCurrentContext->right_hand_side; for(i=0; i<nl; i++) { - b[current_row] -= al->coeff[i].value * xl->coeff[i].value ; + b[current_row] -= al->coeff[i].value * xl->coeff[i].value; } } - __nlCurrentContext->current_row++ ; - __nlCurrentContext->right_hand_side = 0.0 ; - __nlCurrentContext->row_scaling = 1.0 ; + __nlCurrentContext->current_row++; + __nlCurrentContext->right_hand_side = 0.0; + __nlCurrentContext->row_scaling = 1.0; +} + +void nlMatrixAdd(NLuint row, NLuint col, NLfloat value) +{ + __NLSparseMatrix* M = &__nlCurrentContext->M; + __nlCheckState(__NL_STATE_MATRIX); + __nl_range_assert(row, 0, __nlCurrentContext->n - 1); + __nl_range_assert(col, 0, __nlCurrentContext->nb_variables - 1); + __nl_assert(!__nlCurrentContext->least_squares); + + __nlSparseMatrixAdd(M, row, col, value); +} + +void nlRightHandSideAdd(NLuint index, NLfloat value) +{ + NLfloat* b = __nlCurrentContext->b; + + __nlCheckState(__NL_STATE_MATRIX); + __nl_range_assert(index, 0, __nlCurrentContext->n - 1); + __nl_assert(!__nlCurrentContext->least_squares); + + b[index] += value; } void nlCoefficient(NLuint index, NLfloat value) { __NLVariable* v; unsigned int zero= 0; - __nlCheckState(__NL_STATE_ROW) ; - __nl_range_assert(index, zero, __nlCurrentContext->nb_variables - 1) ; - v = &(__nlCurrentContext->variable[index]) ; + __nlCheckState(__NL_STATE_ROW); + __nl_range_assert(index, zero, __nlCurrentContext->nb_variables - 1); + v = &(__nlCurrentContext->variable[index]); if(v->locked) { - __nlRowColumnAppend(&(__nlCurrentContext->al), 0, value) ; - __nlRowColumnAppend(&(__nlCurrentContext->xl), 0, v->value) ; + __nlRowColumnAppend(&(__nlCurrentContext->al), 0, value); + __nlRowColumnAppend(&(__nlCurrentContext->xl), 0, v->value); } else { - __nlRowColumnAppend(&(__nlCurrentContext->af), v->index, value) ; + __nlRowColumnAppend(&(__nlCurrentContext->af), v->index, value); } } void nlBegin(NLenum prim) { switch(prim) { case NL_SYSTEM: { - __nlBeginSystem() ; - } break ; + __nlBeginSystem(); + } break; case NL_MATRIX: { - __nlBeginMatrix() ; - } break ; + __nlBeginMatrix(); + } break; case NL_ROW: { - __nlBeginRow() ; - } break ; + __nlBeginRow(); + } break; default: { - __nl_assert_not_reached ; + __nl_assert_not_reached; } } } @@ -971,16 +995,16 @@ void nlBegin(NLenum prim) { void nlEnd(NLenum prim) { switch(prim) { case NL_SYSTEM: { - __nlEndSystem() ; - } break ; + __nlEndSystem(); + } break; case NL_MATRIX: { - __nlEndMatrix() ; - } break ; + __nlEndMatrix(); + } break; case NL_ROW: { - __nlEndRow() ; - } break ; + __nlEndRow(); + } break; default: { - __nl_assert_not_reached ; + __nl_assert_not_reached; } } } @@ -993,41 +1017,41 @@ void nlEnd(NLenum prim) { static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) { /* OpenNL Context */ - __NLSparseMatrix* M = &(__nlCurrentContext->M) ; - NLfloat* b = __nlCurrentContext->b ; - NLfloat* x = __nlCurrentContext->x ; + __NLSparseMatrix* M = &(__nlCurrentContext->M); + NLfloat* b = __nlCurrentContext->b; + NLfloat* x = __nlCurrentContext->x; /* Compressed Row Storage matrix representation */ - NLuint n = __nlCurrentContext->n ; - NLuint nnz = __nlSparseMatrixNNZ(M) ; /* Number of Non-Zero coeffs */ - NLint* xa = __NL_NEW_ARRAY(NLint, n+1) ; - NLfloat* rhs = __NL_NEW_ARRAY(NLfloat, n) ; - NLfloat* a = __NL_NEW_ARRAY(NLfloat, nnz) ; - NLint* asub = __NL_NEW_ARRAY(NLint, nnz) ; + NLuint n = __nlCurrentContext->n; + NLuint nnz = __nlSparseMatrixNNZ(M); /* Number of Non-Zero coeffs */ + NLint* xa = __NL_NEW_ARRAY(NLint, n+1); + NLfloat* rhs = __NL_NEW_ARRAY(NLfloat, n); + NLfloat* a = __NL_NEW_ARRAY(NLfloat, nnz); + NLint* asub = __NL_NEW_ARRAY(NLint, nnz); /* Permutation vector */ - NLint* perm_r = __NL_NEW_ARRAY(NLint, n) ; - NLint* perm = __NL_NEW_ARRAY(NLint, n) ; + NLint* perm_r = __NL_NEW_ARRAY(NLint, n); + NLint* perm = __NL_NEW_ARRAY(NLint, n); /* SuperLU variables */ - SuperMatrix A, B ; /* System */ - SuperMatrix L, U ; /* Inverse of A */ - NLint info ; /* status code */ - DNformat *vals = NULL ; /* access to result */ - float *rvals = NULL ; /* access to result */ + SuperMatrix A, B; /* System */ + SuperMatrix L, U; /* Inverse of A */ + NLint info; /* status code */ + DNformat *vals = NULL; /* access to result */ + float *rvals = NULL; /* access to result */ /* SuperLU options and stats */ - superlu_options_t options ; - SuperLUStat_t stat ; + superlu_options_t options; + SuperLUStat_t stat; /* Temporary variables */ - __NLRowColumn* Ri = NULL ; - NLuint i,jj,count ; + __NLRowColumn* Ri = NULL; + NLuint i,jj,count; - __nl_assert(!(M->storage & __NL_SYMMETRIC)) ; - __nl_assert(M->storage & __NL_ROWS) ; - __nl_assert(M->m == M->n) ; + __nl_assert(!(M->storage & __NL_SYMMETRIC)); + __nl_assert(M->storage & __NL_ROWS); + __nl_assert(M->m == M->n); /* @@ -1036,20 +1060,20 @@ static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) { * ------------------------------------------------------- */ - count = 0 ; + count = 0; for(i=0; i<n; i++) { - Ri = &(M->row[i]) ; - xa[i] = count ; + Ri = &(M->row[i]); + xa[i] = count; for(jj=0; jj<Ri->size; jj++) { - a[count] = Ri->coeff[jj].value ; - asub[count] = Ri->coeff[jj].index ; - count++ ; + a[count] = Ri->coeff[jj].value; + asub[count] = Ri->coeff[jj].index; + count++; } } - xa[n] = nnz ; + xa[n] = nnz; /* Save memory for SuperLU */ - __nlSparseMatrixClear(M) ; + __nlSparseMatrixClear(M); /* @@ -1079,15 +1103,15 @@ static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) { * 2 -> re-ordering for A^t+A * 3 -> approximate minimum degree ordering */ - get_perm_c(do_perm ? 3 : 0, &A, perm) ; + get_perm_c(do_perm ? 3 : 0, &A, perm); /* Step 4: call SuperLU main routine * --------------------------------- */ - set_default_options(&options) ; - options.ColPerm = MY_PERMC ; - StatInit(&stat) ; + set_default_options(&options); + options.ColPerm = MY_PERMC; + StatInit(&stat); sgssv(&options, &A, perm, perm_r, &L, &U, &B, &stat, &info); @@ -1112,8 +1136,8 @@ static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) { * needs to be deallocated (the arrays have been allocated * by us). */ - Destroy_SuperMatrix_Store(&A) ; - Destroy_SuperMatrix_Store(&B) ; + Destroy_SuperMatrix_Store(&A); + Destroy_SuperMatrix_Store(&B); /* @@ -1123,14 +1147,16 @@ static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) { Destroy_SuperNode_Matrix(&L); Destroy_CompCol_Matrix(&U); - __NL_DELETE_ARRAY(xa) ; - __NL_DELETE_ARRAY(rhs) ; - __NL_DELETE_ARRAY(a) ; - __NL_DELETE_ARRAY(asub) ; - __NL_DELETE_ARRAY(perm_r) ; - __NL_DELETE_ARRAY(perm) ; + StatFree(&stat); + + __NL_DELETE_ARRAY(xa); + __NL_DELETE_ARRAY(rhs); + __NL_DELETE_ARRAY(a); + __NL_DELETE_ARRAY(asub); + __NL_DELETE_ARRAY(perm_r); + __NL_DELETE_ARRAY(perm); - return (info == 0) ; + return (info == 0); } @@ -1138,14 +1164,14 @@ static NLboolean __nlSolve_SUPERLU( NLboolean do_perm) { /* nlSolve() driver routine */ NLboolean nlSolve(void) { - NLboolean result = NL_TRUE ; + NLboolean result = NL_TRUE; - __nlCheckState(__NL_STATE_SYSTEM_CONSTRUCTED) ; - result = __nlSolve_SUPERLU(NL_TRUE) ; + __nlCheckState(__NL_STATE_SYSTEM_CONSTRUCTED); + result = __nlSolve_SUPERLU(NL_TRUE); - __nlVectorToVariables() ; - __nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SOLVED) ; + __nlVectorToVariables(); + __nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SOLVED); - return result ; + return result; } |