Updates to opennl for mesh laplacian matrix building, to make matrix
authorBrecht Van Lommel <brechtvanlommel@pandora.be>
Sat, 28 Jul 2007 14:39:30 +0000 (14:39 +0000)
committerBrecht Van Lommel <brechtvanlommel@pandora.be>
Sat, 28 Jul 2007 14:39:30 +0000 (14:39 +0000)
building with random access work together with variable locking.

intern/opennl/extern/ONL_opennl.h
intern/opennl/intern/opennl.c

index 345cf0dc7170f8f77b9a98b57a6a375277f24c0b..7f501629290215e3520835400987dadf5f8e8b9c 100644 (file)
@@ -131,10 +131,12 @@ void nlBegin(NLenum primitive);
 void nlEnd(NLenum primitive);
 void nlCoefficient(NLuint index, NLfloat value);
 
-/* Setting random elements matrix/vector - not for least squares! */
+/* Setting random elements matrix/vector - not supported for
+   least squares! */
 
 void nlMatrixAdd(NLuint row, NLuint col, NLfloat value);
 void nlRightHandSideAdd(NLuint index, NLfloat value);
+void nlRightHandSideSet(NLuint index, NLfloat value);
 
 /* Solve */
 
index c5518731c6bccf3e106cebbde1ad36d9715ad2f8..836c9d01e60e8815bba76ddd7e095cca91ed14d1 100644 (file)
@@ -441,6 +441,7 @@ typedef struct {
        NLfloat  value;
        NLboolean locked;
        NLuint  index;
+       __NLRowColumn *a;
 } __NLVariable;
 
 #define __NL_STATE_INITIAL                             0
@@ -453,13 +454,13 @@ typedef struct {
 
 typedef struct {
        NLenum             state;
-       __NLVariable*   variable;
        NLuint             n;
+       __NLVariable*   variable;
+       NLfloat*                b;
        __NLSparseMatrix M;
        __NLRowColumn   af;
        __NLRowColumn   al;
        NLfloat*                x;
-       NLfloat*                b;
        NLfloat          right_hand_side;
        NLuint             nb_variables;
        NLuint             current_row;
@@ -472,8 +473,8 @@ typedef struct {
        NLboolean               alloc_variable;
        NLboolean               alloc_x;
        NLboolean               alloc_b;
-       NLfloat          error;
-       __NLMatrixFunc   matrix_vector_prod;
+       NLfloat                 error;
+       __NLMatrixFunc  matrix_vector_prod;
 
        struct __NLSuperLUContext {
                NLboolean alloc_slu;
@@ -503,6 +504,8 @@ static void __nlFree_SUPERLU(__NLContext *context);
 
 void nlDeleteContext(NLContext context_in) {
        __NLContext* context = (__NLContext*)(context_in);
+       int i;
+
        if(__nlCurrentContext == context) {
                __nlCurrentContext = NULL;
        }
@@ -516,14 +519,19 @@ void nlDeleteContext(NLContext context_in) {
                __nlRowColumnDestroy(&context->al);
        }
        if(context->alloc_variable) {
-               __NL_DELETE_ARRAY(context->variable);
-       }
-       if(context->alloc_x) {
-               __NL_DELETE_ARRAY(context->x);
+               for(i=0; i<context->nb_variables; i++) {
+                       if(context->variable[i].a) {
+                               __nlRowColumnDestroy(context->variable[i].a);
+                               __NL_DELETE(context->variable[i].a);
+                       }
+               }
        }
        if(context->alloc_b) {
                __NL_DELETE_ARRAY(context->b);
        }
+       if(context->alloc_x) {
+               __NL_DELETE_ARRAY(context->x);
+       }
        if (context->slu.alloc_slu) {
                __nlFree_SUPERLU(context);
        }
@@ -727,8 +735,10 @@ NLboolean nlVariableIsLocked(NLuint index) {
 
 static void __nlVariablesToVector() {
        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]);
                if(!v->locked) {
@@ -740,8 +750,10 @@ static void __nlVariablesToVector() {
 
 static void __nlVectorToVariables() {
        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]);
                if(!v->locked) {
@@ -760,8 +772,8 @@ static void __nlBeginSystem() {
                __nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM);
 
                __nlCurrentContext->variable = __NL_NEW_ARRAY(
-                       __NLVariable, __nlCurrentContext->nb_variables
-               );
+                       __NLVariable, __nlCurrentContext->nb_variables);
+               
                __nlCurrentContext->alloc_variable = NL_TRUE;
        }
 }
@@ -771,69 +783,93 @@ static void __nlEndSystem() {
 }
 
 static void __nlBeginMatrix() {
-       NLuint i;
+       NLuint i, j;
        NLuint n = 0;
        NLenum storage = __NL_ROWS;
+       __NLContext *context = __nlCurrentContext;
 
        __nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX);
 
-       if (!__nlCurrentContext->solve_again) {
-               for(i=0; i<__nlCurrentContext->nb_variables; i++) {
-                       if(!__nlCurrentContext->variable[i].locked)
-                               __nlCurrentContext->variable[i].index = n++;
+       if (!context->solve_again) {
+               for(i=0; i<context->nb_variables; i++) {
+                       if(context->variable[i].locked) {
+                               context->variable[i].index = ~0;
+                               context->variable[i].a = __NL_NEW(__NLRowColumn);
+                               __nlRowColumnConstruct(context->variable[i].a);
+                       }
                        else
-                               __nlCurrentContext->variable[i].index = ~0;
+                               context->variable[i].index = n++;
                }
 
-               __nlCurrentContext->n = n;
+               context->n = n;
 
                /* a least squares problem results in a symmetric matrix */
-               if(__nlCurrentContext->least_squares)
-                       __nlCurrentContext->symmetric = NL_TRUE;
+               if(context->least_squares)
+                       context->symmetric = NL_TRUE;
 
-               if(__nlCurrentContext->symmetric)
+               if(context->symmetric)
                        storage = (storage | __NL_SYMMETRIC);
 
                /* SuperLU storage does not support symmetric storage */
                storage = (storage & ~__NL_SYMMETRIC);
 
-               __nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage);
-               __nlCurrentContext->alloc_M = NL_TRUE;
+               __nlSparseMatrixConstruct(&context->M, n, n, storage);
+               context->alloc_M = 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;
+               context->b = __NL_NEW_ARRAY(NLfloat, n);
+               context->alloc_b = NL_TRUE;
+
+               context->x = __NL_NEW_ARRAY(NLfloat, n);
+               context->alloc_x = NL_TRUE;
        }
        else {
                /* need to recompute b only, A is not constructed anymore */
-               __NL_CLEAR_ARRAY(NLfloat, __nlCurrentContext->b, __nlCurrentContext->n);
+               __NL_CLEAR_ARRAY(NLfloat, context->b, context->n);
        }
 
        __nlVariablesToVector();
 
-       __nlRowColumnConstruct(&__nlCurrentContext->af);
-       __nlCurrentContext->alloc_af = NL_TRUE;
-       __nlRowColumnConstruct(&__nlCurrentContext->al);
-       __nlCurrentContext->alloc_al = NL_TRUE;
+       __nlRowColumnConstruct(&context->af);
+       context->alloc_af = NL_TRUE;
+       __nlRowColumnConstruct(&context->al);
+       context->alloc_al = NL_TRUE;
 
-       __nlCurrentContext->current_row = 0;
+       context->current_row = 0;
 }
 
 static void __nlEndMatrix() {
+       __NLContext *context = __nlCurrentContext;
+       __NLVariable *variable;
+       __NLRowColumn *a;
+       NLfloat *b;
+       NLuint i, j;
+
        __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(&context->af);
+       context->alloc_af = NL_FALSE;
+       __nlRowColumnDestroy(&context->al);
+       context->alloc_al = NL_FALSE;
        
+       b = context->b;
+
+       for(i=0; i<__nlCurrentContext->nb_variables; i++) {
+               variable = &(context->variable[i]);
+
+               if(variable->locked) {
+                       a = variable->a;
+
+                       for(j=0; j<a->size; j++) {
+                               b[a->coeff[j].index] -= a->coeff[j].value*variable->value;
+                       }
+               }
+       }
+
 #if 0
-       if(!__nlCurrentContext->least_squares) {
+       if(!context->least_squares) {
                __nl_assert(
-                       __nlCurrentContext->current_row == 
-                       __nlCurrentContext->n
+                       context->current_row == 
+                       context->n
                );
        }
 #endif
@@ -895,24 +931,59 @@ static void __nlEndRow() {
 
 void nlMatrixAdd(NLuint row, NLuint col, NLfloat value)
 {
-       __NLSparseMatrix* M  = &__nlCurrentContext->M;
+       __NLContext *context = __nlCurrentContext;
+
        __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);
+       __nl_assert(!context->least_squares);
+
+       if (context->variable[row].locked);
+       else if (context->variable[col].locked) {
+               row = context->variable[row].index;
+               __nlRowColumnAppend(context->variable[col].a, row, value);
+       }
+       else {
+               __NLSparseMatrix* M  = &context->M;
+
+               row = context->variable[row].index;
+               col = context->variable[col].index;
+               
+               __nl_range_assert(row, 0, context->n - 1);
+               __nl_range_assert(col, 0, context->n - 1);
 
-       __nlSparseMatrixAdd(M, row, col, value);
+               __nlSparseMatrixAdd(M, row, col, value);
+       }
 }
 
 void nlRightHandSideAdd(NLuint index, NLfloat value)
 {
-       NLfloat* b = __nlCurrentContext->b;
+       __NLContext *context = __nlCurrentContext;
+       NLfloat* b = context->b;
 
        __nlCheckState(__NL_STATE_MATRIX);
-       __nl_range_assert(index, 0, __nlCurrentContext->n - 1);
-       __nl_assert(!__nlCurrentContext->least_squares);
+       __nl_assert(!context->least_squares);
 
-       b[index] += value;
+       if(!context->variable[index].locked) {
+               index = context->variable[index].index;
+               __nl_range_assert(index, 0, context->n - 1);
+
+               b[index] += value;
+       }
+}
+
+void nlRightHandSideSet(NLuint index, NLfloat value)
+{
+       __NLContext *context = __nlCurrentContext;
+       NLfloat* b = context->b;
+
+       __nlCheckState(__NL_STATE_MATRIX);
+       __nl_assert(!context->least_squares);
+
+       if(!context->variable[index].locked) {
+               index = context->variable[index].index;
+               __nl_range_assert(index, 0, context->n - 1);
+
+               b[index] = value;
+       }
 }
 
 void nlCoefficient(NLuint index, NLfloat value) {
@@ -1049,7 +1120,7 @@ static NLboolean __nlFactorize_SUPERLU(__NLContext *context, NLint *permutation)
        /* Cleanup */
 
        Destroy_SuperMatrix_Store(&At);
-       Destroy_SuperMatrix_Store(&AtP);
+       Destroy_CompCol_Permuted(&AtP);
 
        __NL_DELETE_ARRAY(etree);
        __NL_DELETE_ARRAY(xa);