[tuhopuu-devel] Fw: Help for tuhopuu's uv mapping
Bjornmose
tuhopuu-devel@blender.org
Tue, 11 May 2004 00:26:51 +0200
This is a multi-part message in MIME format.
------=_NextPart_000_005B_01C436EE.A78F01A0
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
hi brecht, luke and all.
just in case you did not get it
jow(bjorn the moose)
----- Original Message -----
From: "Bruno Levy" <Bruno.Levy@loria.fr>
To: <bjornmose@gmx.net>
Sent: Monday, May 10, 2004 12:22 PM
Subject: Help for tuhopuu's uv mapping
> Hello,
>
> I'm currently trying to extract a minimum kernel for numerical stuff
> (sorry for the delay, took
> longer than expected) here is a preliminary version (not tested yet).
I
> send it to you for early
> feedback, if you are interrested, I'll send a version to the mailing
list.
>
> Cheers,
> -- Bruno.
>
>
------------------------------------------------------------------------
--------
> /*
> * OpenNL: Numerical Library
> * Copyright (C) 2004 Bruno Levy
> *
> * This program is free software; you can redistribute it and/or
modify
> * it under the terms of the GNU General Public License as published
by
> * the Free Software Foundation; either version 2 of the License, or
> * (at your option) any later version.
> *
> * This program is distributed in the hope that it will be useful,
> * but WITHOUT ANY WARRANTY; without even the implied warranty of
> * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> * GNU General Public License for more details.
> *
> * You should have received a copy of the GNU General Public License
> * along with this program; if not, write to the Free Software
> * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
> *
> * If you modify this software, you should include a notice giving
the
> * name of the person performing the modification, the date of
modification,
> * and the reason for such modification.
> *
> * Contact: Bruno Levy
> *
> * levy@loria.fr
> *
> * ISA Project
> * LORIA, INRIA Lorraine,
> * Campus Scientifique, BP 239
> * 54506 VANDOEUVRE LES NANCY CEDEX
> * FRANCE
> *
> * Note that the GNU General Public License does not permit
incorporating
> * the Software into proprietary programs.
> */
>
> #ifndef __nl_h__
> #define __nl_h__
>
> #ifdef __cplusplus
> extern "C" {
> #endif
>
> #define NL_VERSION_0_0 1
>
> #define NLAPI
> #define NLAPIENTRY
>
> /*
> *
> * 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 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 void* NLContext ;
>
> /*
> *
> * Constants
> *
> */
>
> #define NL_FALSE 0x0
> #define NL_TRUE 0x1
>
> /* Primitives */
>
> #define NL_SYSTEM 0x0
> #define NL_MATRIX 0x1
> #define NL_ROW 0x2
>
>
> /* Solver Parameters */
>
> #define NL_SOLVER 0x100
> #define NL_NB_VARIABLES 0x101
> #define NL_LEAST_SQUARES 0x102
> #define NL_MAX_ITERATIONS 0x103
> #define NL_THRESHOLD 0x104
> #define NL_OMEGA 0x105
> #define NL_SYMMETRIC 0x106
> #define NL_USED_ITERATIONS 0x107
> #define NL_ERROR 0x108
>
> /* Solvers */
>
> #define NL_CG 0x200
> #define NL_JACOBI_CG 0x201
> #define NL_SSOR_CG 0x202
> #define NL_BICGSTAB 0x203
> #define NL_GMRES 0x204
> #define NL_SUPERLU_EXT 0x210
> #define NL_PERMSUPERLU_EXT 0x211
>
> /* Enable / Disable */
>
> #define NL_NORMALIZE_ROWS 0x300
>
> /* Row parameters */
>
> #define RIGHT_HAND_SIDE 0x400
> #define ROW_SCALING 0x401
>
> /*
> * Contexts
> */
> NLAPI NLContext NLAPIENTRY nlNewContext() ;
> NLAPI void NLAPIENTRY nlDeleteContext(NLContext context) ;
> NLAPI void NLAPIENTRY nlMakeCurrent(NLContext context) ;
> NLAPI NLContext NLAPIENTRY nlGetCurrent() ;
> NLAPI NLboolean NLAPIENTRY nlInitExtension(char* extension) ;
>
> /*
> * State set/get
> */
>
> NLAPI void NLAPIENTRY nlSolverParameterd(NLenum pname, NLdouble
param) ;
> NLAPI void NLAPIENTRY nlSolverParameteri(NLenum pname, NLint
param) ;
>
> NLAPI void NLAPIENTRY nlRowParameterd(NLenum pname, NLdouble
param) ;
> NLAPI void NLAPIENTRY nlRowParameteri(NLenum pname, NLint param) ;
>
> NLAPI void NLAPIENTRY nlGetBooleanv(NLenum pname, NLboolean*
params) ;
> NLAPI void NLAPIENTRY nlGetDoublev(NLenum pname, NLdouble* params)
;
> NLAPI void NLAPIENTRY nlGetIntergerv(NLenum pname, NLint* params)
;
>
> NLAPI void NLAPIENTRY nlEnable(NLenum pname) ;
> NLAPI void NLAPIENTRY nlDisable(NLenum pname) ;
> NLAPI NLboolean nlIsEnabled(NLenum pname) ;
>
> /*
> * Variables
> */
> NLAPI void NLAPIENTRY nlSetVariable(NLuint index, NLdouble value)
;
> NLAPI NLdouble NLAPIENTRY nlGetVariable(NLuint index) ;
> NLAPI void NLAPIENTRY nlLockVariable(NLuint index) ;
> NLAPI void NLAPIENTRY nlUnlockVariable(NLuint index) ;
> NLAPI NLboolean NLAPIENTRY nlVariableIsLocked(NLuint index) ;
>
> /*
> * Begin/End
> */
>
> NLAPI void NLAPIENTRY nlBegin(NLenum primitive) ;
> NLAPI void NLAPIENTRY nlEnd(NLenum primitive) ;
> NLAPI void NLAPIENTRY nlCoefficient(NLuint index, NLdouble value)
;
>
> /*
> * Solve
> */
>
> NLAPI NLboolean NLAPIENTRY nlSolve() ;
>
>
>
> #ifdef __cplusplus
> }
> #endif
>
> #endif
>
------------------------------------------------------------------------
--------
> /*
> * OpenNL: Numerical Library
> * Copyright (C) 2004 Bruno Levy
> *
> * This program is free software; you can redistribute it and/or
modify
> * it under the terms of the GNU General Public License as published
by
> * the Free Software Foundation; either version 2 of the License, or
> * (at your option) any later version.
> *
> * This program is distributed in the hope that it will be useful,
> * but WITHOUT ANY WARRANTY; without even the implied warranty of
> * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> * GNU General Public License for more details.
> *
> * You should have received a copy of the GNU General Public License
> * along with this program; if not, write to the Free Software
> * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
> *
> * If you modify this software, you should include a notice giving
the
> * name of the person performing the modification, the date of
modification,
> * and the reason for such modification.
> *
> * Contact: Bruno Levy
> *
> * levy@loria.fr
> *
> * ISA Project
> * LORIA, INRIA Lorraine,
> * Campus Scientifique, BP 239
> * 54506 VANDOEUVRE LES NANCY CEDEX
> * FRANCE
> *
> * Note that the GNU General Public License does not permit
incorporating
> * the Software into proprietary programs.
> */
>
> #include "nl.h"
>
> #include <stdlib.h>
> #include <stdio.h>
> #include <string.h>
> #include <math.h>
>
>
/***********************************************************************
*************/
> /* Assertions */
>
> static void __nl_assertion_failed(char* cond, 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, 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(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,min_val,max_val) { \
> 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__) ; \
> }
>
> #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
>
>
>
/***********************************************************************
*************/
> /* Error-reporting function */
>
> static void __nl_error(char* function, char* message) {
> fprintf(stderr, "OpenNL error in %s(): %s\n", function, message) ;
> }
>
>
>
/***********************************************************************
*************/
> /* 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((NB),sizeof(T)))
> #define __NL_RENEW_ARRAY(T,x,NB) (T*)(realloc(x,(NB)*sizeof(T)))
> #define __NL_DELETE(x) free(x); x = NULL
> #define __NL_DELETE_ARRAY(x) free(x); x = NULL
>
> #define __NL_CLEAR(T, x) memset(x, 0, sizeof(T))
> #define __NL_CLEAR_ARRAY(T,x,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 == 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 __nlRowColumnZero(__NLRowColumn* c) {
> c->size = 0 ;
> }
>
> 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<n; 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 __nlSparseMatrixZero( __NLSparseMatrix* M) {
> NLuint i ;
> if(M->storage & __NL_ROWS) {
> for(i=0; i<M->m; i++) {
> __nlRowColumnZero(&(M->row[i])) ;
> }
> }
> if(M->storage & __NL_COLUMNS) {
> for(i=0; i<M->n; i++) {
> __nlRowColumnZero(&(M->column[i])) ;
> }
> }
> __NL_CLEAR_ARRAY(NLdouble, M->diag, M->diag_size) ;
> }
>
> 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
> ) {
> __NL_CLEAR_ARRAY(NLdouble, y, A->m) ;
> NLuint n = A->n ;
> NLuint j,ii ;
> __NLRowColumn* Cj = NULL ;
> __NLCoeff* c = NULL ;
> 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 */
>
> 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) ;
> }
> }
> }
>
>
/***********************************************************************
*************/
> /* NLContext data structure */
>
> typedef struct {
> NLdouble value ;
> NLboolean locked ;
> NLuint index ;
> } __NLVariable ;
>
> #define __NL_STATE_INITIAL 0
> #define __NL_STATE_SYSTEM 1
> #define __NL_STATE_MATRIX 2
> #define __NL_STATE_ROW 3
> #define __NL_STATE_MATRIX_CONSTRUCTED 4
> #define __NL_STATE_SYSTEM_CONSTRUCTED 5
> #define __NL_STATE_SOLVED 6
>
> typedef struct {
> NLenum state ;
> __NLVariable* variable ;
> NLuint n ;
> __NLSparseMatrix M ;
> __NLRowColumn af ;
> __NLRowColumn al ;
> __NLRowColumn xl ;
> NLdouble* x ;
> NLdouble* b ;
> NLdouble right_hand_side ;
> NLdouble row_scaling ;
> NLenum solver ;
> NLuint nb_variables ;
> NLuint current_row ;
> NLboolean least_squares ;
> NLboolean symmetric ;
> NLuint max_iterations ;
> NLdouble threshold ;
> NLdouble omega ;
> NLboolean normalize_rows ;
> NLboolean alloc_M ;
> NLboolean alloc_af ;
> NLboolean alloc_al ;
> NLboolean alloc_xl ;
> NLboolean alloc_variable ;
> NLboolean alloc_x ;
> NLboolean alloc_b ;
> NLuint used_iterations ;
> NLdouble error ;
> } __NLContext ;
>
> static __NLContext* __nlCurrentContext = NULL ;
>
> NLContext nlNewContext() {
> __NLContext* result = __NL_NEW(__NLContext) ;
>
> result->state = __NL_STATE_INITIAL ;
> result->solver = NL_BICGSTAB ;
> result->max_iterations = 100 ;
> result->threshold = 1e-6 ;
> result->omega = 1.5 ;
> result->row_scaling = 1.0 ;
> return result ;
> }
>
> void nlDeleteContext(NLContext context_in) {
> __NLContext* context = (__NLContext*)(context_in) ;
> if(__nlCurrentContext == context) {
> __nlCurrentContext = NULL ;
> }
> if(context->alloc_M) {
> __nlSparseMatrixDestroy(&context->M) ;
> }
> if(context->alloc_af) {
> __nlRowColumnDestroy(&context->af) ;
> }
> if(context->alloc_al) {
> __nlRowColumnDestroy(&context->al) ;
> }
> if(context->alloc_xl) {
> __nlRowColumnDestroy(&context->xl) ;
> }
> if(context->alloc_variable) {
> __NL_DELETE_ARRAY(context->variable) ;
> }
> if(context->alloc_x) {
> __NL_DELETE_ARRAY(context->x) ;
> }
> if(context->alloc_b) {
> __NL_DELETE_ARRAY(context->b) ;
> }
>
> #ifdef NL_PARANOID
> __NL_CLEAR(__NLContext, context) ;
> #endif
> __NL_DELETE(context) ;
> }
>
> void nlMakeCurrent(NLContext context) {
> __nlCurrentContext = (__NLContext*)(context) ;
> }
>
> NLContext nlGetCurrent() {
> return __nlCurrentContext ;
> }
>
> NLboolean nlInitExtension(char* extension) {
> return NL_FALSE ;
> }
>
> void __nlCheckState(NLenum state) {
> __nl_assert(__nlCurrentContext->state == state) ;
> }
>
> void __nlTransition(NLenum from_state, NLenum to_state) {
> __nlCheckState(from_state) ;
> __nlCurrentContext->state = to_state ;
> }
>
>
/***********************************************************************
*************/
> /* Get/Set parameters */
>
> void nlSolverParameterd(NLenum pname, NLdouble param) {
> __nlCheckState(__NL_STATE_INITIAL) ;
> switch(pname) {
> case NL_SOLVER: {
> __nlCurrentContext->solver = (NLenum)param ;
> } break ;
> case NL_NB_VARIABLES: {
> __nl_assert(param > 0) ;
> __nlCurrentContext->nb_variables = (NLuint)param ;
> } break ;
> case NL_LEAST_SQUARES: {
> __nlCurrentContext->least_squares = (NLboolean)param ;
> } break ;
> case NL_MAX_ITERATIONS: {
> __nl_assert(param > 0) ;
> __nlCurrentContext->max_iterations = (NLuint)param ;
> } break ;
> case NL_THRESHOLD: {
> __nl_assert(param >= 0) ;
> __nlCurrentContext->threshold = (NLdouble)param ;
> } break ;
> case NL_OMEGA: {
> __nl_range_assert(param,1.0,2.0) ;
> __nlCurrentContext->omega = (NLdouble)param ;
> } break ;
> case NL_SYMMETRIC: {
> __nlCurrentContext->symmetric = (NLboolean)param ;
> }
> default: {
> __nl_assert_not_reached ;
> } break ;
> }
> }
>
> void nlSolverParameteri(NLenum pname, NLint param) {
> __nlCheckState(__NL_STATE_INITIAL) ;
> switch(pname) {
> case NL_SOLVER: {
> __nlCurrentContext->solver = (NLenum)param ;
> } break ;
> case NL_NB_VARIABLES: {
> __nl_assert(param > 0) ;
> __nlCurrentContext->nb_variables = (NLuint)param ;
> } break ;
> case NL_LEAST_SQUARES: {
> __nlCurrentContext->least_squares = (NLboolean)param ;
> } break ;
> case NL_MAX_ITERATIONS: {
> __nl_assert(param > 0) ;
> __nlCurrentContext->max_iterations = (NLuint)param ;
> } break ;
> case NL_THRESHOLD: {
> __nl_assert(param >= 0) ;
> __nlCurrentContext->threshold = (NLdouble)param ;
> } break ;
> case NL_OMEGA: {
> __nl_range_assert(param,1,2) ;
> __nlCurrentContext->omega = (NLdouble)param ;
> } break ;
> case NL_SYMMETRIC: {
> __nlCurrentContext->symmetric = (NLboolean)param ;
> }
> default: {
> __nl_assert_not_reached ;
> } break ;
> }
> }
>
> void nlRowParameterd(NLenum pname, NLdouble param) {
> __nlCheckState(__NL_STATE_MATRIX) ;
> switch(pname) {
> case RIGHT_HAND_SIDE: {
> __nlCurrentContext->right_hand_side = param ;
> } break ;
> case ROW_SCALING: {
> __nlCurrentContext->row_scaling = param ;
> } break ;
> }
> }
>
> void nlRowParameteri(NLenum pname, NLint param) {
> __nlCheckState(__NL_STATE_MATRIX) ;
> switch(pname) {
> case RIGHT_HAND_SIDE: {
> __nlCurrentContext->right_hand_side = (NLdouble)param ;
> } break ;
> case ROW_SCALING: {
> __nlCurrentContext->row_scaling = (NLdouble)param ;
> } break ;
> }
> }
>
> void nlGetBooleanv(NLenum pname, NLboolean* params) {
> switch(pname) {
> case NL_LEAST_SQUARES: {
> *params = __nlCurrentContext->least_squares ;
> } break ;
> case NL_SYMMETRIC: {
> *params = __nlCurrentContext->symmetric ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> } break ;
> }
> }
>
> void nlGetDoublev(NLenum pname, NLdouble* params) {
> switch(pname) {
> case NL_SOLVER: {
> *params = (NLdouble)(__nlCurrentContext->solver) ;
> } break ;
> case NL_NB_VARIABLES: {
> *params = (NLdouble)(__nlCurrentContext->nb_variables) ;
> } break ;
> case NL_LEAST_SQUARES: {
> *params = (NLdouble)(__nlCurrentContext->least_squares) ;
> } break ;
> case NL_MAX_ITERATIONS: {
> *params = (NLdouble)(__nlCurrentContext->max_iterations) ;
> } break ;
> case NL_THRESHOLD: {
> *params = (NLdouble)(__nlCurrentContext->threshold) ;
> } break ;
> case NL_OMEGA: {
> *params = (NLdouble)(__nlCurrentContext->omega) ;
> } break ;
> case NL_SYMMETRIC: {
> *params = (NLdouble)(__nlCurrentContext->symmetric) ;
> } break ;
> case NL_USED_ITERATIONS: {
> *params = (NLdouble)(__nlCurrentContext->used_iterations) ;
> } break ;
> case NL_ERROR: {
> *params = (NLdouble)(__nlCurrentContext->error) ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> } break ;
> }
> }
>
> void nlGetIntergerv(NLenum pname, NLint* params) {
> switch(pname) {
> case NL_SOLVER: {
> *params = (NLint)(__nlCurrentContext->solver) ;
> } break ;
> case NL_NB_VARIABLES: {
> *params = (NLint)(__nlCurrentContext->nb_variables) ;
> } break ;
> case NL_LEAST_SQUARES: {
> *params = (NLint)(__nlCurrentContext->least_squares) ;
> } break ;
> case NL_MAX_ITERATIONS: {
> *params = (NLint)(__nlCurrentContext->max_iterations) ;
> } break ;
> case NL_THRESHOLD: {
> *params = (NLint)(__nlCurrentContext->threshold) ;
> } break ;
> case NL_OMEGA: {
> *params = (NLint)(__nlCurrentContext->omega) ;
> } break ;
> case NL_SYMMETRIC: {
> *params = (NLint)(__nlCurrentContext->symmetric) ;
> } break ;
> case NL_USED_ITERATIONS: {
> *params = (NLint)(__nlCurrentContext->used_iterations) ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> } break ;
> }
> }
>
>
/***********************************************************************
*************/
> /* Enable / Disable */
>
> void nlEnable(NLenum pname) {
> switch(pname) {
> case NL_NORMALIZE_ROWS: {
> __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW) ;
> __nlCurrentContext->normalize_rows = NL_TRUE ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> }
> }
> }
>
> void nlDisable(NLenum pname) {
> switch(pname) {
> case NL_NORMALIZE_ROWS: {
> __nl_assert(__nlCurrentContext->state != __NL_STATE_ROW) ;
> __nlCurrentContext->normalize_rows = NL_FALSE ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> }
> }
> }
>
> NLboolean nlIsEnabled(NLenum pname) {
> switch(pname) {
> case NL_NORMALIZE_ROWS: {
> return __nlCurrentContext->normalize_rows ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> }
> }
> return NL_FALSE ;
> }
>
>
/***********************************************************************
*************/
> /* Get/Set Lock/Unlock variables */
>
> void nlSetVariable(NLuint index, NLdouble value) {
> __nlCheckState(__NL_STATE_SYSTEM) ;
> __nl_parano_range_assert(index, 0,
__nlCurrentContext->nb_variables - 1) ;
> __nlCurrentContext->variable[index].value = value ;
> }
>
> NLdouble 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 ;
> }
>
> 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 nlUnlockVariable(NLuint index) {
> __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 ;
> }
>
>
/***********************************************************************
*************/
> /* System construction */
>
> void __nlVariablesToVector() {
> __nl_assert(__nlCurrentContext->alloc_x) ;
> __nl_assert(__nlCurrentContext->alloc_variable) ;
> NLuint i ;
> for(i=0; i<__nlCurrentContext->nb_variables; i++) {
> __NLVariable* v = &(__nlCurrentContext->variable[i]) ;
> if(!v->locked) {
> __nl_assert(v->index < __nlCurrentContext->n) ;
> __nlCurrentContext->x[v->index] = v->value ;
> }
> }
> }
>
> void __nlVectorToVariables() {
> __nl_assert(__nlCurrentContext->alloc_x) ;
> __nl_assert(__nlCurrentContext->alloc_variable) ;
> NLuint i ;
> for(i=0; i<__nlCurrentContext->nb_variables; i++) {
> __NLVariable* v = &(__nlCurrentContext->variable[i]) ;
> if(!v->locked) {
> __nl_assert(v->index < __nlCurrentContext->n) ;
> v->value = __nlCurrentContext->x[v->index] ;
> }
> }
> }
>
>
> void __nlBeginSystem() {
> __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 ;
> }
>
> void __nlEndSystem() {
> __nlTransition(__NL_STATE_MATRIX_CONSTRUCTED,
__NL_STATE_SYSTEM_CONSTRUCTED) ;
> }
>
> void __nlBeginMatrix() {
> NLuint i ;
> NLuint n = 0 ;
> NLenum storage = __NL_ROWS ;
>
> __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++ ;
> } else {
> __nlCurrentContext->variable[i].index = ~0 ;
> }
> }
>
> __nlCurrentContext->n = n ;
>
> /* SSOR preconditioner requires rows and columns */
> if(__nlCurrentContext->solver == NL_SSOR_CG) {
> storage = (storage | __NL_COLUMNS) ;
> }
>
> /* a least squares problem results in a symmetric matrix */
> if(__nlCurrentContext->least_squares) {
> __nlCurrentContext->symmetric = NL_TRUE ;
> }
>
> if(__nlCurrentContext->symmetric) {
> storage = (storage | __NL_SYMMETRIC) ;
> }
>
> /* SuperLU storage does not support symmetric storage */
> if(
> __nlCurrentContext->solver == NL_SUPERLU_EXT ||
> __nlCurrentContext->solver == NL_PERMSUPERLU_EXT
> ) {
> storage = (storage & ~__NL_SYMMETRIC) ;
> }
>
> __nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage) ;
> __nlCurrentContext->alloc_M = NL_TRUE ;
>
> __nlCurrentContext->x = __NL_NEW_ARRAY(NLdouble, n) ;
> __nlCurrentContext->alloc_x = NL_TRUE ;
>
> __nlCurrentContext->b = __NL_NEW_ARRAY(NLdouble, n) ;
> __nlCurrentContext->alloc_b = NL_TRUE ;
>
> __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 ;
>
> __nlCurrentContext->current_row = 0 ;
> }
>
> void __nlEndMatrix() {
> __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 ;
>
> if(!__nlCurrentContext->least_squares) {
> __nl_assert(
> __nlCurrentContext->current_row ==
> __nlCurrentContext->n
> ) ;
> }
> }
>
> void __nlBeginRow() {
> __nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW) ;
> __nlRowColumnZero(&__nlCurrentContext->af) ;
> __nlRowColumnZero(&__nlCurrentContext->al) ;
> }
>
> void __nlScaleRow(NLdouble s) {
> __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 ;
> }
> for(i=0; i<nl; i++) {
> al->coeff[i].value *= s ;
> }
> __nlCurrentContext->right_hand_side *= s ;
> }
>
> void __nlNormalizeRow(NLdouble weight) {
> __NLRowColumn* af = &__nlCurrentContext->af ;
> __NLRowColumn* al = &__nlCurrentContext->al ;
> NLuint nf = af->size ;
> NLuint nl = al->size ;
> NLuint i ;
> NLdouble norm = 0.0 ;
> for(i=0; i<nf; i++) {
> 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 = sqrt(norm) ;
> __nlScaleRow(weight / norm) ;
> }
>
> void __nlEndRow() {
> __NLRowColumn* af = &__nlCurrentContext->af ;
> __NLRowColumn* al = &__nlCurrentContext->al ;
> __NLRowColumn* xl = &__nlCurrentContext->xl ;
> __NLSparseMatrix* M = &__nlCurrentContext->M ;
> NLdouble* b = __nlCurrentContext->b ;
> NLuint nf = af->size ;
> NLuint nl = al->size ;
> NLuint current_row = __nlCurrentContext->current_row ;
> NLuint i ;
> NLuint j ;
> NLdouble S ;
> __nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX) ;
>
> if(__nlCurrentContext->normalize_rows) {
> __nlNormalizeRow(__nlCurrentContext->row_scaling) ;
> } else {
> __nlScaleRow(__nlCurrentContext->row_scaling) ;
> }
>
> if(__nlCurrentContext->least_squares) {
> for(i=0; i<nf; i++) {
> for(j=0; j<nf; j++) {
> __nlSparseMatrixAdd(
> M, af->coeff[i].index, af->coeff[j].index,
> af->coeff[i].value * af->coeff[j].value
> ) ;
> }
> }
> S = __nlCurrentContext->right_hand_side ;
> for(j=0; j<nl; j++) {
> 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 ;
> }
> } 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 ;
> for(i=0; i<nl; i++) {
> 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 ;
> }
>
> void nlCoefficient(NLuint index, NLdouble value) {
> __nlCheckState(__NL_STATE_ROW) ;
> __nl_range_assert(index, 0, __nlCurrentContext->nb_variables - 1)
;
> __NLVariable* v = &(__nlCurrentContext->variable[index]) ;
> if(v->locked) {
> __nlRowColumnAppend(&(__nlCurrentContext->al), 0, value) ;
> __nlRowColumnAppend(&(__nlCurrentContext->xl), 0, v->value) ;
> } else {
> __nlRowColumnAppend(&(__nlCurrentContext->af), v->index,
value) ;
> }
> }
>
> void nlBegin(NLenum prim) {
> switch(prim) {
> case NL_SYSTEM: {
> __nlBeginSystem() ;
> } break ;
> case NL_MATRIX: {
> __nlBeginMatrix() ;
> } break ;
> case NL_ROW: {
> __nlBeginRow() ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> }
> }
> }
>
> void nlEnd(NLenum prim) {
> switch(prim) {
> case NL_SYSTEM: {
> __nlEndSystem() ;
> } break ;
> case NL_MATRIX: {
> __nlEndMatrix() ;
> } break ;
> case NL_ROW: {
> __nlEndRow() ;
> } break ;
> default: {
> __nl_assert_not_reached ;
> }
> }
> }
>
>
/***********************************************************************
*************/
> /* BLAS routines */
>
> /*
> * daxpy
> * ddot
> * dscal
> * dnrm2
> * dcopy
> */
>
> #ifndef NL_USE_BLAS
>
> typedef NLint integer ;
> typedef NLdouble doublereal ;
>
> /* Subroutine */ int daxpy_(integer *n, doublereal *da, doublereal
*dx,
> integer *incx, doublereal *dy, integer *incy)
> {
>
>
> /* System generated locals */
> integer i__1;
>
> /* Local variables */
> static integer i, m, ix, iy, mp1;
>
>
> /* constant times a vector plus a vector.
> uses unrolled loops for increments equal to one.
> jack dongarra, linpack, 3/11/78.
> modified 12/3/93, array(1) declarations changed to array(*)
>
>
>
> Parameter adjustments
> Function Body */
> #define DY(I) dy[(I)-1]
> #define DX(I) dx[(I)-1]
>
>
> if (*n <= 0) {
> return 0;
> }
> if (*da == 0.) {
> return 0;
> }
> if (*incx == 1 && *incy == 1) {
> goto L20;
> }
>
> /* code for unequal increments or equal increments
> not equal to 1 */
>
> ix = 1;
> iy = 1;
> if (*incx < 0) {
> ix = (-(*n) + 1) * *incx + 1;
> }
> if (*incy < 0) {
> iy = (-(*n) + 1) * *incy + 1;
> }
> i__1 = *n;
> for (i = 1; i <= *n; ++i) {
> DY(iy) += *da * DX(ix);
> ix += *incx;
> iy += *incy;
> /* L10: */
> }
> return 0;
>
> /* code for both increments equal to 1
>
>
> clean-up loop */
>
> L20:
> m = *n % 4;
> if (m == 0) {
> goto L40;
> }
> i__1 = m;
> for (i = 1; i <= m; ++i) {
> DY(i) += *da * DX(i);
> /* L30: */
> }
> if (*n < 4) {
> return 0;
> }
> L40:
> mp1 = m + 1;
> i__1 = *n;
> for (i = mp1; i <= *n; i += 4) {
> DY(i) += *da * DX(i);
> DY(i + 1) += *da * DX(i + 1);
> DY(i + 2) += *da * DX(i + 2);
> DY(i + 3) += *da * DX(i + 3);
> /* L50: */
> }
> return 0;
> } /* daxpy_ */
> #undef DY
> #undef DX
>
>
> doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal
*dy,
> integer *incy)
> {
>
> /* System generated locals */
> integer i__1;
> doublereal ret_val;
>
> /* Local variables */
> static integer i, m;
> static doublereal dtemp;
> static integer ix, iy, mp1;
>
>
> /* forms the dot product of two vectors.
> uses unrolled loops for increments equal to one.
> jack dongarra, linpack, 3/11/78.
> modified 12/3/93, array(1) declarations changed to array(*)
>
>
>
> Parameter adjustments
> Function Body */
> #define DY(I) dy[(I)-1]
> #define DX(I) dx[(I)-1]
>
> ret_val = 0.;
> dtemp = 0.;
> if (*n <= 0) {
> return ret_val;
> }
> if (*incx == 1 && *incy == 1) {
> goto L20;
> }
>
> /* code for unequal increments or equal increments
> not equal to 1 */
>
> ix = 1;
> iy = 1;
> if (*incx < 0) {
> ix = (-(*n) + 1) * *incx + 1;
> }
> if (*incy < 0) {
> iy = (-(*n) + 1) * *incy + 1;
> }
> i__1 = *n;
> for (i = 1; i <= *n; ++i) {
> dtemp += DX(ix) * DY(iy);
> ix += *incx;
> iy += *incy;
> /* L10: */
> }
> ret_val = dtemp;
> return ret_val;
>
> /* code for both increments equal to 1
>
>
> clean-up loop */
>
> L20:
> m = *n % 5;
> if (m == 0) {
> goto L40;
> }
> i__1 = m;
> for (i = 1; i <= m; ++i) {
> dtemp += DX(i) * DY(i);
> /* L30: */
> }
> if (*n < 5) {
> goto L60;
> }
> L40:
> mp1 = m + 1;
> i__1 = *n;
> for (i = mp1; i <= *n; i += 5) {
> dtemp = dtemp + DX(i) * DY(i) + DX(i + 1) * DY(i + 1) + DX(i + 2) *
> DY(i + 2) + DX(i + 3) * DY(i + 3) + DX(i + 4) * DY(i + 4);
> /* L50: */
> }
> L60:
> ret_val = dtemp;
> return ret_val;
> } /* ddot_ */
> #undef DY
> #undef DX
>
> /* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal
*dx,
> integer *incx)
> {
>
>
> /* System generated locals */
> integer i__1, i__2;
>
> /* Local variables */
> static integer i, m, nincx, mp1;
>
>
> /* scales a vector by a constant.
> uses unrolled loops for increment equal to one.
> jack dongarra, linpack, 3/11/78.
> modified 3/93 to return if incx .le. 0.
> modified 12/3/93, array(1) declarations changed to array(*)
>
>
>
> Parameter adjustments
> Function Body */
> #ifdef DX
> #undef DX
> #endif
> #define DX(I) dx[(I)-1]
>
>
> if (*n <= 0 || *incx <= 0) {
> return 0;
> }
> if (*incx == 1) {
> goto L20;
> }
>
> /* code for increment not equal to 1 */
>
> nincx = *n * *incx;
> i__1 = nincx;
> i__2 = *incx;
> for (i = 1; *incx < 0 ? i >= nincx : i <= nincx; i += *incx) {
> DX(i) = *da * DX(i);
> /* L10: */
> }
> return 0;
>
> /* code for increment equal to 1
>
>
> clean-up loop */
>
> L20:
> m = *n % 5;
> if (m == 0) {
> goto L40;
> }
> i__2 = m;
> for (i = 1; i <= m; ++i) {
> DX(i) = *da * DX(i);
> /* L30: */
> }
> if (*n < 5) {
> return 0;
> }
> L40:
> mp1 = m + 1;
> i__2 = *n;
> for (i = mp1; i <= *n; i += 5) {
> DX(i) = *da * DX(i);
> DX(i + 1) = *da * DX(i + 1);
> DX(i + 2) = *da * DX(i + 2);
> DX(i + 3) = *da * DX(i + 3);
> DX(i + 4) = *da * DX(i + 4);
> /* L50: */
> }
> return 0;
> } /* dscal_ */
> #undef DX
>
> doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
> {
>
>
> /* System generated locals */
> integer i__1, i__2;
> doublereal ret_val, d__1;
>
> /* Builtin functions */
> double sqrt(doublereal);
>
> /* Local variables */
> static doublereal norm, scale, absxi;
> static integer ix;
> static doublereal ssq;
>
>
> /* DNRM2 returns the euclidean norm of a vector via the function
> name, so that
>
> DNRM2 := sqrt( x'*x )
>
>
>
> -- This version written on 25-October-1982.
> Modified on 14-October-1993 to inline the call to DLASSQ.
> Sven Hammarling, Nag Ltd.
>
>
>
> Parameter adjustments
> Function Body */
> #ifdef X
> #undef X
> #endif
> #define X(I) x[(I)-1]
>
>
> if (*n < 1 || *incx < 1) {
> norm = 0.;
> } else if (*n == 1) {
> norm = abs(X(1));
> } else {
> scale = 0.;
> ssq = 1.;
> /* The following loop is equivalent to this call to the LAPACK
>
> auxiliary routine:
> CALL DLASSQ( N, X, INCX, SCALE, SSQ ) */
>
> i__1 = (*n - 1) * *incx + 1;
> i__2 = *incx;
> for (ix = 1; *incx < 0 ? ix >= (*n-1)**incx+1 : ix <= (*n-1)**incx+1;
ix += *incx) {
> if (X(ix) != 0.) {
> absxi = (d__1 = X(ix), abs(d__1));
> if (scale < absxi) {
> /* Computing 2nd power */
> d__1 = scale / absxi;
> ssq = ssq * (d__1 * d__1) + 1.;
> scale = absxi;
> } else {
> /* Computing 2nd power */
> d__1 = absxi / scale;
> ssq += d__1 * d__1;
> }
> }
> /* L10: */
> }
> norm = scale * sqrt(ssq);
> }
>
> ret_val = norm;
> return ret_val;
>
> /* End of DNRM2. */
>
> } /* dnrm2_ */
> #undef X
>
> /* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx,
> doublereal *dy, integer *incy)
> {
>
> /* System generated locals */
> integer i__1;
>
> /* Local variables */
> static integer i, m, ix, iy, mp1;
>
>
> /* copies a vector, x, to a vector, y.
> uses unrolled loops for increments equal to one.
> jack dongarra, linpack, 3/11/78.
> modified 12/3/93, array(1) declarations changed to array(*)
>
>
>
> Parameter adjustments
> Function Body */
> #define DY(I) dy[(I)-1]
> #define DX(I) dx[(I)-1]
>
>
> if (*n <= 0) {
> return 0;
> }
> if (*incx == 1 && *incy == 1) {
> goto L20;
> }
>
> /* code for unequal increments or equal increments
> not equal to 1 */
>
> ix = 1;
> iy = 1;
> if (*incx < 0) {
> ix = (-(*n) + 1) * *incx + 1;
> }
> if (*incy < 0) {
> iy = (-(*n) + 1) * *incy + 1;
> }
> i__1 = *n;
> for (i = 1; i <= *n; ++i) {
> DY(iy) = DX(ix);
> ix += *incx;
> iy += *incy;
> /* L10: */
> }
> return 0;
>
> /* code for both increments equal to 1
>
>
> clean-up loop */
>
> L20:
> m = *n % 7;
> if (m == 0) {
> goto L40;
> }
> i__1 = m;
> for (i = 1; i <= m; ++i) {
> DY(i) = DX(i);
> /* L30: */
> }
> if (*n < 7) {
> return 0;
> }
> L40:
> mp1 = m + 1;
> i__1 = *n;
> for (i = mp1; i <= *n; i += 7) {
> DY(i) = DX(i);
> DY(i + 1) = DX(i + 1);
> DY(i + 2) = DX(i + 2);
> DY(i + 3) = DX(i + 3);
> DY(i + 4) = DX(i + 4);
> DY(i + 5) = DX(i + 5);
> DY(i + 6) = DX(i + 6);
> /* L50: */
> }
> return 0;
> } /* dcopy_ */
>
> #undef DX
> #undef DY
>
> #else
>
> extern void daxpy_(
> int *n, double *alpha, double *x,
> int *incx, double *y, int *incy
> ) ;
>
> extern double ddot_(
> int *n, double *x, int *incx, double *y,
> int *incy
> ) ;
>
> extern double dnrm2_( int *n, double *x, int *incx ) ;
>
> extern int dcopy_(int* n, double* dx, int* incx, double* dy, int*
incy) ;
>
> extern void dscal_(int* n, double* alpha, double *x, int* incx) ;
>
> #endif
>
>
/***********************************************************************
************/
> /* C wrappers for BLAS routines */
>
> /* x <- a*x */
> static void dscal( int n, double alpha, double *x, int incx ) {
> dscal_(&n,&alpha,x,&incx);
> }
>
> /* y <- x */
> static void dcopy(
> int n, double *x, int incx, double *y, int incy
> ) {
> dcopy_(&n,x,&incx,y,&incy);
> }
>
> /* y <- a*x+y */
> static void daxpy(
> int n, double alpha, double *x, int incx, double *y,
> int incy
> ) {
> daxpy_(&n,&alpha,x,&incx,y,&incy);
> }
>
> /* returns x^T*y */
> static double ddot(
> int n, double *x, int incx, double *y, int incy
> ) {
> return ddot_(&n,x,&incx,y,&incy);
> }
>
> /** returns |x|_2 */
> static double dnrm2( int n, double *x, int incx ) {
> double d=0.;
> while ( n-- ) {
> d+=(*x)*(*x);
> x+=incx ;
> }
> return sqrt(d);
> // return dnrm2_(&n,x,&incx);
> }
>
>
>
/***********************************************************************
*************/
> /* Solvers */
>
> static NLuint __nlSolve_CG(
> __NLSparseMatrix* A,
> NLdouble* b, NLdouble* x,
> NLdouble eps, NLuint max_iter
> ) {
> NLint N = A->n ;
> NLdouble *g = __NL_NEW_ARRAY(NLdouble, N) ;
> NLdouble *r = __NL_NEW_ARRAY(NLdouble, N) ;
> NLdouble *p = __NL_NEW_ARRAY(NLdouble, N) ;
> NLuint its=0;
> NLdouble t, tau, sig, rho, gam;
> NLdouble err=eps*eps*ddot(N,b,1,b,1);
>
> __nlSparseMatrixMult(A,x,g);
> daxpy(N,-1.,b,1,g,1);
> dscal(N,-1.,g,1);
> dcopy(N,g,1,r,1);
> while ( ddot(N,g,1,g,1)>err && its < max_iter) {
> __nlSparseMatrixMult(A,r,p);
> rho=ddot(N,p,1,p,1);
> sig=ddot(N,r,1,p,1);
> tau=ddot(N,g,1,r,1);
> t=tau/sig;
> daxpy(N,t,r,1,x,1);
> daxpy(N,-t,p,1,g,1);
> gam=(t*t*rho-tau)/tau;
> dscal(N,gam,r,1);
> daxpy(N,1.,g,1,r,1);
> ++its;
> }
> __NL_DELETE_ARRAY(g) ;
> __NL_DELETE_ARRAY(r) ;
> __NL_DELETE_ARRAY(p) ;
> return its;
> }
>
> NLboolean nlSolve() {
> __nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SOLVED) ;
> switch(__nlCurrentContext->solver) {
> case NL_CG: {
> __nlCurrentContext->used_iterations =
> __nlSolve_CG(
> &(__nlCurrentContext->M),
> __nlCurrentContext->b,
> __nlCurrentContext->x,
> __nlCurrentContext->threshold,
> __nlCurrentContext->max_iterations
> ) ;
> } break ;
> case NL_JACOBI_CG: {
> __nl_assert_not_reached ;
> } break ;
> case NL_SSOR_CG: {
> __nl_assert_not_reached ;
> } break ;
> case NL_BICGSTAB: {
> __nl_assert_not_reached ;
> } break ;
> case NL_GMRES: {
> __nl_assert_not_reached ;
> } break ;
> case NL_SUPERLU_EXT: {
> __nl_assert_not_reached ;
> } break ;
> case NL_PERMSUPERLU_EXT: {
> __nl_assert_not_reached ;
> } break ;
> default:
> __nl_assert_not_reached ;
> }
> __nlVectorToVariables() ;
> return NL_TRUE ;
> }
>
>
------=_NextPart_000_005B_01C436EE.A78F01A0
Content-Type: application/octet-stream;
name="nl.h"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="nl.h"
/*
* OpenNL: Numerical Library
* Copyright (C) 2004 Bruno Levy
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* If you modify this software, you should include a notice giving the
* name of the person performing the modification, the date of =
modification,
* and the reason for such modification.
*
* Contact: Bruno Levy
*
* levy@loria.fr
*
* ISA Project
* LORIA, INRIA Lorraine,=20
* Campus Scientifique, BP 239
* 54506 VANDOEUVRE LES NANCY CEDEX=20
* FRANCE
*
* Note that the GNU General Public License does not permit =
incorporating
* the Software into proprietary programs.=20
*/
#ifndef __nl_h__
#define __nl_h__
#ifdef __cplusplus
extern "C" {
#endif
#define NL_VERSION_0_0 1
#define NLAPI
#define NLAPIENTRY
/*
*
* 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 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 void* NLContext ;
/*
*
* Constants
*
*/
#define NL_FALSE 0x0
#define NL_TRUE 0x1
/* Primitives */
#define NL_SYSTEM 0x0
#define NL_MATRIX 0x1
#define NL_ROW 0x2
/* Solver Parameters */
#define NL_SOLVER 0x100
#define NL_NB_VARIABLES 0x101
#define NL_LEAST_SQUARES 0x102
#define NL_MAX_ITERATIONS 0x103
#define NL_THRESHOLD 0x104
#define NL_OMEGA 0x105
#define NL_SYMMETRIC 0x106
#define NL_USED_ITERATIONS 0x107
#define NL_ERROR 0x108
/* Solvers */
#define NL_CG 0x200
#define NL_JACOBI_CG 0x201
#define NL_SSOR_CG 0x202
#define NL_BICGSTAB 0x203
#define NL_GMRES 0x204
#define NL_SUPERLU_EXT 0x210
#define NL_PERMSUPERLU_EXT 0x211
/* Enable / Disable */
#define NL_NORMALIZE_ROWS 0x300
/* Row parameters */
#define RIGHT_HAND_SIDE 0x400
#define ROW_SCALING 0x401
/*
* Contexts
*/
NLAPI NLContext NLAPIENTRY nlNewContext() ;
NLAPI void NLAPIENTRY nlDeleteContext(NLContext context) ;
NLAPI void NLAPIENTRY nlMakeCurrent(NLContext context) ;
NLAPI NLContext NLAPIENTRY nlGetCurrent() ;
NLAPI NLboolean NLAPIENTRY nlInitExtension(char* extension) ;
/*
* State set/get
*/
NLAPI void NLAPIENTRY nlSolverParameterd(NLenum pname, NLdouble =
param) ;
NLAPI void NLAPIENTRY nlSolverParameteri(NLenum pname, NLint param) =
;
NLAPI void NLAPIENTRY nlRowParameterd(NLenum pname, NLdouble param) =
;
NLAPI void NLAPIENTRY nlRowParameteri(NLenum pname, NLint param) ;
NLAPI void NLAPIENTRY nlGetBooleanv(NLenum pname, NLboolean* params) =
;
NLAPI void NLAPIENTRY nlGetDoublev(NLenum pname, NLdouble* params) ;
NLAPI void NLAPIENTRY nlGetIntergerv(NLenum pname, NLint* params) ;
NLAPI void NLAPIENTRY nlEnable(NLenum pname) ;
NLAPI void NLAPIENTRY nlDisable(NLenum pname) ;
NLAPI NLboolean nlIsEnabled(NLenum pname) ;
/*
* Variables
*/
NLAPI void NLAPIENTRY nlSetVariable(NLuint index, NLdouble value) ;
NLAPI NLdouble NLAPIENTRY nlGetVariable(NLuint index) ;
NLAPI void NLAPIENTRY nlLockVariable(NLuint index) ;
NLAPI void NLAPIENTRY nlUnlockVariable(NLuint index) ;
NLAPI NLboolean NLAPIENTRY nlVariableIsLocked(NLuint index) ;
/*
* Begin/End
*/
NLAPI void NLAPIENTRY nlBegin(NLenum primitive) ;
NLAPI void NLAPIENTRY nlEnd(NLenum primitive) ;
NLAPI void NLAPIENTRY nlCoefficient(NLuint index, NLdouble value) ;
/*
* Solve
*/
NLAPI NLboolean NLAPIENTRY nlSolve() ;
=20
#ifdef __cplusplus
}
#endif
#endif
------=_NextPart_000_005B_01C436EE.A78F01A0
Content-Type: application/octet-stream;
name="nl.c"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="nl.c"
/*
* OpenNL: Numerical Library
* Copyright (C) 2004 Bruno Levy
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
* If you modify this software, you should include a notice giving the
* name of the person performing the modification, the date of =
modification,
* and the reason for such modification.
*
* Contact: Bruno Levy
*
* levy@loria.fr
*
* ISA Project
* LORIA, INRIA Lorraine,=20
* Campus Scientifique, BP 239
* 54506 VANDOEUVRE LES NANCY CEDEX=20
* FRANCE
*
* Note that the GNU General Public License does not permit =
incorporating
* the Software into proprietary programs.=20
*/
#include "nl.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
/************************************************************************=
************/
/* Assertions */
static void __nl_assertion_failed(char* cond, char* file, int line) {
fprintf(
stderr,=20
"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, char* file, int line
) {
fprintf(
stderr,=20
"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(char* file, int line) {
fprintf(
stderr,=20
"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__) ; \
} \
}=20
#define __nl_range_assert(x,min_val,max_val) { \
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__) ; \
}
#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)=20
#define __nl_debug_range_assert(x,min_val,max_val)=20
#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)=20
#define __nl_parano_range_assert(x,min_val,max_val)=20
#endif
/************************************************************************=
************/
/* Error-reporting function */
static void __nl_error(char* function, char* message) {
fprintf(stderr, "OpenNL error in %s(): %s\n", function, message) ;=20
}
/************************************************************************=
************/
/* classic macros */
#ifndef MIN
#define MIN(x,y) (((x) < (y)) ? (x) : (y))=20
#endif
#ifndef MAX
#define MAX(x,y) (((x) > (y)) ? (x) : (y))=20
#endif
/************************************************************************=
************/
/* memory management */
#define __NL_NEW(T) (T*)(calloc(1, sizeof(T)))=20
#define __NL_NEW_ARRAY(T,NB) (T*)(calloc((NB),sizeof(T)))=20
#define __NL_RENEW_ARRAY(T,x,NB) (T*)(realloc(x,(NB)*sizeof(T)))=20
#define __NL_DELETE(x) free(x); x =3D NULL=20
#define __NL_DELETE_ARRAY(x) free(x); x =3D NULL
#define __NL_CLEAR(T, x) memset(x, 0, sizeof(T))=20
#define __NL_CLEAR_ARRAY(T,x,NB) memset(x, 0, (NB)*sizeof(T))=20
/************************************************************************=
************/
/* 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 =3D 0 ;
c->capacity =3D 0 ;
c->coeff =3D NULL ;
}
static void __nlRowColumnDestroy(__NLRowColumn* c) {
__NL_DELETE_ARRAY(c->coeff) ;
#ifdef NL_PARANOID
__NL_CLEAR(__NLRowColumn, c) ;=20
#endif
}
static void __nlRowColumnGrow(__NLRowColumn* c) {
if(c->capacity !=3D 0) {
c->capacity =3D 2 * c->capacity ;
c->coeff =3D __NL_RENEW_ARRAY(__NLCoeff, c->coeff, c->capacity) =
;
} else {
c->capacity =3D 4 ;
c->coeff =3D __NL_NEW_ARRAY(__NLCoeff, c->capacity) ;
}
}
static void __nlRowColumnAdd(__NLRowColumn* c, NLint index, NLdouble =
value) {
NLuint i ;
for(i=3D0; i<c->size; i++) {
if(c->coeff[i].index =3D=3D index) {
c->coeff[i].value +=3D value ;
return ;
}
}
if(c->size =3D=3D c->capacity) {
__nlRowColumnGrow(c) ;
}
c->coeff[c->size].index =3D index ;
c->coeff[c->size].value =3D value ;
c->size++ ;
}
/* Does not check whether the index already exists */
static void __nlRowColumnAppend(__NLRowColumn* c, NLint index, NLdouble =
value) {
if(c->size =3D=3D c->capacity) {
__nlRowColumnGrow(c) ;
}
c->coeff[c->size].index =3D index ;
c->coeff[c->size].value =3D value ;
c->size++ ;
}
static void __nlRowColumnZero(__NLRowColumn* c) {
c->size =3D 0 ;
}
static void __nlRowColumnClear(__NLRowColumn* c) {
c->size =3D 0 ;
c->capacity =3D 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 =3D m ;
M->n =3D n ;
M->storage =3D storage ;
if(storage & __NL_ROWS) {
M->row =3D __NL_NEW_ARRAY(__NLRowColumn, m) ;
for(i=3D0; i<n; i++) {
__nlRowColumnConstruct(&(M->row[i])) ;
}
} else {
M->row =3D NULL ;
}
if(storage & __NL_COLUMNS) {
M->column =3D __NL_NEW_ARRAY(__NLRowColumn, n) ;
for(i=3D0; i<n; i++) {
__nlRowColumnConstruct(&(M->column[i])) ;
}
} else {
M->column =3D NULL ;
}
M->diag_size =3D MIN(m,n) ;
M->diag =3D __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=3D0; i<M->m; i++) {
__nlRowColumnDestroy(&(M->row[i])) ;
}
__NL_DELETE_ARRAY(M->row) ;
}
if(M->storage & __NL_COLUMNS) {
for(i=3D0; 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 =3D=3D j) {
M->diag[i] +=3D 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 __nlSparseMatrixZero( __NLSparseMatrix* M) {
NLuint i ;
if(M->storage & __NL_ROWS) {
for(i=3D0; i<M->m; i++) {
__nlRowColumnZero(&(M->row[i])) ;
}
}
if(M->storage & __NL_COLUMNS) {
for(i=3D0; i<M->n; i++) {
__nlRowColumnZero(&(M->column[i])) ;
}
}
__NL_CLEAR_ARRAY(NLdouble, M->diag, M->diag_size) ; =20
}
static void __nlSparseMatrixClear( __NLSparseMatrix* M) {
NLuint i ;
if(M->storage & __NL_ROWS) {
for(i=3D0; i<M->m; i++) {
__nlRowColumnClear(&(M->row[i])) ;
}
}
if(M->storage & __NL_COLUMNS) {
for(i=3D0; i<M->n; i++) {
__nlRowColumnClear(&(M->column[i])) ;
}
}
__NL_CLEAR_ARRAY(NLdouble, M->diag, M->diag_size) ; =20
}
/************************************************************************=
************/
/* SparseMatrix x Vector routines, internal helper routines */
static void __nlSparseMatrix_mult_rows_symmetric(
__NLSparseMatrix* A, NLdouble* x, NLdouble* y
) {
NLuint m =3D A->m ;
NLuint i,ij ;
__NLRowColumn* Ri =3D NULL ;
__NLCoeff* c =3D NULL ;
for(i=3D0; i<m; i++) {
y[i] =3D 0 ;
Ri =3D &(A->row[i]) ;
for(ij=3D0; ij<Ri->size; ij++) {
c =3D &(Ri->coeff[ij]) ;
y[i] +=3D c->value * x[c->index] ;
if(i !=3D c->index) {
y[c->index] +=3D c->value * x[i] ;
}
}
}
}
static void __nlSparseMatrix_mult_rows(
__NLSparseMatrix* A, NLdouble* x, NLdouble* y
) {
NLuint m =3D A->m ;
NLuint i,ij ;
__NLRowColumn* Ri =3D NULL ;
__NLCoeff* c =3D NULL ;
for(i=3D0; i<m; i++) {
y[i] =3D 0 ;
Ri =3D &(A->row[i]) ;
for(ij=3D0; ij<Ri->size; ij++) {
c =3D &(Ri->coeff[ij]) ;
y[i] +=3D c->value * x[c->index] ;
}
}
}
static void __nlSparseMatrix_mult_cols_symmetric(
__NLSparseMatrix* A, NLdouble* x, NLdouble* y
) {
NLuint n =3D A->n ;
NLuint j,ii ;
__NLRowColumn* Cj =3D NULL ;
__NLCoeff* c =3D NULL ;
for(j=3D0; j<n; j++) {
y[j] =3D 0 ;
Cj =3D &(A->column[j]) ;
for(ii=3D0; ii<Cj->size; ii++) {
c =3D &(Cj->coeff[ii]) ;
y[c->index] +=3D c->value * x[j] ;
if(j !=3D c->index) {
y[j] +=3D c->value * x[c->index] ;
}
}
}
}
static void __nlSparseMatrix_mult_cols(
__NLSparseMatrix* A, NLdouble* x, NLdouble* y
) {
__NL_CLEAR_ARRAY(NLdouble, y, A->m) ;
NLuint n =3D A->n ;
NLuint j,ii ;=20
__NLRowColumn* Cj =3D NULL ;
__NLCoeff* c =3D NULL ;
for(j=3D0; j<n; j++) {
Cj =3D &(A->column[j]) ;
for(ii=3D0; ii<Cj->size; ii++) {
c =3D &(Cj->coeff[ii]) ;
y[c->index] +=3D c->value * x[j] ;
}
}
}
/************************************************************************=
************/
/* SparseMatrix x Vector routines, main driver routine */
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) ;
}
}
}
/************************************************************************=
************/
/* NLContext data structure */
typedef struct {
NLdouble value ;
NLboolean locked ;
NLuint index ;
} __NLVariable ;
#define __NL_STATE_INITIAL 0
#define __NL_STATE_SYSTEM 1
#define __NL_STATE_MATRIX 2
#define __NL_STATE_ROW 3
#define __NL_STATE_MATRIX_CONSTRUCTED 4
#define __NL_STATE_SYSTEM_CONSTRUCTED 5
#define __NL_STATE_SOLVED 6
typedef struct {
NLenum state ;
__NLVariable* variable ;
NLuint n ;
__NLSparseMatrix M ;
__NLRowColumn af ;
__NLRowColumn al ;
__NLRowColumn xl ;
NLdouble* x ;
NLdouble* b ;
NLdouble right_hand_side ;
NLdouble row_scaling ;
NLenum solver ;
NLuint nb_variables ;
NLuint current_row ;
NLboolean least_squares ;
NLboolean symmetric ;
NLuint max_iterations ;
NLdouble threshold ;
NLdouble omega ;
NLboolean normalize_rows ;
NLboolean alloc_M ;
NLboolean alloc_af ;
NLboolean alloc_al ;
NLboolean alloc_xl ;
NLboolean alloc_variable ;
NLboolean alloc_x ;
NLboolean alloc_b ;
NLuint used_iterations ;
NLdouble error ;
} __NLContext ;
static __NLContext* __nlCurrentContext =3D NULL ;
NLContext nlNewContext() {
__NLContext* result =3D __NL_NEW(__NLContext) ;
result->state =3D __NL_STATE_INITIAL ;
result->solver =3D NL_BICGSTAB ;
result->max_iterations =3D 100 ;
result->threshold =3D 1e-6 ;
result->omega =3D 1.5 ;
result->row_scaling =3D 1.0 ;
return result ;
}
void nlDeleteContext(NLContext context_in) {
__NLContext* context =3D (__NLContext*)(context_in) ;
if(__nlCurrentContext =3D=3D context) {
__nlCurrentContext =3D NULL ;
}
if(context->alloc_M) {
__nlSparseMatrixDestroy(&context->M) ;
}
if(context->alloc_af) {
__nlRowColumnDestroy(&context->af) ;
}
if(context->alloc_al) {
__nlRowColumnDestroy(&context->al) ;
}
if(context->alloc_xl) {
__nlRowColumnDestroy(&context->xl) ;
}
if(context->alloc_variable) {
__NL_DELETE_ARRAY(context->variable) ;
}
if(context->alloc_x) {
__NL_DELETE_ARRAY(context->x) ;
}
if(context->alloc_b) {
__NL_DELETE_ARRAY(context->b) ;
}
#ifdef NL_PARANOID
__NL_CLEAR(__NLContext, context) ;
#endif
__NL_DELETE(context) ;
}
void nlMakeCurrent(NLContext context) {
__nlCurrentContext =3D (__NLContext*)(context) ;
}
NLContext nlGetCurrent() {
return __nlCurrentContext ;
}
NLboolean nlInitExtension(char* extension) {
return NL_FALSE ;
}
void __nlCheckState(NLenum state) {
__nl_assert(__nlCurrentContext->state =3D=3D state) ;
}
void __nlTransition(NLenum from_state, NLenum to_state) {
__nlCheckState(from_state) ;
__nlCurrentContext->state =3D to_state ;
}
/************************************************************************=
************/
/* Get/Set parameters */
void nlSolverParameterd(NLenum pname, NLdouble param) {
__nlCheckState(__NL_STATE_INITIAL) ;
switch(pname) {
case NL_SOLVER: {
__nlCurrentContext->solver =3D (NLenum)param ;
} break ;
case NL_NB_VARIABLES: {
__nl_assert(param > 0) ;
__nlCurrentContext->nb_variables =3D (NLuint)param ;
} break ;
case NL_LEAST_SQUARES: {
__nlCurrentContext->least_squares =3D (NLboolean)param ;
} break ;
case NL_MAX_ITERATIONS: {
__nl_assert(param > 0) ;
__nlCurrentContext->max_iterations =3D (NLuint)param ;
} break ;
case NL_THRESHOLD: {
__nl_assert(param >=3D 0) ;
__nlCurrentContext->threshold =3D (NLdouble)param ;
} break ;
case NL_OMEGA: {
__nl_range_assert(param,1.0,2.0) ;
__nlCurrentContext->omega =3D (NLdouble)param ;
} break ;
case NL_SYMMETRIC: {
__nlCurrentContext->symmetric =3D (NLboolean)param ; =20
}
default: {
__nl_assert_not_reached ;
} break ;
}
}
void nlSolverParameteri(NLenum pname, NLint param) {
__nlCheckState(__NL_STATE_INITIAL) ;
switch(pname) {
case NL_SOLVER: {
__nlCurrentContext->solver =3D (NLenum)param ;
} break ;
case NL_NB_VARIABLES: {
__nl_assert(param > 0) ;
__nlCurrentContext->nb_variables =3D (NLuint)param ;
} break ;
case NL_LEAST_SQUARES: {
__nlCurrentContext->least_squares =3D (NLboolean)param ;
} break ;
case NL_MAX_ITERATIONS: {
__nl_assert(param > 0) ;
__nlCurrentContext->max_iterations =3D (NLuint)param ;
} break ;
case NL_THRESHOLD: {
__nl_assert(param >=3D 0) ;
__nlCurrentContext->threshold =3D (NLdouble)param ;
} break ;
case NL_OMEGA: {
__nl_range_assert(param,1,2) ;
__nlCurrentContext->omega =3D (NLdouble)param ;
} break ;
case NL_SYMMETRIC: {
__nlCurrentContext->symmetric =3D (NLboolean)param ; =20
}
default: {
__nl_assert_not_reached ;
} break ;
}
}
void nlRowParameterd(NLenum pname, NLdouble param) {
__nlCheckState(__NL_STATE_MATRIX) ;
switch(pname) {
case RIGHT_HAND_SIDE: {
__nlCurrentContext->right_hand_side =3D param ;
} break ;
case ROW_SCALING: {
__nlCurrentContext->row_scaling =3D param ;
} break ;
}
}
void nlRowParameteri(NLenum pname, NLint param) {
__nlCheckState(__NL_STATE_MATRIX) ;
switch(pname) {
case RIGHT_HAND_SIDE: {
__nlCurrentContext->right_hand_side =3D (NLdouble)param ;
} break ;
case ROW_SCALING: {
__nlCurrentContext->row_scaling =3D (NLdouble)param ;
} break ;
}
}
void nlGetBooleanv(NLenum pname, NLboolean* params) {
switch(pname) {
case NL_LEAST_SQUARES: {
*params =3D __nlCurrentContext->least_squares ;
} break ;
case NL_SYMMETRIC: {
*params =3D __nlCurrentContext->symmetric ;
} break ;
default: {
__nl_assert_not_reached ;
} break ;
}
}
void nlGetDoublev(NLenum pname, NLdouble* params) {
switch(pname) {
case NL_SOLVER: {
*params =3D (NLdouble)(__nlCurrentContext->solver) ;
} break ;
case NL_NB_VARIABLES: {
*params =3D (NLdouble)(__nlCurrentContext->nb_variables) ;
} break ;
case NL_LEAST_SQUARES: {
*params =3D (NLdouble)(__nlCurrentContext->least_squares) ;
} break ;
case NL_MAX_ITERATIONS: {
*params =3D (NLdouble)(__nlCurrentContext->max_iterations) ;
} break ;
case NL_THRESHOLD: {
*params =3D (NLdouble)(__nlCurrentContext->threshold) ;
} break ;
case NL_OMEGA: {
*params =3D (NLdouble)(__nlCurrentContext->omega) ;
} break ;
case NL_SYMMETRIC: {
*params =3D (NLdouble)(__nlCurrentContext->symmetric) ;
} break ;
case NL_USED_ITERATIONS: {
*params =3D (NLdouble)(__nlCurrentContext->used_iterations) ;
} break ;
case NL_ERROR: {
*params =3D (NLdouble)(__nlCurrentContext->error) ;
} break ;
default: {
__nl_assert_not_reached ;
} break ;
}
}
void nlGetIntergerv(NLenum pname, NLint* params) {
switch(pname) {
case NL_SOLVER: {
*params =3D (NLint)(__nlCurrentContext->solver) ;
} break ;
case NL_NB_VARIABLES: {
*params =3D (NLint)(__nlCurrentContext->nb_variables) ;
} break ;
case NL_LEAST_SQUARES: {
*params =3D (NLint)(__nlCurrentContext->least_squares) ;
} break ;
case NL_MAX_ITERATIONS: {
*params =3D (NLint)(__nlCurrentContext->max_iterations) ;
} break ;
case NL_THRESHOLD: {
*params =3D (NLint)(__nlCurrentContext->threshold) ;
} break ;
case NL_OMEGA: {
*params =3D (NLint)(__nlCurrentContext->omega) ;
} break ;
case NL_SYMMETRIC: {
*params =3D (NLint)(__nlCurrentContext->symmetric) ;
} break ;
case NL_USED_ITERATIONS: {
*params =3D (NLint)(__nlCurrentContext->used_iterations) ;
} break ;
default: {
__nl_assert_not_reached ;
} break ;
}
}
/************************************************************************=
************/
/* Enable / Disable */
void nlEnable(NLenum pname) {
switch(pname) {
case NL_NORMALIZE_ROWS: {
__nl_assert(__nlCurrentContext->state !=3D __NL_STATE_ROW) ;
__nlCurrentContext->normalize_rows =3D NL_TRUE ;
} break ;
default: {
__nl_assert_not_reached ;
}
}
}
void nlDisable(NLenum pname) {
switch(pname) {
case NL_NORMALIZE_ROWS: {
__nl_assert(__nlCurrentContext->state !=3D __NL_STATE_ROW) ;
__nlCurrentContext->normalize_rows =3D NL_FALSE ;
} break ;
default: {
__nl_assert_not_reached ;
}
}
}
NLboolean nlIsEnabled(NLenum pname) {
switch(pname) {
case NL_NORMALIZE_ROWS: {
return __nlCurrentContext->normalize_rows ;
} break ;
default: {
__nl_assert_not_reached ;
}
}
return NL_FALSE ;
}
/************************************************************************=
************/
/* Get/Set Lock/Unlock variables */
void nlSetVariable(NLuint index, NLdouble value) {
__nlCheckState(__NL_STATE_SYSTEM) ;
__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables =
- 1) ;
__nlCurrentContext->variable[index].value =3D value ; =20
}
NLdouble nlGetVariable(NLuint index) {
__nl_assert(__nlCurrentContext->state !=3D __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 =3D NL_TRUE ;
}
void nlUnlockVariable(NLuint index) {
__nlCheckState(__NL_STATE_SYSTEM) ;
__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables =
- 1) ;
__nlCurrentContext->variable[index].locked =3D NL_FALSE ;
}
NLboolean nlVariableIsLocked(NLuint index) {
__nl_assert(__nlCurrentContext->state !=3D __NL_STATE_INITIAL) ;
__nl_parano_range_assert(index, 0, __nlCurrentContext->nb_variables =
- 1) ;
return __nlCurrentContext->variable[index].locked ;
}
/************************************************************************=
************/
/* System construction */
void __nlVariablesToVector() {
__nl_assert(__nlCurrentContext->alloc_x) ;
__nl_assert(__nlCurrentContext->alloc_variable) ;
NLuint i ;
for(i=3D0; i<__nlCurrentContext->nb_variables; i++) {
__NLVariable* v =3D &(__nlCurrentContext->variable[i]) ;
if(!v->locked) {
__nl_assert(v->index < __nlCurrentContext->n) ;
__nlCurrentContext->x[v->index] =3D v->value ;
}
}
}
void __nlVectorToVariables() {
__nl_assert(__nlCurrentContext->alloc_x) ;
__nl_assert(__nlCurrentContext->alloc_variable) ;
NLuint i ;
for(i=3D0; i<__nlCurrentContext->nb_variables; i++) {
__NLVariable* v =3D &(__nlCurrentContext->variable[i]) ;
if(!v->locked) {
__nl_assert(v->index < __nlCurrentContext->n) ;
v->value =3D __nlCurrentContext->x[v->index] ;
}
}
}
void __nlBeginSystem() {
__nlTransition(__NL_STATE_INITIAL, __NL_STATE_SYSTEM) ;
__nl_assert(__nlCurrentContext->nb_variables > 0) ;
__nlCurrentContext->variable =3D __NL_NEW_ARRAY(
__NLVariable, __nlCurrentContext->nb_variables
) ;
__nlCurrentContext->alloc_variable =3D NL_TRUE ;
}
void __nlEndSystem() {
__nlTransition(__NL_STATE_MATRIX_CONSTRUCTED, =
__NL_STATE_SYSTEM_CONSTRUCTED) ; =20
}
void __nlBeginMatrix() {
NLuint i ;
NLuint n =3D 0 ;
NLenum storage =3D __NL_ROWS ;
__nlTransition(__NL_STATE_SYSTEM, __NL_STATE_MATRIX) ;
for(i=3D0; i<__nlCurrentContext->nb_variables; i++) {
if(!__nlCurrentContext->variable[i].locked) {
__nlCurrentContext->variable[i].index =3D n ;
n++ ;
} else {
__nlCurrentContext->variable[i].index =3D ~0 ;
}
}
__nlCurrentContext->n =3D n ;
/* SSOR preconditioner requires rows and columns */
if(__nlCurrentContext->solver =3D=3D NL_SSOR_CG) {
storage =3D (storage | __NL_COLUMNS) ;
}
/* a least squares problem results in a symmetric matrix */
if(__nlCurrentContext->least_squares) {
__nlCurrentContext->symmetric =3D NL_TRUE ;
}
if(__nlCurrentContext->symmetric) {
storage =3D (storage | __NL_SYMMETRIC) ;
}
/* SuperLU storage does not support symmetric storage */
if(
__nlCurrentContext->solver =3D=3D NL_SUPERLU_EXT ||
__nlCurrentContext->solver =3D=3D NL_PERMSUPERLU_EXT=20
) {
storage =3D (storage & ~__NL_SYMMETRIC) ;
}
__nlSparseMatrixConstruct(&__nlCurrentContext->M, n, n, storage) ;
__nlCurrentContext->alloc_M =3D NL_TRUE ;
__nlCurrentContext->x =3D __NL_NEW_ARRAY(NLdouble, n) ;
__nlCurrentContext->alloc_x =3D NL_TRUE ;
=20
__nlCurrentContext->b =3D __NL_NEW_ARRAY(NLdouble, n) ;
__nlCurrentContext->alloc_b =3D NL_TRUE ;
__nlVariablesToVector() ;
__nlRowColumnConstruct(&__nlCurrentContext->af) ;
__nlCurrentContext->alloc_af =3D NL_TRUE ;
__nlRowColumnConstruct(&__nlCurrentContext->al) ;
__nlCurrentContext->alloc_al =3D NL_TRUE ;
__nlRowColumnConstruct(&__nlCurrentContext->xl) ;
__nlCurrentContext->alloc_xl =3D NL_TRUE ;
__nlCurrentContext->current_row =3D 0 ;
}
void __nlEndMatrix() {
__nlTransition(__NL_STATE_MATRIX, __NL_STATE_MATRIX_CONSTRUCTED) ; =
=20
=20
__nlRowColumnDestroy(&__nlCurrentContext->af) ;
__nlCurrentContext->alloc_af =3D NL_FALSE ;
__nlRowColumnDestroy(&__nlCurrentContext->al) ;
__nlCurrentContext->alloc_al =3D NL_FALSE ;
__nlRowColumnDestroy(&__nlCurrentContext->xl) ;
__nlCurrentContext->alloc_al =3D NL_FALSE ;
=20
if(!__nlCurrentContext->least_squares) {
__nl_assert(
__nlCurrentContext->current_row =3D=3D=20
__nlCurrentContext->n
) ;
}
}
void __nlBeginRow() {
__nlTransition(__NL_STATE_MATRIX, __NL_STATE_ROW) ;
__nlRowColumnZero(&__nlCurrentContext->af) ;
__nlRowColumnZero(&__nlCurrentContext->al) ;
}
void __nlScaleRow(NLdouble s) {
__NLRowColumn* af =3D &__nlCurrentContext->af ;
__NLRowColumn* al =3D &__nlCurrentContext->al ;
NLuint nf =3D af->size ;
NLuint nl =3D al->size ;
NLuint i ;
for(i=3D0; i<nf; i++) {
af->coeff[i].value *=3D s ;
}
for(i=3D0; i<nl; i++) {
al->coeff[i].value *=3D s ;
}
__nlCurrentContext->right_hand_side *=3D s ;
}
void __nlNormalizeRow(NLdouble weight) {
__NLRowColumn* af =3D &__nlCurrentContext->af ;
__NLRowColumn* al =3D &__nlCurrentContext->al ;
NLuint nf =3D af->size ;
NLuint nl =3D al->size ;
NLuint i ;
NLdouble norm =3D 0.0 ;
for(i=3D0; i<nf; i++) {
norm +=3D af->coeff[i].value * af->coeff[i].value ;
}
for(i=3D0; i<nl; i++) {
norm +=3D al->coeff[i].value * al->coeff[i].value ;
}
norm =3D sqrt(norm) ;
__nlScaleRow(weight / norm) ;
}
void __nlEndRow() {
__NLRowColumn* af =3D &__nlCurrentContext->af ;
__NLRowColumn* al =3D &__nlCurrentContext->al ;
__NLRowColumn* xl =3D &__nlCurrentContext->xl ;
__NLSparseMatrix* M =3D &__nlCurrentContext->M ;
NLdouble* b =3D __nlCurrentContext->b ;
NLuint nf =3D af->size ;
NLuint nl =3D al->size ;
NLuint current_row =3D __nlCurrentContext->current_row ;
NLuint i ;
NLuint j ;
NLdouble S ;
__nlTransition(__NL_STATE_ROW, __NL_STATE_MATRIX) ;
if(__nlCurrentContext->normalize_rows) {
__nlNormalizeRow(__nlCurrentContext->row_scaling) ;
} else {
__nlScaleRow(__nlCurrentContext->row_scaling) ;
}
if(__nlCurrentContext->least_squares) {
for(i=3D0; i<nf; i++) {
for(j=3D0; j<nf; j++) {
__nlSparseMatrixAdd(
M, af->coeff[i].index, af->coeff[j].index,
af->coeff[i].value * af->coeff[j].value
) ;
}
}
S =3D __nlCurrentContext->right_hand_side ;
for(j=3D0; j<nl; j++) {
S +=3D al->coeff[j].value * xl->coeff[j].value ;
}
for(i=3D0; i<nf; i++) {
b[ af->coeff[i].index ] -=3D af->coeff[i].value * S ;
}
} else {
for(i=3D0; i<nf; i++) {
__nlSparseMatrixAdd(M, current_row, af->coeff[i].index, =
af->coeff[i].value) ;
}
b[current_row] =3D __nlCurrentContext->right_hand_side ;
for(i=3D0; i<nl; i++) {
b[current_row] -=3D al->coeff[i].value * xl->coeff[i].value =
;
}
}
__nlCurrentContext->current_row++ ;
__nlCurrentContext->right_hand_side =3D 0.0 ; =20
__nlCurrentContext->row_scaling =3D 1.0 ;
}
void nlCoefficient(NLuint index, NLdouble value) {
__nlCheckState(__NL_STATE_ROW) ;
__nl_range_assert(index, 0, __nlCurrentContext->nb_variables - 1) ;
__NLVariable* v =3D &(__nlCurrentContext->variable[index]) ;
if(v->locked) {
__nlRowColumnAppend(&(__nlCurrentContext->al), 0, value) ;
__nlRowColumnAppend(&(__nlCurrentContext->xl), 0, v->value) ; =
=20
} else {
__nlRowColumnAppend(&(__nlCurrentContext->af), v->index, value) =
;
}
}
void nlBegin(NLenum prim) {
switch(prim) {
case NL_SYSTEM: {
__nlBeginSystem() ;
} break ;
case NL_MATRIX: {
__nlBeginMatrix() ;
} break ;
case NL_ROW: {
__nlBeginRow() ;
} break ;
default: {
__nl_assert_not_reached ;
}
}
}
void nlEnd(NLenum prim) {
switch(prim) {
case NL_SYSTEM: {
__nlEndSystem() ;
} break ;
case NL_MATRIX: {
__nlEndMatrix() ;
} break ;
case NL_ROW: {
__nlEndRow() ;
} break ;
default: {
__nl_assert_not_reached ;
}
}
}
/************************************************************************=
************/
/* BLAS routines */
/*
* daxpy
* ddot
* dscal
* dnrm2
* dcopy
*/
#ifndef NL_USE_BLAS
typedef NLint integer ;
typedef NLdouble doublereal ;
/* Subroutine */ int daxpy_(integer *n, doublereal *da, doublereal *dx,=20
integer *incx, doublereal *dy, integer *incy)
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i, m, ix, iy, mp1;
/* constant times a vector plus a vector. =20
uses unrolled loops for increments equal to one. =20
jack dongarra, linpack, 3/11/78. =20
modified 12/3/93, array(1) declarations changed to array(*) =20
=20
Parameter adjustments =20
Function Body */
#define DY(I) dy[(I)-1]
#define DX(I) dx[(I)-1]
if (*n <=3D 0) {
return 0;
}
if (*da =3D=3D 0.) {
return 0;
}
if (*incx =3D=3D 1 && *incy =3D=3D 1) {
goto L20;
}
/* code for unequal increments or equal increments =20
not equal to 1 */
ix =3D 1;
iy =3D 1;
if (*incx < 0) {
ix =3D (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy =3D (-(*n) + 1) * *incy + 1;
}
i__1 =3D *n;
for (i =3D 1; i <=3D *n; ++i) {
DY(iy) +=3D *da * DX(ix);
ix +=3D *incx;
iy +=3D *incy;
/* L10: */
}
return 0;
/* code for both increments equal to 1 =20
clean-up loop */
L20:
m =3D *n % 4;
if (m =3D=3D 0) {
goto L40;
}
i__1 =3D m;
for (i =3D 1; i <=3D m; ++i) {
DY(i) +=3D *da * DX(i);
/* L30: */
}
if (*n < 4) {
return 0;
}
L40:
mp1 =3D m + 1;
i__1 =3D *n;
for (i =3D mp1; i <=3D *n; i +=3D 4) {
DY(i) +=3D *da * DX(i);
DY(i + 1) +=3D *da * DX(i + 1);
DY(i + 2) +=3D *da * DX(i + 2);
DY(i + 3) +=3D *da * DX(i + 3);
/* L50: */
}
return 0;
} /* daxpy_ */
#undef DY
#undef DX
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal =
*dy,=20
integer *incy)
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
static integer i, m;
static doublereal dtemp;
static integer ix, iy, mp1;
/* forms the dot product of two vectors. =20
uses unrolled loops for increments equal to one. =20
jack dongarra, linpack, 3/11/78. =20
modified 12/3/93, array(1) declarations changed to array(*) =20
=20
Parameter adjustments =20
Function Body */
#define DY(I) dy[(I)-1]
#define DX(I) dx[(I)-1]
ret_val =3D 0.;
dtemp =3D 0.;
if (*n <=3D 0) {
return ret_val;
}
if (*incx =3D=3D 1 && *incy =3D=3D 1) {
goto L20;
}
/* code for unequal increments or equal increments =20
not equal to 1 */
ix =3D 1;
iy =3D 1;
if (*incx < 0) {
ix =3D (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy =3D (-(*n) + 1) * *incy + 1;
}
i__1 =3D *n;
for (i =3D 1; i <=3D *n; ++i) {
dtemp +=3D DX(ix) * DY(iy);
ix +=3D *incx;
iy +=3D *incy;
/* L10: */
}
ret_val =3D dtemp;
return ret_val;
/* code for both increments equal to 1 =20
clean-up loop */
L20:
m =3D *n % 5;
if (m =3D=3D 0) {
goto L40;
}
i__1 =3D m;
for (i =3D 1; i <=3D m; ++i) {
dtemp +=3D DX(i) * DY(i);
/* L30: */
}
if (*n < 5) {
goto L60;
}
L40:
mp1 =3D m + 1;
i__1 =3D *n;
for (i =3D mp1; i <=3D *n; i +=3D 5) {
dtemp =3D dtemp + DX(i) * DY(i) + DX(i + 1) * DY(i + 1) + DX(i + 2) *=20
DY(i + 2) + DX(i + 3) * DY(i + 3) + DX(i + 4) * DY(i + 4);
/* L50: */
}
L60:
ret_val =3D dtemp;
return ret_val;
} /* ddot_ */
#undef DY
#undef DX
/* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx,=20
integer *incx)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
static integer i, m, nincx, mp1;
/* scales a vector by a constant. =20
uses unrolled loops for increment equal to one. =20
jack dongarra, linpack, 3/11/78. =20
modified 3/93 to return if incx .le. 0. =20
modified 12/3/93, array(1) declarations changed to array(*) =20
=20
Parameter adjustments =20
Function Body */
#ifdef DX
#undef DX
#endif
#define DX(I) dx[(I)-1]
if (*n <=3D 0 || *incx <=3D 0) {
return 0;
}
if (*incx =3D=3D 1) {
goto L20;
}
/* code for increment not equal to 1 */
nincx =3D *n * *incx;
i__1 =3D nincx;
i__2 =3D *incx;
for (i =3D 1; *incx < 0 ? i >=3D nincx : i <=3D nincx; i +=3D *incx) =
{
DX(i) =3D *da * DX(i);
/* L10: */
}
return 0;
/* code for increment equal to 1 =20
clean-up loop */
L20:
m =3D *n % 5;
if (m =3D=3D 0) {
goto L40;
}
i__2 =3D m;
for (i =3D 1; i <=3D m; ++i) {
DX(i) =3D *da * DX(i);
/* L30: */
}
if (*n < 5) {
return 0;
}
L40:
mp1 =3D m + 1;
i__2 =3D *n;
for (i =3D mp1; i <=3D *n; i +=3D 5) {
DX(i) =3D *da * DX(i);
DX(i + 1) =3D *da * DX(i + 1);
DX(i + 2) =3D *da * DX(i + 2);
DX(i + 3) =3D *da * DX(i + 3);
DX(i + 4) =3D *da * DX(i + 4);
/* L50: */
}
return 0;
} /* dscal_ */
#undef DX
doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
{
/* System generated locals */
integer i__1, i__2;
doublereal ret_val, d__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static doublereal norm, scale, absxi;
static integer ix;
static doublereal ssq;
/* DNRM2 returns the euclidean norm of a vector via the function =20
name, so that =20
DNRM2 :=3D sqrt( x'*x ) =20
-- This version written on 25-October-1982. =20
Modified on 14-October-1993 to inline the call to DLASSQ. =20
Sven Hammarling, Nag Ltd. =20
=20
Parameter adjustments =20
Function Body */
#ifdef X
#undef X
#endif
#define X(I) x[(I)-1]
if (*n < 1 || *incx < 1) {
norm =3D 0.;
} else if (*n =3D=3D 1) {
norm =3D abs(X(1));
} else {
scale =3D 0.;
ssq =3D 1.;
/* The following loop is equivalent to this call to the LAPACK=20
=20
auxiliary routine: =20
CALL DLASSQ( N, X, INCX, SCALE, SSQ ) */
i__1 =3D (*n - 1) * *incx + 1;
i__2 =3D *incx;
for (ix =3D 1; *incx < 0 ? ix >=3D (*n-1)**incx+1 : ix <=3D =
(*n-1)**incx+1; ix +=3D *incx) {
if (X(ix) !=3D 0.) {
absxi =3D (d__1 =3D X(ix), abs(d__1));
if (scale < absxi) {
/* Computing 2nd power */
d__1 =3D scale / absxi;
ssq =3D ssq * (d__1 * d__1) + 1.;
scale =3D absxi;
} else {
/* Computing 2nd power */
d__1 =3D absxi / scale;
ssq +=3D d__1 * d__1;
}
}
/* L10: */
}
norm =3D scale * sqrt(ssq);
}
ret_val =3D norm;
return ret_val;
/* End of DNRM2. */
} /* dnrm2_ */
#undef X
/* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx,=20
doublereal *dy, integer *incy)
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i, m, ix, iy, mp1;
/* copies a vector, x, to a vector, y. =20
uses unrolled loops for increments equal to one. =20
jack dongarra, linpack, 3/11/78. =20
modified 12/3/93, array(1) declarations changed to array(*) =20
=20
Parameter adjustments =20
Function Body */
#define DY(I) dy[(I)-1]
#define DX(I) dx[(I)-1]
if (*n <=3D 0) {
return 0;
}
if (*incx =3D=3D 1 && *incy =3D=3D 1) {
goto L20;
}
/* code for unequal increments or equal increments =20
not equal to 1 */
ix =3D 1;
iy =3D 1;
if (*incx < 0) {
ix =3D (-(*n) + 1) * *incx + 1;
}
if (*incy < 0) {
iy =3D (-(*n) + 1) * *incy + 1;
}
i__1 =3D *n;
for (i =3D 1; i <=3D *n; ++i) {
DY(iy) =3D DX(ix);
ix +=3D *incx;
iy +=3D *incy;
/* L10: */
}
return 0;
/* code for both increments equal to 1 =20
clean-up loop */
L20:
m =3D *n % 7;
if (m =3D=3D 0) {
goto L40;
}
i__1 =3D m;
for (i =3D 1; i <=3D m; ++i) {
DY(i) =3D DX(i);
/* L30: */
}
if (*n < 7) {
return 0;
}
L40:
mp1 =3D m + 1;
i__1 =3D *n;
for (i =3D mp1; i <=3D *n; i +=3D 7) {
DY(i) =3D DX(i);
DY(i + 1) =3D DX(i + 1);
DY(i + 2) =3D DX(i + 2);
DY(i + 3) =3D DX(i + 3);
DY(i + 4) =3D DX(i + 4);
DY(i + 5) =3D DX(i + 5);
DY(i + 6) =3D DX(i + 6);
/* L50: */
}
return 0;
} /* dcopy_ */
#undef DX
#undef DY
#else
extern void daxpy_(=20
int *n, double *alpha, double *x,
int *incx, double *y, int *incy=20
) ;
extern double ddot_(=20
int *n, double *x, int *incx, double *y,
int *incy=20
) ;
extern double dnrm2_( int *n, double *x, int *incx ) ;
extern int dcopy_(int* n, double* dx, int* incx, double* dy, int* incy) =
;
extern void dscal_(int* n, double* alpha, double *x, int* incx) ;
#endif
/************************************************************************=
***********/
/* C wrappers for BLAS routines */
/* x <- a*x */
static void dscal( int n, double alpha, double *x, int incx ) {
dscal_(&n,&alpha,x,&incx);
}
/* y <- x */
static void dcopy(=20
int n, double *x, int incx, double *y, int incy=20
) {
dcopy_(&n,x,&incx,y,&incy);
}
/* y <- a*x+y */
static void daxpy(=20
int n, double alpha, double *x, int incx, double *y,
int incy=20
) {
daxpy_(&n,&alpha,x,&incx,y,&incy);
}
/* returns x^T*y */
static double ddot(=20
int n, double *x, int incx, double *y, int incy=20
) {
return ddot_(&n,x,&incx,y,&incy);
}
/** returns |x|_2 */
static double dnrm2( int n, double *x, int incx ) {
double d=3D0.;
while ( n-- ) {
d+=3D(*x)*(*x);=20
x+=3Dincx ;
}
return sqrt(d);
// return dnrm2_(&n,x,&incx);
}
/************************************************************************=
************/
/* Solvers */
static NLuint __nlSolve_CG(
__NLSparseMatrix* A,=20
NLdouble* b, NLdouble* x,=20
NLdouble eps, NLuint max_iter
) {
NLint N =3D A->n ;
NLdouble *g =3D __NL_NEW_ARRAY(NLdouble, N) ;
NLdouble *r =3D __NL_NEW_ARRAY(NLdouble, N) ;=20
NLdouble *p =3D __NL_NEW_ARRAY(NLdouble, N) ;
NLuint its=3D0;
NLdouble t, tau, sig, rho, gam;
NLdouble err=3Deps*eps*ddot(N,b,1,b,1);
__nlSparseMatrixMult(A,x,g);
daxpy(N,-1.,b,1,g,1);
dscal(N,-1.,g,1);
dcopy(N,g,1,r,1);
while ( ddot(N,g,1,g,1)>err && its < max_iter) {
__nlSparseMatrixMult(A,r,p);
rho=3Dddot(N,p,1,p,1);
sig=3Dddot(N,r,1,p,1);
tau=3Dddot(N,g,1,r,1);
t=3Dtau/sig;
daxpy(N,t,r,1,x,1);
daxpy(N,-t,p,1,g,1);
gam=3D(t*t*rho-tau)/tau;
dscal(N,gam,r,1);
daxpy(N,1.,g,1,r,1);
++its;
}
__NL_DELETE_ARRAY(g) ;
__NL_DELETE_ARRAY(r) ;
__NL_DELETE_ARRAY(p) ;
return its;
}=20
NLboolean nlSolve() {
__nlTransition(__NL_STATE_SYSTEM_CONSTRUCTED, __NL_STATE_SOLVED) ;
switch(__nlCurrentContext->solver) {
case NL_CG: {
__nlCurrentContext->used_iterations =3D=20
__nlSolve_CG(
&(__nlCurrentContext->M),
__nlCurrentContext->b,
__nlCurrentContext->x,
__nlCurrentContext->threshold,
__nlCurrentContext->max_iterations
) ;
} break ;
case NL_JACOBI_CG: {
__nl_assert_not_reached ;
} break ;
case NL_SSOR_CG: {
__nl_assert_not_reached ;
} break ;
case NL_BICGSTAB: {
__nl_assert_not_reached ;
} break ;
case NL_GMRES: {
__nl_assert_not_reached ;
} break ;
case NL_SUPERLU_EXT: {
__nl_assert_not_reached ;
} break ;
case NL_PERMSUPERLU_EXT: {
__nl_assert_not_reached ;
} break ;
default:
__nl_assert_not_reached ;
}
__nlVectorToVariables() ;
return NL_TRUE ;
}
------=_NextPart_000_005B_01C436EE.A78F01A0--