00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <NL/nl_private.h>
00038 #include <NL/nl_matrix.h>
00039 #include <NL/nl_context.h>
00040 #include <NL/nl_iterative_solvers.h>
00041 #include <NL/nl_preconditioners.h>
00042 #include <NL/nl_superlu.h>
00043
00044 NLboolean nlInitExtension(char* extension) {
00045 #ifdef NL_USE_SUPERLU
00046 if(!strcmp(extension, "SUPERLU")) {
00047 return NL_TRUE ;
00048 }
00049 #endif
00050 return NL_FALSE ;
00051 }
00052
00053
00054
00055
00056
00057 void nlSolverParameterd(NLenum pname, NLdouble param) {
00058 nlCheckState(NL_STATE_INITIAL) ;
00059 switch(pname) {
00060 case NL_SOLVER: {
00061 nlCurrentContext->solver = (NLenum)param ;
00062 } break ;
00063 case NL_NB_VARIABLES: {
00064 nl_assert(param > 0) ;
00065 nlCurrentContext->nb_variables = (NLuint)param ;
00066 } break ;
00067 case NL_LEAST_SQUARES: {
00068 nlCurrentContext->least_squares = (NLboolean)param ;
00069 } break ;
00070 case NL_MAX_ITERATIONS: {
00071 nl_assert(param > 0) ;
00072 nlCurrentContext->max_iterations = (NLuint)param ;
00073 } break ;
00074 case NL_THRESHOLD: {
00075 nl_assert(param >= 0) ;
00076 nlCurrentContext->threshold = (NLdouble)param ;
00077 } break ;
00078 case NL_OMEGA: {
00079 nl_range_assert(param,1.0,2.0) ;
00080 nlCurrentContext->omega = (NLdouble)param ;
00081 } break ;
00082 case NL_SYMMETRIC: {
00083 nlCurrentContext->symmetric = (NLboolean)param ;
00084 }
00085 case NL_INNER_ITERATIONS: {
00086 nl_assert(param > 0) ;
00087 nlCurrentContext->inner_iterations = (NLuint)param ;
00088 } break ;
00089 case NL_PRECONDITIONER: {
00090 nlCurrentContext->preconditioner = (NLuint)param ;
00091 } break ;
00092 default: {
00093 nl_assert_not_reached ;
00094 } break ;
00095 }
00096 }
00097
00098 void nlSolverParameteri(NLenum pname, NLint param) {
00099 nlCheckState(NL_STATE_INITIAL) ;
00100 switch(pname) {
00101 case NL_SOLVER: {
00102 nlCurrentContext->solver = (NLenum)param ;
00103 } break ;
00104 case NL_NB_VARIABLES: {
00105 nl_assert(param > 0) ;
00106 nlCurrentContext->nb_variables = (NLuint)param ;
00107 } break ;
00108 case NL_LEAST_SQUARES: {
00109 nlCurrentContext->least_squares = (NLboolean)param ;
00110 } break ;
00111 case NL_MAX_ITERATIONS: {
00112 nl_assert(param > 0) ;
00113 nlCurrentContext->max_iterations = (NLuint)param ;
00114 } break ;
00115 case NL_THRESHOLD: {
00116 nl_assert(param >= 0) ;
00117 nlCurrentContext->threshold = (NLdouble)param ;
00118 } break ;
00119 case NL_OMEGA: {
00120 nl_range_assert(param,1,2) ;
00121 nlCurrentContext->omega = (NLdouble)param ;
00122 } break ;
00123 case NL_SYMMETRIC: {
00124 nlCurrentContext->symmetric = (NLboolean)param ;
00125 }
00126 case NL_INNER_ITERATIONS: {
00127 nl_assert(param > 0) ;
00128 nlCurrentContext->inner_iterations = (NLuint)param ;
00129 } break ;
00130 case NL_PRECONDITIONER: {
00131 nlCurrentContext->preconditioner = (NLuint)param ;
00132 } break ;
00133 default: {
00134 nl_assert_not_reached ;
00135 } break ;
00136 }
00137 }
00138
00139 void nlRowParameterd(NLenum pname, NLdouble param) {
00140 nlCheckState(NL_STATE_MATRIX) ;
00141 switch(pname) {
00142 case NL_RIGHT_HAND_SIDE: {
00143 nlCurrentContext->right_hand_side = param ;
00144 } break ;
00145 case NL_ROW_SCALING: {
00146 nlCurrentContext->row_scaling = param ;
00147 } break ;
00148 }
00149 }
00150
00151 void nlRowParameteri(NLenum pname, NLint param) {
00152 nlCheckState(NL_STATE_MATRIX) ;
00153 switch(pname) {
00154 case NL_RIGHT_HAND_SIDE: {
00155 nlCurrentContext->right_hand_side = (NLdouble)param ;
00156 } break ;
00157 case NL_ROW_SCALING: {
00158 nlCurrentContext->row_scaling = (NLdouble)param ;
00159 } break ;
00160 }
00161 }
00162
00163 void nlGetBooleanv(NLenum pname, NLboolean* params) {
00164 switch(pname) {
00165 case NL_LEAST_SQUARES: {
00166 *params = nlCurrentContext->least_squares ;
00167 } break ;
00168 case NL_SYMMETRIC: {
00169 *params = nlCurrentContext->symmetric ;
00170 } break ;
00171 default: {
00172 nl_assert_not_reached ;
00173 } break ;
00174 }
00175 }
00176
00177 void nlGetDoublev(NLenum pname, NLdouble* params) {
00178 switch(pname) {
00179 case NL_SOLVER: {
00180 *params = (NLdouble)(nlCurrentContext->solver) ;
00181 } break ;
00182 case NL_NB_VARIABLES: {
00183 *params = (NLdouble)(nlCurrentContext->nb_variables) ;
00184 } break ;
00185 case NL_LEAST_SQUARES: {
00186 *params = (NLdouble)(nlCurrentContext->least_squares) ;
00187 } break ;
00188 case NL_MAX_ITERATIONS: {
00189 *params = (NLdouble)(nlCurrentContext->max_iterations) ;
00190 } break ;
00191 case NL_THRESHOLD: {
00192 *params = (NLdouble)(nlCurrentContext->threshold) ;
00193 } break ;
00194 case NL_OMEGA: {
00195 *params = (NLdouble)(nlCurrentContext->omega) ;
00196 } break ;
00197 case NL_SYMMETRIC: {
00198 *params = (NLdouble)(nlCurrentContext->symmetric) ;
00199 } break ;
00200 case NL_USED_ITERATIONS: {
00201 *params = (NLdouble)(nlCurrentContext->used_iterations) ;
00202 } break ;
00203 case NL_ERROR: {
00204 *params = (NLdouble)(nlCurrentContext->error) ;
00205 } break ;
00206 case NL_ELAPSED_TIME: {
00207 *params = (NLdouble)(nlCurrentContext->elapsed_time) ;
00208 } break ;
00209 case NL_PRECONDITIONER: {
00210 *params = (NLdouble)(nlCurrentContext->preconditioner) ;
00211 } break ;
00212 default: {
00213 nl_assert_not_reached ;
00214 } break ;
00215 }
00216 }
00217
00218 void nlGetIntergerv(NLenum pname, NLint* params) {
00219 switch(pname) {
00220 case NL_SOLVER: {
00221 *params = (NLint)(nlCurrentContext->solver) ;
00222 } break ;
00223 case NL_NB_VARIABLES: {
00224 *params = (NLint)(nlCurrentContext->nb_variables) ;
00225 } break ;
00226 case NL_LEAST_SQUARES: {
00227 *params = (NLint)(nlCurrentContext->least_squares) ;
00228 } break ;
00229 case NL_MAX_ITERATIONS: {
00230 *params = (NLint)(nlCurrentContext->max_iterations) ;
00231 } break ;
00232 case NL_THRESHOLD: {
00233 *params = (NLint)(nlCurrentContext->threshold) ;
00234 } break ;
00235 case NL_OMEGA: {
00236 *params = (NLint)(nlCurrentContext->omega) ;
00237 } break ;
00238 case NL_SYMMETRIC: {
00239 *params = (NLint)(nlCurrentContext->symmetric) ;
00240 } break ;
00241 case NL_USED_ITERATIONS: {
00242 *params = (NLint)(nlCurrentContext->used_iterations) ;
00243 } break ;
00244 case NL_PRECONDITIONER: {
00245 *params = (NLint)(nlCurrentContext->preconditioner) ;
00246 } break ;
00247 default: {
00248 nl_assert_not_reached ;
00249 } break ;
00250 }
00251 }
00252
00253
00254
00255
00256 void nlEnable(NLenum pname) {
00257 switch(pname) {
00258 case NL_NORMALIZE_ROWS: {
00259 nl_assert(nlCurrentContext->state != NL_STATE_ROW) ;
00260 nlCurrentContext->normalize_rows = NL_TRUE ;
00261 } break ;
00262 default: {
00263 nl_assert_not_reached ;
00264 }
00265 }
00266 }
00267
00268 void nlDisable(NLenum pname) {
00269 switch(pname) {
00270 case NL_NORMALIZE_ROWS: {
00271 nl_assert(nlCurrentContext->state != NL_STATE_ROW) ;
00272 nlCurrentContext->normalize_rows = NL_FALSE ;
00273 } break ;
00274 default: {
00275 nl_assert_not_reached ;
00276 }
00277 }
00278 }
00279
00280 NLboolean nlIsEnabled(NLenum pname) {
00281 switch(pname) {
00282 case NL_NORMALIZE_ROWS: {
00283 return nlCurrentContext->normalize_rows ;
00284 } break ;
00285 default: {
00286 nl_assert_not_reached ;
00287 }
00288 }
00289 return NL_FALSE ;
00290 }
00291
00292
00293
00294
00295 void nlSetFunction(NLenum pname, void* param) {
00296 switch(pname) {
00297 case NL_FUNC_SOLVER:
00298 nlCurrentContext->solver_func = (NLSolverFunc)(param) ;
00299 break ;
00300 case NL_FUNC_MATRIX:
00301 nlCurrentContext->matrix_vector_prod = (NLMatrixFunc)(param) ;
00302 nlCurrentContext->solver = NL_SOLVER_USER ;
00303 break ;
00304 case NL_FUNC_PRECONDITIONER:
00305 nlCurrentContext->precond_vector_prod = (NLMatrixFunc)(param) ;
00306 nlCurrentContext->preconditioner = NL_PRECOND_USER ;
00307 break ;
00308 default:
00309 nl_assert_not_reached ;
00310 }
00311 }
00312
00313 void nlGetFunction(NLenum pname, void** param) {
00314 switch(pname) {
00315 case NL_FUNC_SOLVER:
00316 *param = nlCurrentContext->solver_func ;
00317 break ;
00318 case NL_FUNC_MATRIX:
00319 *param = nlCurrentContext->matrix_vector_prod ;
00320 break ;
00321 case NL_FUNC_PRECONDITIONER:
00322 *param = nlCurrentContext->precond_vector_prod ;
00323 break ;
00324 default:
00325 nl_assert_not_reached ;
00326 }
00327 }
00328
00329
00330
00331
00332 void nlSetVariable(NLuint index, NLdouble value) {
00333 nlCheckState(NL_STATE_SYSTEM) ;
00334 nl_debug_range_assert(index, 0, nlCurrentContext->nb_variables - 1) ;
00335 nlCurrentContext->variable[index].value = value ;
00336 }
00337
00338 NLdouble nlGetVariable(NLuint index) {
00339 nl_assert(nlCurrentContext->state != NL_STATE_INITIAL) ;
00340 nl_debug_range_assert(index, 0, nlCurrentContext->nb_variables - 1) ;
00341 return nlCurrentContext->variable[index].value ;
00342 }
00343
00344 void nlLockVariable(NLuint index) {
00345 nlCheckState(NL_STATE_SYSTEM) ;
00346 nl_debug_range_assert(index, 0, nlCurrentContext->nb_variables - 1) ;
00347 nlCurrentContext->variable[index].locked = NL_TRUE ;
00348 }
00349
00350 void nlUnlockVariable(NLuint index) {
00351 nlCheckState(NL_STATE_SYSTEM) ;
00352 nl_debug_range_assert(index, 0, nlCurrentContext->nb_variables - 1) ;
00353 nlCurrentContext->variable[index].locked = NL_FALSE ;
00354 }
00355
00356 NLboolean nlVariableIsLocked(NLuint index) {
00357 nl_assert(nlCurrentContext->state != NL_STATE_INITIAL) ;
00358 nl_debug_range_assert(index, 0, nlCurrentContext->nb_variables - 1) ;
00359 return nlCurrentContext->variable[index].locked ;
00360 }
00361
00362
00363
00364
00365 void nlVariablesToVector() {
00366 NLuint i ;
00367 nl_assert(nlCurrentContext->alloc_x) ;
00368 nl_assert(nlCurrentContext->alloc_variable) ;
00369 for(i=0; i<nlCurrentContext->nb_variables; i++) {
00370 NLVariable* v = &(nlCurrentContext->variable[i]) ;
00371 if(!v->locked) {
00372 nl_assert(v->index < nlCurrentContext->n) ;
00373 nlCurrentContext->x[v->index] = v->value ;
00374 }
00375 }
00376 }
00377
00378 void nlVectorToVariables() {
00379 NLuint i ;
00380 nl_assert(nlCurrentContext->alloc_x) ;
00381 nl_assert(nlCurrentContext->alloc_variable) ;
00382 for(i=0; i<nlCurrentContext->nb_variables; i++) {
00383 NLVariable* v = &(nlCurrentContext->variable[i]) ;
00384 if(!v->locked) {
00385 nl_assert(v->index < nlCurrentContext->n) ;
00386 v->value = nlCurrentContext->x[v->index] ;
00387 }
00388 }
00389 }
00390
00391
00392 void nlBeginSystem() {
00393 nlTransition(NL_STATE_INITIAL, NL_STATE_SYSTEM) ;
00394 nl_assert(nlCurrentContext->nb_variables > 0) ;
00395 nlCurrentContext->variable = NL_NEW_ARRAY(
00396 NLVariable, nlCurrentContext->nb_variables
00397 ) ;
00398 nlCurrentContext->alloc_variable = NL_TRUE ;
00399 }
00400
00401 void nlEndSystem() {
00402 nlTransition(NL_STATE_MATRIX_CONSTRUCTED, NL_STATE_SYSTEM_CONSTRUCTED) ;
00403 }
00404
00405 void nlBeginMatrix() {
00406 NLuint i ;
00407 NLuint n = 0 ;
00408 NLenum storage = NL_MATRIX_STORE_ROWS ;
00409
00410 nlTransition(NL_STATE_SYSTEM, NL_STATE_MATRIX) ;
00411
00412 for(i=0; i<nlCurrentContext->nb_variables; i++) {
00413 if(!nlCurrentContext->variable[i].locked) {
00414 nlCurrentContext->variable[i].index = n ;
00415 n++ ;
00416 } else {
00417 nlCurrentContext->variable[i].index = ~0 ;
00418 }
00419 }
00420
00421 nlCurrentContext->n = n ;
00422
00423
00424 if(nlCurrentContext->preconditioner == NL_PRECOND_SSOR) {
00425 storage = (storage | NL_MATRIX_STORE_COLUMNS) ;
00426 }
00427
00428
00429 if(nlCurrentContext->least_squares) {
00430 nlCurrentContext->symmetric = NL_TRUE ;
00431 }
00432
00433 if(nlCurrentContext->symmetric) {
00434 storage = (storage | NL_MATRIX_STORE_SYMMETRIC) ;
00435 }
00436
00437
00438 if(
00439 nlCurrentContext->solver == NL_SUPERLU_EXT ||
00440 nlCurrentContext->solver == NL_PERM_SUPERLU_EXT ||
00441 nlCurrentContext->solver == NL_SYMMETRIC_SUPERLU_EXT
00442 ) {
00443 storage = (storage & ~NL_SYMMETRIC) ;
00444 }
00445
00446 nlSparseMatrixConstruct(&nlCurrentContext->M, n, n, storage) ;
00447 nlCurrentContext->alloc_M = NL_TRUE ;
00448
00449 nlCurrentContext->x = NL_NEW_ARRAY(NLdouble, n) ;
00450 nlCurrentContext->alloc_x = NL_TRUE ;
00451
00452 nlCurrentContext->b = NL_NEW_ARRAY(NLdouble, n) ;
00453 nlCurrentContext->alloc_b = NL_TRUE ;
00454
00455 nlVariablesToVector() ;
00456
00457 nlRowColumnConstruct(&nlCurrentContext->af) ;
00458 nlCurrentContext->alloc_af = NL_TRUE ;
00459 nlRowColumnConstruct(&nlCurrentContext->al) ;
00460 nlCurrentContext->alloc_al = NL_TRUE ;
00461 nlRowColumnConstruct(&nlCurrentContext->xl) ;
00462 nlCurrentContext->alloc_xl = NL_TRUE ;
00463
00464 nlCurrentContext->current_row = 0 ;
00465 }
00466
00467 void nlEndMatrix() {
00468 nlTransition(NL_STATE_MATRIX, NL_STATE_MATRIX_CONSTRUCTED) ;
00469
00470 nlRowColumnDestroy(&nlCurrentContext->af) ;
00471 nlCurrentContext->alloc_af = NL_FALSE ;
00472 nlRowColumnDestroy(&nlCurrentContext->al) ;
00473 nlCurrentContext->alloc_al = NL_FALSE ;
00474 nlRowColumnDestroy(&nlCurrentContext->xl) ;
00475 nlCurrentContext->alloc_al = NL_FALSE ;
00476
00477 if(!nlCurrentContext->least_squares) {
00478 nl_assert(
00479 nlCurrentContext->current_row ==
00480 nlCurrentContext->n
00481 ) ;
00482 }
00483 }
00484
00485 void nlBeginRow() {
00486 nlTransition(NL_STATE_MATRIX, NL_STATE_ROW) ;
00487 nlRowColumnZero(&nlCurrentContext->af) ;
00488 nlRowColumnZero(&nlCurrentContext->al) ;
00489 nlRowColumnZero(&nlCurrentContext->xl) ;
00490 }
00491
00492 void nlScaleRow(NLdouble s) {
00493 NLRowColumn* af = &nlCurrentContext->af ;
00494 NLRowColumn* al = &nlCurrentContext->al ;
00495 NLuint nf = af->size ;
00496 NLuint nl = al->size ;
00497 NLuint i ;
00498 for(i=0; i<nf; i++) {
00499 af->coeff[i].value *= s ;
00500 }
00501 for(i=0; i<nl; i++) {
00502 al->coeff[i].value *= s ;
00503 }
00504 nlCurrentContext->right_hand_side *= s ;
00505 }
00506
00507 void nlNormalizeRow(NLdouble weight) {
00508 NLRowColumn* af = &nlCurrentContext->af ;
00509 NLRowColumn* al = &nlCurrentContext->al ;
00510 NLuint nf = af->size ;
00511 NLuint nl = al->size ;
00512 NLuint i ;
00513 NLdouble norm = 0.0 ;
00514 for(i=0; i<nf; i++) {
00515 norm += af->coeff[i].value * af->coeff[i].value ;
00516 }
00517 for(i=0; i<nl; i++) {
00518 norm += al->coeff[i].value * al->coeff[i].value ;
00519 }
00520 norm = sqrt(norm) ;
00521 nlScaleRow(weight / norm) ;
00522 }
00523
00524 void nlEndRow() {
00525 NLRowColumn* af = &nlCurrentContext->af ;
00526 NLRowColumn* al = &nlCurrentContext->al ;
00527 NLRowColumn* xl = &nlCurrentContext->xl ;
00528 NLSparseMatrix* M = &nlCurrentContext->M ;
00529 NLdouble* b = nlCurrentContext->b ;
00530 NLuint nf = af->size ;
00531 NLuint nl = al->size ;
00532 NLuint current_row = nlCurrentContext->current_row ;
00533 NLuint i ;
00534 NLuint j ;
00535 NLdouble S ;
00536 nlTransition(NL_STATE_ROW, NL_STATE_MATRIX) ;
00537
00538 if(nlCurrentContext->normalize_rows) {
00539 nlNormalizeRow(nlCurrentContext->row_scaling) ;
00540 } else {
00541 nlScaleRow(nlCurrentContext->row_scaling) ;
00542 }
00543
00544 if(nlCurrentContext->least_squares) {
00545 for(i=0; i<nf; i++) {
00546 for(j=0; j<nf; j++) {
00547 nlSparseMatrixAdd(
00548 M, af->coeff[i].index, af->coeff[j].index,
00549 af->coeff[i].value * af->coeff[j].value
00550 ) ;
00551 }
00552 }
00553 S = -nlCurrentContext->right_hand_side ;
00554 for(j=0; j<nl; j++) {
00555 S += al->coeff[j].value * xl->coeff[j].value ;
00556 }
00557 for(i=0; i<nf; i++) {
00558 b[ af->coeff[i].index ] -= af->coeff[i].value * S ;
00559 }
00560 } else {
00561 for(i=0; i<nf; i++) {
00562 nlSparseMatrixAdd(
00563 M, current_row, af->coeff[i].index, af->coeff[i].value
00564 ) ;
00565 }
00566 b[current_row] = -nlCurrentContext->right_hand_side ;
00567 for(i=0; i<nl; i++) {
00568 b[current_row] -= al->coeff[i].value * xl->coeff[i].value ;
00569 }
00570 }
00571 nlCurrentContext->current_row++ ;
00572 nlCurrentContext->right_hand_side = 0.0 ;
00573 nlCurrentContext->row_scaling = 1.0 ;
00574 }
00575
00576 void nlCoefficient(NLuint index, NLdouble value) {
00577 NLVariable* v = NULL ;
00578 nlCheckState(NL_STATE_ROW) ;
00579 nl_debug_range_assert(index, 0, nlCurrentContext->nb_variables - 1) ;
00580 v = &(nlCurrentContext->variable[index]) ;
00581 if(v->locked) {
00582 nlRowColumnAppend(&(nlCurrentContext->al), 0, value) ;
00583 nlRowColumnAppend(&(nlCurrentContext->xl), 0, v->value) ;
00584 } else {
00585 nlRowColumnAppend(&(nlCurrentContext->af), v->index, value) ;
00586 }
00587 }
00588
00589 void nlBegin(NLenum prim) {
00590 switch(prim) {
00591 case NL_SYSTEM: {
00592 nlBeginSystem() ;
00593 } break ;
00594 case NL_MATRIX: {
00595 nlBeginMatrix() ;
00596 } break ;
00597 case NL_ROW: {
00598 nlBeginRow() ;
00599 } break ;
00600 default: {
00601 nl_assert_not_reached ;
00602 }
00603 }
00604 }
00605
00606 void nlEnd(NLenum prim) {
00607 switch(prim) {
00608 case NL_SYSTEM: {
00609 nlEndSystem() ;
00610 } break ;
00611 case NL_MATRIX: {
00612 nlEndMatrix() ;
00613 } break ;
00614 case NL_ROW: {
00615 nlEndRow() ;
00616 } break ;
00617 default: {
00618 nl_assert_not_reached ;
00619 }
00620 }
00621 }
00622
00623
00624
00625
00626 NLboolean nlSolve() {
00627 NLboolean result ;
00628 nlCheckState(NL_STATE_SYSTEM_CONSTRUCTED) ;
00629 NLdouble start_time = nlCurrentTime() ;
00630 nlCurrentContext->elapsed_time = 0 ;
00631 result = nlCurrentContext->solver_func() ;
00632 nlVectorToVariables() ;
00633 nlCurrentContext->elapsed_time = nlCurrentTime() - start_time ;
00634 nlTransition(NL_STATE_SYSTEM_CONSTRUCTED, NL_STATE_SOLVED) ;
00635 return result ;
00636 }
00637
