diff options
Diffstat (limited to 'r839.patch')
| -rw-r--r-- | r839.patch | 5467 |
1 files changed, 0 insertions, 5467 deletions
diff --git a/r839.patch b/r839.patch deleted file mode 100644 index 9fde401d3862..000000000000 --- a/r839.patch +++ /dev/null @@ -1,5467 +0,0 @@ -Index: setup.py -=================================================================== ---- setup.py (revision 838) -+++ setup.py (revision 839) -@@ -387,39 +387,7 @@ - if self.platform == "mac": - el.extra_compile_args += ["-Wno-error=return-type"] - el.extra_compile_args += self.flag_32bit + self.extra_c_flags -- -- if self.with_SUNDIALS: -- cordir_KINSOL_wSLU = os.path.join(self.assimulo_lib,'sundials_kinsol_core_wSLU.pyx') -- cordir_KINSOL = os.path.join(self.assimulo_lib,'sundials_kinsol_core.pyx') -- -- cordir_KINSOL_jmod_wSLU = os.path.join(self.assimulo_lib,'kinsol_jmod_wSLU.c') -- cordir_KINSOL_jmod = os.path.join(self.assimulo_lib,'kinsol_jmod.c') -- -- cordir_kinpinv = os.path.join(self.assimulo_lib,'kinpinv.c') -- cordir_kinslug = os.path.join(self.assimulo_lib,'kinslug.c') -- cordir_reg_routines = os.path.join(self.assimulo_lib,'reg_routines.c') -- if self.with_SLU: -- ext_list = ext_list + cythonize([cordir_KINSOL_wSLU], include_path=[".","assimulo",self.assimulo_lib]) -- ext_list[-1].sources += [cordir_KINSOL_jmod_wSLU,cordir_kinpinv,cordir_kinslug,cordir_reg_routines] -- ext_list[-1].include_dirs = [np.get_include(), self.SLUincdir, self.incdirs] -- ext_list[-1].library_dirs = [self.libdirs,self.SLUlibdir,self.BLASdir] -- ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial", "superlu_4.1",self.BLASname,'gfortran'] -- else: -- ext_list = ext_list + cythonize([cordir_KINSOL])#, include_path=[".","assimulo",self.assimulo_lib]) -- ext_list[-1].sources += [cordir_KINSOL_jmod,cordir_kinpinv] -- ext_list[-1].include_dirs = [np.get_include(), self.incdirs] -- ext_list[-1].library_dirs = [self.libdirs] -- ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial"] -- if self.SUNDIALS_version > (2,5,0): -- ext_list[-1].define_macros.append(("SUNDIALS_26", 1)) -- if self.debug_flag: -- ext_list[-1].extra_compile_args = ["-g", "-fno-strict-aliasing"] -- else: -- ext_list[-1].extra_compile_args = ["-O2", "-fno-strict-aliasing"] -- if self.platform == "mac": -- ext_list[-1].extra_compile_args += ["-Wno-error=return-type"] -- ext_list[-1].extra_compile_args += self.flag_32bit + self.extra_c_flags -- -+ - for el in ext_list: - if self.is_python3: - el.cython_directives = {"language_level": 3} -Index: src/kinsol.py -=================================================================== ---- assimulo/kinsol.py (revision 838) -+++ assimulo/kinsol.py (nonexistent) -@@ -1,350 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# 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 Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. --try: -- from .lib.sundials_kinsol_core_wSLU import KINSOL_wrap, KINError --except: -- from .lib.sundials_kinsol_core import KINSOL_wrap, KINError -- --from scipy.linalg import pinv2 --from scipy.optimize import fminbound --from numpy.linalg import solve --import numpy as N --import pylab as P --import operator as O --import re --from timeit import default_timer as timer --from assimulo.problem_algebraic import * -- --class KINSOL_Exception(Exception): -- pass -- -- --class KINSOL: -- -- def __init__(self,problem): -- """ -- Create the solver -- -- Parameters:: -- problem-- -- instance of ProblemAlgebraic found in problem_algebraic.py -- """ -- -- self.solver = KINSOL_wrap() -- -- # extract info from problem -- self.problem = problem -- if hasattr(self.problem,'_x0'): -- try: -- _x0 = self.problem.get_x0() -- except ProblemAlg_Exception: -- raise KINSOL_Exception("Problem has not implemented method 'get_x0'") -- else: -- raise KINSOL_Exception("Problem has no instance '_x0'") -- -- # calculate dimension -- try: -- if isinstance(_x0, int) or isinstance(_x0, float): -- self.x0 = [_x0] -- else: -- self.x0 = _x0 -- self.dim = len([N.array(self.x0, dtype=float)][0]) -- except ValueError: -- raise KINSOL_Exception("Initial guess must be a Numpy.array with either ints or floats.") -- -- # check for functions and test them -- try: -- tmp = self.problem.f(self.x0) -- self.norm_of_res = 10000*N.linalg.norm(self.x0) -- #self.norm_of_res = 1e20 -- self.func = self.problem.f -- except ProblemAlg_Exception: -- raise KINSOL_Exception("Problem has not implemented method 'f'") -- except IndexError: -- raise KINSOL_Exception("Problem has mismatching dimensions of f and initial guess") -- -- # Calculate 'stupid' scaling -- res = self.problem.f(self.x0) -- self.fscale = N.ones(self.dim) -- -- for i,val in zip(N.arange(self.dim),res): -- if val > 10: -- self.fscale[i] = 1.0/val -- -- # check for constraints and test them -- broken_constraints = [] -- if hasattr(self.problem, 'get_constraints'): -- self.constraints = self.problem.get_constraints() -- -- if self.constraints is not None: -- # test if constraints are of correct type -- if type(self.constraints).__name__ != 'ndarray': -- raise KINSOL_Exception("Constraints must be of type numpy.ndarray") -- -- if len(self.constraints) != len(self.x0): -- raise KINSOL_Exception("Constraints must have same length as x0") -- # Test if initial guess x0 is consistant with constraints -- for c,xi,i in zip(self.constraints,self.x0,N.arange(0,self.x0.__len__())): -- if re.search('float',type(c).__name__) is None: -- print("Type problem with: ", c, type(c).__name__) -- raise KINSOL_Exception("Constraints must contain floats.") -- if abs(c) > 2: -- raise KINSOL_Exception("Entries in constraint vector must be between -2 and 2, see documentation.") -- -- if c != 0.0: -- if c == 1.0: -- if xi < 0.0: -- broken_constraints.append(i) -- elif c == 2.0: -- if xi <= 0.0: -- broken_constraints.append(i) -- elif c == -1.0: -- if xi > 0.0: -- broken_constraints.append(i) -- elif c == -2.0: -- if xi >= 0.0: -- broken_constraints.append(i) -- else: -- raise KINSOL_Exception("Constraint vector contains illegal elements.") -- -- else: -- self.constraints = None -- -- if broken_constraints != []: -- print("Variables breaking initial constraint: ") -- for i in broken_constraints: -- self.problem.print_var_info(i) -- -- raise KINSOL_Exception("Initial guess does not fulfill applied constraints.") -- -- -- if hasattr(self.problem, 'check_constraints'): -- self.check_with_model = True -- else: -- self.check_with_model = False -- -- self._use_jac = True -- self.reg_count = 0 -- self.lin_count = 0 -- self.verbosity = 0 -- self.max_reg = 2.0 -- self.use_sparse = False -- self.reg_param = 0.0 -- self.exec_time = 0 -- self._use_ls = True -- self._use_fscale = False -- -- def set_jac_usage(self,use_jac): -- """ -- Set whether to use the jacobian supplied by the model -- or if we are to calculate it numericaly -- -- Parameters:: -- -- use_jac -- -- Boolean set to True if the jacobian is to be -- supplied by the model -- -- """ -- if type(use_jac).__name__ == 'bool': -- self._use_jac = use_jac -- else: -- raise KINSOL_Exception("The variable sent to 'set_jac_usage' must be a boolean.") -- -- def use_LineSearch(self,use_ls): -- """ -- Set whether the solver starts using the Linesearch Algorithm -- or not -- -- Parameters:: -- -- use_ls -- -- Boolean set to True if LineSearch is to -- be used. -- -- """ -- if type(use_ls).__name__ == 'bool': -- self._use_ls = use_ls -- else: -- raise KINSOL_Exception("The variable sent to 'use_LineSearch' must be a boolean.") -- -- def use_fscale(self,use_fscale): -- """ -- Set whether the solver starts using the Linesearch Algorithm -- or not -- -- Parameters:: -- -- use_ls -- -- Boolean set to True if LineSearch is to -- be used. -- -- """ -- if type(use_fscale).__name__ == 'bool': -- self._use_fscale = use_fscale -- else: -- raise KINSOL_Exception("The variable sent to 'use_fscale' must be a boolean.") -- -- -- def set_verbosity(self,verbosity): -- """ -- Method used to set the verbosity of KINSOL -- -- Parameters:: -- -- verbosity -- -- Integer set to one of 0,1,2,3 where 0 is no info and 3 the -- most info. For more information please see the documentation for KINSOL. -- """ -- type_name = type(verbosity).__name__ -- if re.search('int',type_name) is not None: -- -- # It is an integer, tes bounds -- if verbosity < 4 and verbosity > -1: -- self.verbosity = verbosity -- else: -- raise KINSOL_Exception("The variable sent to 'set_verbosity' must be either 0, 1, 2 or 3.") -- else: -- raise KINSOL_Exception("The variable sent to 'set_verbosity' must be an integer.") -- -- def set_sparsity(self,use_sparse): -- """ -- Method used to set if the problem should be treated as sparse by -- the linear solver in KINSOL. If the problem supplied has not implemented -- a method sparse_jac an exception will be thrown -- -- Parameters:: -- -- use_sparse -- -- Boolean set to True if the problem is to be treated as -- sparse and False otherwise. -- -- """ -- -- if hasattr(self.problem,'sparse_jac'): -- self.use_sparse = use_sparse -- else: -- raise KINSOL_Exception("The problem must have implemented a method 'sparse_jac' for sparsity to by used.") -- -- def set_reg_param(self,reg_param): -- """ -- Method used to set the regularization parameter. -- If set to zero the parameter will be set by the solver. -- -- Parameters:: -- -- reg_param -- -- Float larger or equal to zero. -- -- """ -- -- type_name = type(reg_param).__name__ -- if re.search('float',type_name) is not None: -- if reg_param < 0: -- raise KINSOL_Exception("Value sent to set_reg_param must be equal to or larger than zero.") -- else: -- self.reg_param = reg_param -- else: -- raise KINSOL_Exception("The variable sent to 'set_reg_param' must be a float.") -- -- def solve(self): -- """ -- Function called when solving function rhs_fct -- -- """ -- # check for jacobian and set it if present and to be used -- if self.use_sparse: -- if self._use_jac and hasattr(self.problem,'sparse_jac'): -- jac = self.problem.sparse_jac -- else: -- jac = None -- else: -- if self._use_jac and hasattr(self.problem,'jac'): -- jac = self.problem.jac -- else: -- jac = None -- -- # Initialize solver and solve -- -- solved = False -- local_min = False -- -- res = N.zeros(self.x0.__len__()) -- while (not solved) and self.reg_count < 2: -- try: -- if self._use_fscale: -- self.solver.KINSOL_init(self.func,self.x0,self.dim,jac,self.constraints,self.use_sparse,self.verbosity,self.norm_of_res,self.reg_param,self.fscale) -- else: -- self.solver.KINSOL_init(self.func,self.x0,self.dim,jac,self.constraints,self.use_sparse,self.verbosity,self.norm_of_res,self.reg_param,None) -- start = timer() -- res = self.solver.KINSOL_solve(not self._use_ls) -- stop = timer() -- self.exec_time += (stop - start) -- solved = True -- except KINError as error: -- if error.value == 42: -- # Try the heuristic -- if hasattr(self.problem, 'get_heuristic_x0'): -- print("----------------------------------------------------") -- print(" Solver stuck with zero step-length.") -- print("----------------------------------------------------") -- print("The following variables have start value zero") -- print("and min set to zero causing the zero step-lenght.") -- print("These settings are either set by default or by user.") -- print("") -- -- self.x0 = self.problem.get_heuristic_x0() -- self.reg_count += 1 -- -- print("") -- print("This setting (start and min to zero) can often") -- print("cause problem when initializing the system. ") -- print("") -- print("To avoid this the above variables have") -- print("their start attributes reset to one.") -- print("") -- print("Trying to solve the system again...") -- else: -- raise KINSOL_Exception("Regularization failed due to constraints, tried getting heuristic initial guess but failed.") -- -- -- elif (error.value == 2): -- print("---------------------------------------------------------") -- print("") -- print(" !!! WARNING !!!") -- print("") -- print(" KINSOL has returned a result but the algorithm has converged") -- print(" to a local minima, the initial values are NOT consistant!") -- print("") -- print("---------------------------------------------------------") -- solved = True -- local_min = True -- else: -- # Other error, send onward as exception -- self.problem.check_constraints(res) -- raise KINSOL_Exception(error.msg[error.value]) -- -- if not solved: -- self.solver.Free_KINSOL() -- raise KINSOL_Exception("Algorithm exited solution loop without finding a solution, please contact Assimulo support.") -- -- if self.check_with_model: -- self.problem.check_constraints(res) -- if not local_min: -- print("Problem sent to KINSOL solved.") -- -- return res -Index: src/lib/sundials_kinsol_core_wSLU.pxi -=================================================================== ---- assimulo/lib/sundials_kinsol_core_wSLU.pxi (revision 838) -+++ assimulo/lib/sundials_kinsol_core_wSLU.pxi (nonexistent) -@@ -1,293 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# 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 Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" -- --import numpy as np --from numpy cimport NPY_DOUBLE, npy_intp, NPY_INT -- --#============================================== --#External definitions from numpy headers --#============================================== --cdef extern from "numpy/arrayobject.h": -- -- ctypedef int npy_intp -- -- ctypedef extern class numpy.dtype [object PyArray_Descr]: -- cdef int type_num, elsize, alignment -- cdef char type, kind, byteorder, hasobject -- cdef object fields, typeobj -- -- ctypedef int intp -- ctypedef extern class numpy.ndarray [object PyArrayObject]: -- cdef char *data -- cdef int nd -- cdef intp *dimensions -- cdef intp *strides -- cdef int flags -- -- cdef object PyArray_SimpleNew(int nd, npy_intp* dims, int typenum) -- cdef object PyArray_SimpleNewFromData(int nd, npy_intp *dims, -- int typenum, void *data) -- void import_array() -- void* PyArray_GetPtr(ndarray aobj, npy_intp* ind) -- void *PyArray_DATA(ndarray aobj) -- --import_array() -- --#============================================== --#External definitions from SuperLU headers --#============================================== --cdef extern from "slu_ddefs.h": -- ctypedef int int_t -- -- ctypedef enum Stype_t: -- SLU_NC -- ctypedef enum Dtype_t: -- SLU_D -- ctypedef enum Mtype_t: -- SLU_GE -- ctypedef struct SuperMatrix: -- Stype_t Stype -- Dtype_t Dtype -- Mtype_t Mtype -- int_t nrow -- int_t ncol -- void *Store -- void dCreate_CompCol_Matrix(SuperMatrix *, int, int, int, double *,int *, int *, Stype_t, Dtype_t, Mtype_t) --#============================================== --# C headers --#============================================== --cdef extern from "string.h": -- void *memcpy(void *s1, void *s2, int n) --cdef extern from "stdio.h": -- int printf (char *__format, ...) --cdef extern from "stdlib.h": -- ctypedef unsigned long size_t -- void *calloc(size_t nelem, size_t elsize) -- --#============================================== --#External definitions from Sundials headers --#============================================== -- --# import data types --# basic data types --cdef extern from "sundials/sundials_types.h": -- ctypedef double realtype -- ctypedef bint booleantype -- --# sundials handling of vectors, so far only serial, --# parallel to be implemented later? --cdef extern from "sundials/sundials_nvector.h": -- cdef struct _generic_N_Vector: -- void* content -- ctypedef _generic_N_Vector* N_Vector -- --cdef extern from "nvector/nvector_serial.h": -- cdef struct _N_VectorContent_Serial: -- long int length -- realtype* data -- ctypedef _N_VectorContent_Serial* N_VectorContent_Serial -- N_Vector N_VNew_Serial(long int vec_length) -- void N_VDestroy_Serial(N_Vector v) -- void N_VPrint_Serial(N_Vector v) -- cdef N_Vector N_VMake_Serial(long int vec_length, realtype *v_data) -- N_Vector *N_VCloneVectorArray_Serial(int count, N_Vector w) -- N_Vector *N_VCloneVectorArrayEmpty_Serial(int count, N_Vector w) -- void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v) -- void N_VConst_Serial(realtype c, N_Vector z) -- --#Struct for handling the Jacobian data --cdef extern from "sundials/sundials_direct.h": -- cdef struct _DlsMat: -- int type -- int M -- int N -- int ldim -- int mu -- int ml -- int s_mu -- realtype *data -- int ldata -- realtype **cols -- ctypedef _DlsMat* DlsMat -- cdef realtype* DENSE_COL(DlsMat A, int j) -- --# KINSOL functions and routines --cdef extern from "kinsol/kinsol.h": -- # user defined functions -- ctypedef int (*KINSysFn)(N_Vector uu, N_Vector fval, void *user_data ) -- ctypedef void (*KINErrHandlerFn)(int error_code, char *module, char *function, char *msg, void *user_data) -- ctypedef void (*KINInfoHandlerFn)(char *module, char *function, char *msg, void *user_data) -- # initialization routines -- void *KINCreate() -- int KINInit(void *kinmem, KINSysFn func, N_Vector tmpl) -- -- # optional input spec. functions, -- # for specificationsdocumentation cf. kinsol.h line 218-449 -- int KINSetErrHandlerFn(void *kinmem, KINErrHandlerFn ehfun, void *eh_data) -- int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_data) -- int KINSetUserData(void *kinmem, void *user_data) -- int KINSetPrintLevel(void *kinmemm, int printfl) -- int KINSetNumMaxIters(void *kinmem, long int mxiter) -- int KINSetNoInitSetup(void *kinmem, booleantype noInitSetup) -- int KINSetNoResMon(void *kinmem, booleantype noNNIResMon) -- int KINSetMaxSetupCalls(void *kinmem, long int msbset) -- int KINSetMaxSubSetupCalls(void *kinmem, long int msbsetsub) -- int KINSetEtaForm(void *kinmem, int etachoice) -- int KINSetEtaConstValue(void *kinmem, realtype eta) -- int KINSetEtaParams(void *kinmem, realtype egamma, realtype ealpha) -- int KINSetResMonParams(void *kinmem, realtype omegamin, realtype omegamax) -- int KINSetResMonConstValue(void *kinmem, realtype omegaconst) -- int KINSetNoMinEps(void *kinmem, booleantype noMinEps) -- int KINSetMaxNewtonStep(void *kinmem, realtype mxnewtstep) -- int KINSetMaxBetaFails(void *kinmem, long int mxnbcf) -- int KINSetRelErrFunc(void *kinmem, realtype relfunc) -- int KINSetFuncNormTol(void *kinmem, realtype fnormtol) -- int KINSetScaledStepTol(void *kinmem, realtype scsteptol) -- int KINSetConstraints(void *kinmem, N_Vector constraints) -- int KINSetSysFunc(void *kinmem, KINSysFn func) -- -- # solver routine -- int KINSol(void *kinmem, N_Vector uu, int strategy, N_Vector u_scale, N_Vector f_scale) -- -- # optional output routines. -- # Documentation see kinsol.h line 670-735 -- int KINGetWorkSpace(void *kinmem, long int *lenrw, long int *leniw) -- int KINGetNumNonlinSolvIters(void *kinmem, long int *nniters) -- int KINGetNumFuncEvals(void *kinmem, long int *nfevals) -- int KINGetNumBetaCondFails(void *kinmem, long int *nbcfails) -- int KINGetNumBacktrackOps(void *kinmem, long int *nbacktr) -- int KINGetFuncNorm(void *kinmem, realtype *fnorm) -- int KINGetStepLength(void *kinmem, realtype *steplength) -- char *KINGetReturnFlagName(int flag) -- -- # fuction used to deallocate memory used by KINSOL -- void KINFree(void **kinmem) -- --# linear solver modules --cdef extern from "kinsol/kinsol_dense.h": -- int KINDense(void *kinmem, int N) -- --cdef extern from "kinpinv.h": -- int KINPinv(void *kinmem, int N) -- --cdef extern from "kinslug.h": -- int KINSLUG(void *kinmem, int N) -- --# functions used for supplying jacobian, and receiving info from linear solver --cdef extern from "kinsol/kinsol_direct.h": -- # user functions -- ctypedef int (*KINDlsDenseJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user functions to KINSOL -- int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac) -- -- # optional output fcts for linear direct solver -- int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINDlsGetLastFlag(void *kinmem, int *flag) -- char *KINDlsGetReturnFlagName(int flag) -- --cdef extern from "kinsol_jmod_wSLU.h": -- # user functions -- ctypedef int (*KINPinvJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- ctypedef int (*KINSLUGJacFn)(int N, N_Vector u, N_Vector fu, SuperMatrix *J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user jacobian to KINSOL -- int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) -- int KINSLUGSetJacFn(void *kinmem, KINSLUGJacFn jac) -- -- # functions used to set regularizationparameter -- int KINSLUGSetRegParam(void *kinmem, realtype reg_p) -- int KINPinvSetRegParam(void *kinmem, realtype reg_p) -- -- # optional output fcts for linear direct solver -- int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINPinvGetLastFlag(void *kinmem, int *flag) -- char *KINPinvGetReturnFlagName(int flag) -- -- --#========================= --# END SUNDIALS DEFINITIONS --#========================= -- --#=========================== --# Problem information struct --#=========================== --cdef class ProblemData: -- cdef: -- void *RHS # Residual or right-hand-side -- void *JAC -- int dim # Dimension of the problem -- --#================= --# Module functions --#================= -- --cdef inline N_Vector arr2nv(x): -- x=np.array(x) -- cdef long int n = len(x) -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- cdef N_Vector v=N_VNew_Serial(n) -- memcpy((<N_VectorContent_Serial>v.content).data, data_ptr, n*sizeof(double)) -- return v -- --cdef inline nv2arr(N_Vector v): -- cdef long int n = (<N_VectorContent_Serial>v.content).length -- cdef realtype* v_data = (<N_VectorContent_Serial>v.content).data -- cdef long int i -- import_array() -- cdef ndarray x=np.empty(n) -- #cdef ndarray x = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE,v_data) -- memcpy(x.data, v_data, n*sizeof(double)) -- return x -- --cdef inline realtype2arr(realtype *data, int n): -- """ Create new numpy array from realtype*""" -- import_array() -- cdef ndarray[realtype, ndim=1, mode='c'] x=np.empty(n) -- memcpy(x.data, data, n*sizeof(double)) -- return x -- --cdef inline arr2realtype(ndarray x, realtype* out_data, int n): -- """ Create new realtype* from numpy """ -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(double)) -- -- --cdef inline arr2int(ndarray x, int* out_data, int n): -- """ Create new int* from numpy array """ -- cdef ndarray[int, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(int)) -- -Index: src/lib/sundials_kinsol_includes.pxi -=================================================================== ---- assimulo/lib/sundials_kinsol_includes.pxi (revision 838) -+++ assimulo/lib/sundials_kinsol_includes.pxi (nonexistent) -@@ -1,267 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# 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 Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" -- --import numpy as np --from numpy cimport NPY_DOUBLE, npy_intp, NPY_INT -- --#============================================== --#External definitions from numpy headers --#============================================== --cdef extern from "numpy/arrayobject.h": -- -- ctypedef int npy_intp -- -- ctypedef extern class numpy.dtype [object PyArray_Descr]: -- cdef int type_num, elsize, alignment -- cdef char type, kind, byteorder, hasobject -- cdef object fields, typeobj -- -- ctypedef int intp -- ctypedef extern class numpy.ndarray [object PyArrayObject]: -- cdef char *data -- cdef int nd -- cdef intp *dimensions -- cdef intp *strides -- cdef int flags -- -- cdef object PyArray_SimpleNew(int nd, npy_intp* dims, int typenum) -- cdef object PyArray_SimpleNewFromData(int nd, npy_intp *dims, -- int typenum, void *data) -- void import_array() -- void* PyArray_GetPtr(ndarray aobj, npy_intp* ind) -- void *PyArray_DATA(ndarray aobj) -- --import_array() -- --#============================================== --# C headers --#============================================== --cdef extern from "string.h": -- void *memcpy(void *s1, void *s2, int n) --cdef extern from "stdio.h": -- int printf (char *__format, ...) --cdef extern from "stdlib.h": -- ctypedef unsigned long size_t -- void *calloc(size_t nelem, size_t elsize) -- --#============================================== --#External definitions from Sundials headers --#============================================== -- --# import data types --# basic data types --cdef extern from "sundials/sundials_types.h": -- ctypedef double realtype -- ctypedef bint booleantype -- --# sundials handling of vectors, so far only serial, --# parallel to be implemented later? --cdef extern from "sundials/sundials_nvector.h": -- cdef struct _generic_N_Vector: -- void* content -- ctypedef _generic_N_Vector* N_Vector -- --cdef extern from "nvector/nvector_serial.h": -- cdef struct _N_VectorContent_Serial: -- long int length -- realtype* data -- ctypedef _N_VectorContent_Serial* N_VectorContent_Serial -- N_Vector N_VNew_Serial(long int vec_length) -- void N_VDestroy_Serial(N_Vector v) -- void N_VPrint_Serial(N_Vector v) -- cdef N_Vector N_VMake_Serial(long int vec_length, realtype *v_data) -- N_Vector *N_VCloneVectorArray_Serial(int count, N_Vector w) -- N_Vector *N_VCloneVectorArrayEmpty_Serial(int count, N_Vector w) -- void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v) -- void N_VConst_Serial(realtype c, N_Vector z) -- --#Struct for handling the Jacobian data --cdef extern from "sundials/sundials_direct.h": -- cdef struct _DlsMat: -- int type -- int M -- int N -- int ldim -- int mu -- int ml -- int s_mu -- realtype *data -- int ldata -- realtype **cols -- ctypedef _DlsMat* DlsMat -- cdef realtype* DENSE_COL(DlsMat A, int j) -- --# KINSOL functions and routines --cdef extern from "kinsol/kinsol.h": -- # user defined functions -- ctypedef int (*KINSysFn)(N_Vector uu, N_Vector fval, void *user_data ) -- ctypedef void (*KINErrHandlerFn)(int error_code, char *module, char *function, char *msg, void *user_data) -- ctypedef void (*KINInfoHandlerFn)(char *module, char *function, char *msg, void *user_data) -- # initialization routines -- void *KINCreate() -- int KINInit(void *kinmem, KINSysFn func, N_Vector tmpl) -- -- # optional input spec. functions, -- # for specificationsdocumentation cf. kinsol.h line 218-449 -- int KINSetErrHandlerFn(void *kinmem, KINErrHandlerFn ehfun, void *eh_data) -- int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_data) -- int KINSetUserData(void *kinmem, void *user_data) -- int KINSetPrintLevel(void *kinmemm, int printfl) -- int KINSetNumMaxIters(void *kinmem, long int mxiter) -- int KINSetNoInitSetup(void *kinmem, booleantype noInitSetup) -- int KINSetNoResMon(void *kinmem, booleantype noNNIResMon) -- int KINSetMaxSetupCalls(void *kinmem, long int msbset) -- int KINSetMaxSubSetupCalls(void *kinmem, long int msbsetsub) -- int KINSetEtaForm(void *kinmem, int etachoice) -- int KINSetEtaConstValue(void *kinmem, realtype eta) -- int KINSetEtaParams(void *kinmem, realtype egamma, realtype ealpha) -- int KINSetResMonParams(void *kinmem, realtype omegamin, realtype omegamax) -- int KINSetResMonConstValue(void *kinmem, realtype omegaconst) -- int KINSetNoMinEps(void *kinmem, booleantype noMinEps) -- int KINSetMaxNewtonStep(void *kinmem, realtype mxnewtstep) -- int KINSetMaxBetaFails(void *kinmem, long int mxnbcf) -- int KINSetRelErrFunc(void *kinmem, realtype relfunc) -- int KINSetFuncNormTol(void *kinmem, realtype fnormtol) -- int KINSetScaledStepTol(void *kinmem, realtype scsteptol) -- int KINSetConstraints(void *kinmem, N_Vector constraints) -- int KINSetSysFunc(void *kinmem, KINSysFn func) -- -- # solver routine -- int KINSol(void *kinmem, N_Vector uu, int strategy, N_Vector u_scale, N_Vector f_scale) -- -- # optional output routines. -- # Documentation see kinsol.h line 670-735 -- int KINGetWorkSpace(void *kinmem, long int *lenrw, long int *leniw) -- int KINGetNumNonlinSolvIters(void *kinmem, long int *nniters) -- int KINGetNumFuncEvals(void *kinmem, long int *nfevals) -- int KINGetNumBetaCondFails(void *kinmem, long int *nbcfails) -- int KINGetNumBacktrackOps(void *kinmem, long int *nbacktr) -- int KINGetFuncNorm(void *kinmem, realtype *fnorm) -- int KINGetStepLength(void *kinmem, realtype *steplength) -- char *KINGetReturnFlagName(int flag) -- -- # fuction used to deallocate memory used by KINSOL -- void KINFree(void **kinmem) -- --# linear solver modules --cdef extern from "kinsol/kinsol_dense.h": -- int KINDense(void *kinmem, int N) -- --cdef extern from "kinpinv.h": -- int KINPinv(void *kinmem, int N) -- --# functions used for supplying jacobian, and receiving info from linear solver --cdef extern from "kinsol/kinsol_direct.h": -- # user functions -- ctypedef int (*KINDlsDenseJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user functions to KINSOL -- int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac) -- -- # optional output fcts for linear direct solver -- int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINDlsGetLastFlag(void *kinmem, int *flag) -- char *KINDlsGetReturnFlagName(int flag) -- --cdef extern from "kinsol_jmod.h": -- # user functions -- ctypedef int (*KINPinvJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user jacobian to KINSOL -- int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) -- -- # functions used to set regularizationparameter -- int KINPinvSetRegParam(void *kinmem, realtype reg_p) -- -- # optional output fcts for linear direct solver -- int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINPinvGetLastFlag(void *kinmem, int *flag) -- char *KINPinvGetReturnFlagName(int flag) -- -- --#========================= --# END SUNDIALS DEFINITIONS --#========================= -- --#=========================== --# Problem information struct --#=========================== --cdef class ProblemData: -- cdef: -- void *RHS # Residual or right-hand-side -- void *JAC -- int dim # Dimension of the problem -- --#================= --# Module functions --#================= -- --cdef inline N_Vector arr2nv(x): -- x=np.array(x) -- cdef long int n = len(x) -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- cdef N_Vector v=N_VNew_Serial(n) -- memcpy((<N_VectorContent_Serial>v.content).data, data_ptr, n*sizeof(double)) -- return v -- --cdef inline nv2arr(N_Vector v): -- cdef long int n = (<N_VectorContent_Serial>v.content).length -- cdef realtype* v_data = (<N_VectorContent_Serial>v.content).data -- cdef long int i -- import_array() -- cdef ndarray x=np.empty(n) -- #cdef ndarray x = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE,v_data) -- memcpy(x.data, v_data, n*sizeof(double)) -- return x -- --cdef inline realtype2arr(realtype *data, int n): -- """ Create new numpy array from realtype*""" -- import_array() -- cdef ndarray[realtype, ndim=1, mode='c'] x=np.empty(n) -- memcpy(x.data, data, n*sizeof(double)) -- return x -- --cdef inline arr2realtype(ndarray x, realtype* out_data, int n): -- """ Create new realtype* from numpy """ -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(double)) -- -- --cdef inline arr2int(ndarray x, int* out_data, int n): -- """ Create new int* from numpy array """ -- cdef ndarray[int, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(int)) -- -Index: src/lib/reg_routines.c -=================================================================== ---- assimulo/lib/reg_routines.c (revision 838) -+++ assimulo/lib/reg_routines.c (nonexistent) -@@ -1,202 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * 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 Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- * Johan Ylikiiskilä - johan.ylikiiskila@gmail.com -- * -- */ -- --#include "reg_routines.h" --#include <math.h> -- --SuperMatrix* regSparseMatrix(SuperMatrix *jac, double h) --{ -- SuperMatrix *res = NULL; -- int i=0,j=0,k=0,eli=0,elj=0,rowi=0,rowj=0,nnz=0,jnnz=0,jn=0,is_zero=0,max_nz=0; -- NCformat *Jstore = NULL; -- double *nzt = NULL, *jnzv=NULL,*nzval=NULL; -- double reg_param = h*h; -- int *rit = NULL, *colptr=NULL, *jri=NULL, *jcp=NULL,*rowind=NULL; -- Stype_t S = SLU_NC; -- Dtype_t D = SLU_D; -- Mtype_t M = SLU_GE; -- -- double *tmp_nzv = NULL; -- int *tmp_ri = NULL; -- -- -- /* Extract data */ -- Jstore = jac->Store; -- jn = jac->nrow; -- jnnz = Jstore->nnz; -- jnzv = Jstore->nzval; -- jri = Jstore->rowind; -- jcp = Jstore->colptr; -- -- /* Initialize data and allocate memory for result */ -- max_nz = 2*jnnz; -- res = calloc(1,sizeof(SuperMatrix)); -- nnz = 0; -- nzt = calloc(max_nz,sizeof(double)); -- rit = calloc(max_nz,sizeof(int)); -- colptr = calloc(jn+1,sizeof(int)); -- -- /* Calculate nnz,nzval,rowind and colptr */ -- colptr[0] = 0; -- -- /* Iterate over the columns of both J^T and J */ -- for (i = 0;i < jn;i++) { -- for (j = 0;j < jn;j++) { -- /* Get first index of first elt in both columns */ -- rowi = jcp[i]; -- rowj = jcp[j]; -- -- if (i==j) { -- /* Diagonal elt, perform regularization */ -- is_zero = 0; -- rit[nnz] = j; -- nzt[nnz++] += reg_param; -- } else { -- is_zero = 1; -- } -- -- /* Chack for matching and perform multiplication */ -- while ((rowi < jcp[i+1])&&(rowj < jcp[j+1])) { -- /* Get row numbers of the two current elts */ -- eli = jri[rowi]; -- elj = jri[rowj]; -- if (eli == elj) { -- /* Match! Non-zero elt to be inserted! */ -- /* If it is the first time, add to the non_zero structure */ -- if (is_zero) { -- is_zero = 0; -- rit[nnz++] = j; -- } -- nzt[nnz-1] += jnzv[rowi]*jnzv[rowj]; -- rowi++; -- rowj++; -- } else if (eli < elj) { -- rowi++; -- } else { -- rowj++; -- } -- } -- } -- /* New column, add info to colptr */ -- colptr[i+1] = nnz; -- -- /* -- * Check if we need to expand the allocated memory -- * and expand if necessary -- */ -- -- if (nnz > (max_nz-jn)) { -- tmp_nzv = calloc(2*max_nz,sizeof(double)); -- tmp_ri = calloc(2*max_nz,sizeof(int)); -- -- for (k = 0; k < nnz; k++) { -- tmp_nzv[k] = nzt[k]; -- tmp_ri[k] = rit[k]; -- } -- -- free(nzt); -- free(rit); -- -- nzt = tmp_nzv; -- rit = tmp_ri; -- -- tmp_nzv = NULL; -- tmp_ri = NULL; -- -- max_nz = 2*max_nz; -- } -- } -- -- /* Reallocate to the needed size of storage */ -- nzval = realloc(nzt,nnz*sizeof(double)); -- rowind = realloc(rit,nnz*sizeof(int)); -- -- /* Now create the matrix */ -- dCreate_CompCol_Matrix(res, jn, jn, nnz, nzval, rowind, colptr, S, D, M); -- return res ; --} -- --SuperMatrix* getRegRHS( SuperMatrix *jac, SuperMatrix *B) { -- SuperMatrix *res; -- NCformat *Jstore; -- DNformat *Bstore; -- double *data, *jac_data, *rhs_data; -- int *jac_ri, *jac_cp; -- int n,i,j,k ; -- -- /* Extract data from Jacobian */ -- Jstore = jac->Store; -- jac_data = Jstore->nzval; -- jac_ri = Jstore->rowind; -- jac_cp = Jstore->colptr; -- -- /* Extract data from right hand side */ -- n = B->nrow; -- Bstore = B->Store; -- rhs_data = Bstore->nzval; -- -- /* Allocate data for result */ -- data = calloc(n,sizeof(double)); -- -- /* Calculate data for result */ -- for (i = 0; i< n; i++) { -- for (j = jac_cp[i]; j < jac_cp[i+1]; j++) { -- k = jac_ri[j] ; -- data[i] += jac_data[j]*rhs_data[k]; -- } -- } -- -- /* Create the Matrix in SuperMatrix format */ -- res = calloc(1,sizeof(SuperMatrix)); -- dCreate_Dense_Matrix(res,n,1,data,n,SLU_DN,SLU_D,SLU_GE); -- -- return res; --} -- --double getRegParam(SuperMatrix *jac, SuperMatrix *B) { -- double res = 0; -- NCformat *Jstore; -- DNformat *Bstore; -- double data, *jac_data, *rhs_data; -- int *jac_ri, *jac_cp; -- int n,i,j,k ; -- -- /* Extract data from Jacobian */ -- Jstore = jac->Store; -- jac_data = Jstore->nzval; -- jac_ri = Jstore->rowind; -- jac_cp = Jstore->colptr; -- -- /* Extract data from right hand side */ -- n = B->nrow; -- Bstore = B->Store; -- rhs_data = Bstore->nzval; -- -- /* Calculate data for result */ -- for (i = 0; i< n; i++) { -- data = 0; -- for (j = jac_cp[i]; j < jac_cp[i+1]; j++) { -- k = jac_ri[j] ; -- data += jac_data[j]*rhs_data[k]; -- } -- res += data*data; -- } -- -- return sqrt(res); --} -Index: src/lib/sundials_kinsol_core.pyx -=================================================================== ---- assimulo/lib/sundials_kinsol_core.pyx (revision 838) -+++ assimulo/lib/sundials_kinsol_core.pyx (nonexistent) -@@ -1,371 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# 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 Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" -- --from __future__ import division -- --include "sundials_kinsol_includes.pxi" # Includes fcts from other header files --include "sundials_kinsol_core.pxi" # Includes the constants (textual include) -- --cdef int kin_res(N_Vector xv, N_Vector fval, void *problem_data): -- """ -- Residual fct called by KINSOL -- """ -- -- cdef: -- int i -- ProblemData pData = <ProblemData>problem_data -- ndarray[realtype, ndim = 1, mode = 'c'] rhs # Container holding return from user fct -- -- ndarray x = nv2arr(xv) -- realtype* resptr = (<N_VectorContent_Serial>fval.content).data -- -- try: -- rhs = (<object>pData.RHS)(x) -- for i in range(pData.dim): -- -- resptr[i] = rhs[i] -- -- return KIN_SUCCESS -- except: -- return KIN_REC_ERR -- --cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian, -- void *problem_data, N_Vector tmp1, N_Vector tmp2): -- """ -- This method is used to connect the Assimulo.Problem.jac to the Sundials -- Jacobian function. -- """ -- cdef: -- ProblemData pData = <ProblemData>problem_data -- realtype* col_i=DENSE_COL(Jacobian,0) -- ndarray x = nv2arr(xv) -- -- try: -- -- jac=(<object>pData.JAC)(x) -- #This needs further investigations: -- #memcpy(Jacobian.data,<realtype*>jac.data, pData.memSizeJac) -- -- for i in range(Neq): -- col_i = DENSE_COL(Jacobian, i) -- for j in range(Neq): -- col_i[j] = jac[j,i] -- -- return KINDLS_SUCCESS -- except: -- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description) -- -- --class KINError(Exception): -- """ -- Kinsol exception -- """ -- -- msg = {KIN_MEM_FAIL : 'Memory allocation failed.', -- KIN_MEM_NULL : 'KINSOL not properly created.', -- KIN_ILL_INPUT: 'There was an input error.', -- KIN_NO_MALLOC: 'Memory not allocated, call KINSOL_init(...)', -- KIN_LINESEARCH_NONCONV: 'Linesearch could not find a step big enough.', -- KIN_MAXITER_REACHED: 'Max number of iterations reached', -- KIN_MXNEWT_5X_EXCEEDED: 'Max size of newton step exceeded 5 times.', -- KIN_LINESEARCH_BCFAIL: 'Linesearch beta-test failed, probably because of too poor progress.', -- KIN_LINSOLV_NO_RECOVERY: 'psolve encountered an error, but preconditioner is current', -- KIN_LINIT_FAIL: 'Init of linear solver failed.', -- KIN_LSETUP_FAIL: 'pset (preconditioner setup fct) encountered an error', -- KIN_LSOLVE_FAIL: 'Error in either psolve or in linear solver routine.', -- KIN_SYSFUNC_FAIL: 'Call to RHS failed.', -- KIN_FIRST_SYSFUNC_ERR: 'Call to RHS failed on first call', -- KIN_REPTD_SYSFUNC_ERR: 'Call to RHS failed multiple times.', -- KIN_STEP_LT_STPTOL: 'Scaled step length too small. Either an approximate solution or a local minimum is reached. Check value of residual.', -- 234: 'KINSOL not compiled with SuperLU, sparse functionality not available.'} -- value = 0 -- def __init__(self, value): -- self.value = value -- -- def __str__(self): -- try: -- return repr(self.msg[self.value]) -- except KeyError: -- return repr('KINSOL failed with flag %s.'%(self.value)) -- --cdef class KINSOL_wrap: -- -- cdef: -- void* solver -- ProblemData pData # A struct containing problem data -- #ProblemData *ppData # Pointer to above mentioned struct -- N_Vector x_cur,x_init, x_scale, f_scale, con_nv -- booleantype noInitSetup, sparse -- realtype norm_of_res,reg_param -- ndarray con -- int print_level -- ndarray fscale -- -- def __cinit__(self): -- self.pData = ProblemData() #Create a new problem struct -- #self.ppData = &self.pData -- -- cpdef KINSOL_set_problem_info(self, RHS, dim, JAC= None): -- """ -- Sets the problem info to the desired struct -- """ -- -- # Sets residual function -- self.pData.RHS = <void*>RHS -- self.pData.dim = dim -- if JAC != None: -- self.pData.JAC = <void*>JAC -- else: -- self.pData.JAC = NULL -- -- def KINSOL_init(self,RHS,x0,dim, JAC = None, con = None,sparse = False, print_level = 0,norm_of_res = 0,reg_param=0,fscale=None): -- """ -- Initializes solver -- """ -- cdef int flag -- # Set problem info -- self.noInitSetup = False -- self.print_level = print_level -- self.norm_of_res = norm_of_res -- self.reg_param = reg_param -- self.KINSOL_set_problem_info(RHS,dim,JAC) -- self.con = con -- self.sparse = sparse -- self.fscale = fscale -- -- -- # Create initial guess from the supplied numpy array -- # print "x0 got in KINSOL wrapper: ", x0 -- self.x_cur = arr2nv(x0) -- self.x_init = arr2nv(x0) -- -- -- if self.solver == NULL: # solver runs for the first time -- -- # Create solver -- self.solver = KINCreate() -- if self.solver == NULL: -- raise KINError(KIN_MEM_FAIL) -- if self.print_level >0: -- print "KINSOL created" -- -- # Stop solver from performing precond setup -- flag = KINSetNoInitSetup(self.solver, self.noInitSetup) -- if flag < 0: -- raise KINError(flag) -- -- flag = KINSetPrintLevel(self.solver, self.print_level) -- if flag < 0: -- raise KINError(flag) -- -- if self.norm_of_res >0: -- flag = KINSetMaxNewtonStep(self.solver, self.norm_of_res) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetMaxNewtonStep(self.solver, 1e20) -- if flag < 0: -- raise KINError(flag) -- -- # If the user has specified constraints, connect them -- if self.con is not None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(self.con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- -- # Allocate internal memory -- flag = KINInit(self.solver,kin_res,self.x_cur) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "KINSOL initialized" -- -- # Link to linear solver -- if self.sparse: -- raise KINError(234) -- else: -- flag = KINPinv(self.solver,self.pData.dim) -- #flag = KINDense(self.solver,self.pData.dim) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Linear solver connected" -- -- # Set regularization parameter, if necessary -- if self.reg_param != 0.0: -- if self.sparse: -- raise KINError(234) -- else: -- flag = KINPinvSetRegParam(self.solver, self.reg_param) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Reg param ", self.reg_param, " set" -- else: -- # If the user has specified constraints, connect them -- if con is not None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetConstraints(self.solver, NULL) -- if flag < 0: -- raise KINError(flag) -- -- # If the user supplied a Jacobien, link it to the solver -- if self.pData.JAC != NULL: -- if sparse: -- raise KINError(234) -- else: -- flag = KINPinvSetJacFn(self.solver,kin_jac) -- #flag = KINDlsSetDenseJacFn(self.solver,kin_jac) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Jacobian supplied by user connected" -- else: -- if sparse: -- raise KINError(234) -- else: -- flag = KINPinvSetJacFn(self.solver,NULL) -- #flag = KINDlsSetDenseJacFn(self.solver,NULL) -- -- if flag < 0: -- raise KINError(flag) -- -- # Link the solver to the supplied problem data -- flag = KINSetUserData(self.solver,<void*>self.pData) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "User data set" -- -- # Create scaling vectors filled with ones -- # since the problem is assumed to be scaled -- self.x_scale = arr2nv(np.ones(self.pData.dim)) -- if self.fscale is not None: -- self.f_scale = arr2nv(self.fscale) -- else: -- self.f_scale = arr2nv(np.ones(self.pData.dim)) -- -- -- def KINSOL_solve(self,wo_LineSearch): -- """ -- Function that should be called after a call to KINSOL_init -- solves the function supplied as RHS -- """ -- -- cdef: -- int flag, DLSflag, lsflag, count,ind -- long int nb_iters -- ndarray x0 -- if self.print_level >0: -- print "Calling solver..." -- -- if not wo_LineSearch: -- # Use LineSearch -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_LINESEARCH,self.x_scale,self.f_scale) -- count = 0 -- if flag > 1: -- # Probably stuck at a minimum -- self.x_cur = self.x_init -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- -- if flag > 1: -- # Stuck at minimum again -- raise KINError(flag) -- else: -- # Don't use linesearch -- flag = -7 -- -- while flag <0: -- # Get nbr of non-linear iterations -- lsflag = KINGetNumNonlinSolvIters(<void*>self.solver, &nb_iters) -- if lsflag != 0: -- raise KINError(lsflag) -- -- -- if (flag == -5 and nb_iters == 1): -- # Cannot take a step bigger than zero on the first iteration -- # Try with heuristic initial guess -- raise KINError(42) -- -- if (flag == -7 or flag == -6 or (flag == -5 and nb_iters > 1)): -- if count == 0: -- # First time -- print "Calling solver again as simple newton." -- self.x_cur = self.x_init -- DLSflag = KINSetMaxSetupCalls(self.solver, 1); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 1: -- # second time -- print "Turning of residual monitoring" -- DLSflag = KINSetNoResMon(self.solver, True); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 3: -- # fourth time, release the bound on the step -- -- print "Remove limit on stepsize" -- DLSflag = KINSetMaxNewtonStep(self.solver, 1.0e34) -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- count += 1 -- if count >= 10: -- if self.print_level >0: -- print "Tried with simple newton ten times" -- raise KINError(-7) -- -- else: -- -- print "Error in solve, flag: ", flag -- lsflag = KINPinvGetLastFlag(self.solver, &lsflag) -- #lsflag = KINDlsGetLastFlag(self.solver, &lsflag) -- if lsflag != 0: -- if lsflag <0: -- print "Last flag from Linear solver: ", lsflag -- else: -- print "Jacobian singular at column ", lsflag -- raise KINError(flag) -- -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- if self.print_level >0: -- print "Problem solved, returning result" -- -- return nv2arr(self.x_cur) -Index: src/lib/kinsol_jmod.c -=================================================================== ---- assimulo/lib/kinsol_jmod.c (revision 838) -+++ assimulo/lib/kinsol_jmod.c (nonexistent) -@@ -1,481 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --/* -- * ================================================================= -- * IMPORTED HEADER FILES -- * ================================================================= -- */ -- --#include <stdio.h> --#include <stdlib.h> -- --#include "kinsol/kinsol_impl.h" --#include "kinsol_jmod_impl.h" --#include <sundials/sundials_math.h> -- --/* -- * ================================================================= -- * FUNCTION SPECIFIC CONSTANTS -- * ================================================================= -- */ -- --/* Constant for DQ Jacobian approximation */ --#define MIN_INC_MULT RCONST(1000.0) -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define uround (kin_mem->kin_uround) --#define func (kin_mem->kin_func) --#define user_data (kin_mem->kin_user_data) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#define sfdotJp (kin_mem->kin_sfdotJp) --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define bjac (kinpinv_mem->d_bjac) --#define J (kinpinv_mem->d_J) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define RTR (kinpinv_mem->d_RTR) --#define reg_param (kinpinv_mem->d_reg_param) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPinvSetJacFn -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (jac != NULL) { -- jacDQ = FALSE; -- djac = jac; -- } else { -- jacDQ = TRUE; -- } -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvSetRegParam -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetRegParam(void *kinmem, realtype reg_p) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- reg_param = reg_p; -- -- return(KINPINV_SUCCESS); --} -- -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetWorkSpace -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetWorkSpace(void *kinmem, long int *lenrwLS, long int *leniwLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (mtype == SUNDIALS_DENSE) { -- *lenrwLS = n*n; -- *leniwLS = n; -- } else if (mtype == SUNDIALS_BAND) { -- *lenrwLS = n*(smu + mu + 2*ml + 2); -- *leniwLS = n; -- } -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumJacEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumJacEvals(void *kinmem, long int *njevals) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *njevals = nje; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumFuncEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumFuncEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumGuncEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *nfevalsLS = nfeDQ; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetLastFlag -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetLastFlag(void *kinmem, int *flag) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *flag = last_flag; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetReturnFlagName -- * ----------------------------------------------------------------- -- */ -- --char *KINPinvGetReturnFlagName(int flag) --{ -- char *name; -- -- name = (char *)malloc(30*sizeof(char)); -- -- switch(flag) { -- case KINPINV_SUCCESS: -- sprintf(name, "KINPINV_SUCCESS"); -- break; -- case KINPINV_MEM_NULL: -- sprintf(name, "KINPINV_MEM_NULL"); -- break; -- case KINPINV_LMEM_NULL: -- sprintf(name, "KINPINV_LMEM_NULL"); -- break; -- case KINPINV_ILL_INPUT: -- sprintf(name, "KINPINV_ILL_INPUT"); -- break; -- case KINPINV_MEM_FAIL: -- sprintf(name, "KINPINV_MEM_FAIL"); -- break; -- default: -- sprintf(name, "NONE"); -- } -- -- return(name); --} -- -- --/* -- * ================================================================= -- * DQ JACOBIAN APPROXIMATIONS -- * ================================================================= -- */ -- -- -- -- --/* -- * ----------------------------------------------------------------- -- * kinPinvDQJac -- * ----------------------------------------------------------------- -- * This routine generates a dense difference quotient approximation to -- * the Jacobian of F(u). It assumes that a dense matrix of type -- * DlsMat is stored column-wise, and that elements within each column -- * are contiguous. The address of the jth column of J is obtained via -- * the macro DENSE_COL and this pointer is associated with an N_Vector -- * using the N_VGetArrayPointer/N_VSetArrayPointer functions. -- * Finally, the actual computation of the jth column of the Jacobian is -- * done with a call to N_VLinearSum. -- * -- * The increment used in the finite-difference approximation -- * J_ij = ( F_i(u+sigma_j * e_j) - F_i(u) ) / sigma_j -- * is -- * sigma_j = max{|u_j|, |1/uscale_j|} * sqrt(uround) -- * -- * Note: uscale_j = 1/typ(u_j) -- * -- * NOTE: Any type of failure of the system function her leads to an -- * unrecoverable failure of the Jacobian function and thus -- * of the linear solver setup function, stopping KINSOL. -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2) --{ -- realtype inc, inc_inv, ujsaved, ujscale, sign; -- realtype *tmp2_data, *u_data, *uscale_data; -- N_Vector ftemp, jthCol; -- long int j; -- int retval = 0; -- -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* data points to kin_mem */ -- kin_mem = (KINMem) data; -- kinpinv_mem = (KINPinvMem) lmem; -- -- /* Save pointer to the array in tmp2 */ -- tmp2_data = N_VGetArrayPointer(tmp2); -- -- /* Rename work vectors for readibility */ -- ftemp = tmp1; -- jthCol = tmp2; -- -- /* Obtain pointers to the data for u and uscale */ -- u_data = N_VGetArrayPointer(u); -- uscale_data = N_VGetArrayPointer(uscale); -- -- /* This is the only for loop for 0..N-1 in KINSOL */ -- -- for (j = 0; j < N; j++) { -- -- /* Generate the jth col of Jac(u) */ -- -- N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol); -- -- ujsaved = u_data[j]; -- ujscale = ONE/uscale_data[j]; -- sign = (ujsaved >= ZERO) ? ONE : -ONE; --#if SUNDIALS_26 -- inc = sqrt_relfunc*SUNMAX(SUNRabs(ujsaved), ujscale)*sign; --#else -- inc = sqrt_relfunc*MAX(ABS(ujsaved), ujscale)*sign; --#endif -- u_data[j] += inc; -- -- retval = func(u, ftemp, user_data); -- nfeDQ++; -- if(retval > 0) { -- /* try to recover by stepping in the opposite direction */ -- inc = -inc; -- u_data[j] = ujsaved + inc; -- -- retval = func(u, ftemp, user_data); -- nfeDQ++; -- } -- if (retval != 0) break; -- -- u_data[j] = ujsaved; -- -- inc_inv = ONE/inc; -- N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol); -- -- } -- -- /* Restore original array pointer in tmp2 */ -- N_VSetArrayPointer(tmp2_data, tmp2); -- -- return(retval); --} -Index: src/lib/kinsol_jmod_impl_wSLU.h -=================================================================== ---- assimulo/lib/kinsol_jmod_impl_wSLU.h (revision 838) -+++ assimulo/lib/kinsol_jmod_impl_wSLU.h (nonexistent) -@@ -1,168 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifndef _KINJMOD_IMPL_H --#define _KINJMOD_IMPL_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include "kinsol_jmod_wSLU.h" --#include "slu_ddefs.h" -- -- --/* -- * ----------------------------------------------------------------- -- * Types: KINPinvMemRec, KINPinvMem -- * ----------------------------------------------------------------- -- * The type KINPinvMem is pointer to a KINPinvMemRec. -- * This structure contains KINPinv/KINSLUG solver-specific data. -- * ----------------------------------------------------------------- -- */ -- --typedef struct KINPinvMemRec { -- -- int d_type; /* SUNDIALS_DENSE or SUNDIALS_BAND */ -- -- int d_n; /* problem dimension */ -- -- int d_ml; /* lower bandwidth of Jacobian */ -- int d_mu; /* upper bandwidth of Jacobian */ -- int d_smu; /* upper bandwith of M = MIN(N-1,d_mu+d_ml) */ -- -- booleantype d_jacDQ; /* TRUE if using internal DQ Jacobian approx. */ -- KINPinvJacFn d_djac; /* dense Jacobian routine to be called */ -- KINSLUGJacFn d_spjac; /* sparse Jacobian routine to be called */ -- -- void *d_J_data; /* J_data is passed to djac or bjac */ -- -- DlsMat d_J; /* problem Jacobian */ -- -- SuperMatrix *sp_J; /* sparse problem jacobian */ -- SuperMatrix *sp_L; /* L in the sparse LU factorization */ -- SuperMatrix *sp_U; /* U in the sparse LU factorization */ -- SuperMatrix *sp_B; /* sparse right hand side */ -- int *sp_perm_r; /* row permutations from partial pivoting */ -- int *sp_perm_c; /* column permutation vector */ -- -- SuperMatrix *sp_JTJ; /* Matrix needed for regularisation */ -- -- superlu_options_t *sp_options; /* options struct for SuperLU */ -- SuperLUStat_t *sp_stat; /* statistcis struct for SuperLU */ -- -- int *d_pivots; /* pivot array for PM = LU */ -- realtype *d_beta; -- realtype d_reg_param; /* Regularization parameter */ -- -- long int d_nje; /* no. of calls to jac */ -- -- long int d_nfeDQ; /* no. of calls to F due to DQ Jacobian approx. */ -- -- int d_last_flag; /* last error return flag */ -- -- DlsMat d_JTJ; -- booleantype d_regularized; /* Boolean set to true if problem is regularized*/ -- booleantype d_redojac; -- --} *KINPinvMem; -- -- --/* -- * ----------------------------------------------------------------- -- * Prototypes of internal functions -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2); -- --/* -- * ----------------------------------------------------------------- -- * Error Messages -- * ----------------------------------------------------------------- -- */ -- --#define MSGD_KINMEM_NULL "KINSOL memory is NULL." --#define MSGD_BAD_NVECTOR "A required vector operation is not implemented." --#define MSGD_MEM_FAIL "A memory request failed." --#define MSGD_LMEM_NULL "Linear solver memory is NULL." --#define MSGD_BAD_SIZES "Illegal bandwidth parameter(s). Must have 0 <= ml, mu <= N-1." --#define MSGD_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner." -- --#ifdef __cplusplus --} --#endif -- --#endif -Index: src/lib/reg_routines.h -=================================================================== ---- assimulo/lib/reg_routines.h (revision 838) -+++ assimulo/lib/reg_routines.h (nonexistent) -@@ -1,43 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * 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 Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- * Johan Ylikiiskilä - johan.ylikiiskila@gmail.com -- * -- */ -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#ifndef _REG_ROUTINES_H --#define _REG_ROUTINES_H -- -- --#include "slu_ddefs.h" -- -- /* -- * Routines used in regularisation -- */ -- --SuperMatrix* regSparseMatrix( SuperMatrix *jac, double h); -- --SuperMatrix* getRegRHS( SuperMatrix *jac, SuperMatrix *B); -- --double getRegParam(SuperMatrix *jac, SuperMatrix *B); --#endif -- --#ifdef __cplusplus --} --#endif -Index: src/lib/kinslug.c -=================================================================== ---- assimulo/lib/kinslug.c (revision 838) -+++ assimulo/lib/kinslug.c (nonexistent) -@@ -1,678 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * 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 Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- */ -- --/* -- * This is a linear solver using Tikhonov regularization and the sparse -- * solver SuperLU -- */ --#include <stdio.h> --#include <stdlib.h> --#include <math.h> -- --#include "slu_ddefs.h" --#include "supermatrix.h" --#include "reg_routines.h" --#include "kinsol_jmod_impl_wSLU.h" --#include "kinslug.h" --#include "kinsol/kinsol_impl.h" -- --#include <nvector/nvector_serial.h> --#include <sundials/sundials_direct.h> --#include <sundials/sundials_math.h> -- --/* Constants */ -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --#define DEBUG -- --/* -- * ================================================================= -- * PROTOTYPES FOR PRIVATE FUNCTIONS -- * ================================================================= -- */ -- -- --/* help functions */ --void DenSLU2Vector(N_Vector res, SuperMatrix *x); --void Vector2DenSLU(N_Vector x, SuperMatrix *res); -- -- --/* KINPinv linit, lsetup, lsolve, and lfree routines */ -- --static int kinSLUGInit(KINMem kin_mem); --static int kinSLUGSetup(KINMem kin_mem); --static int kinSLUGSolve(KINMem kin_mem, N_Vector x, N_Vector b, -- realtype *res_norm); --static void kinSLUGFree(KINMem kin_mem); -- -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define func (kin_mem->kin_func) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#define sfdotJp (kin_mem->kin_sfdotJp) --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define spjac (kinpinv_mem->d_spjac) --#define J (kinpinv_mem->d_J) --#define spJ (kinpinv_mem->sp_J) --#define L (kinpinv_mem->sp_L) --#define U (kinpinv_mem->sp_U) --#define B (kinpinv_mem->sp_B) --#define perm_r (kinpinv_mem->sp_perm_r) --#define perm_c (kinpinv_mem->sp_perm_c) --#define spJTJ (kinpinv_mem->sp_JTJ) --#define options (kinpinv_mem->sp_options) --#define stat (kinpinv_mem->sp_stat) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define JTJ (kinpinv_mem->d_JTJ) --#define regularized (kinpinv_mem->d_regularized) --#define redojac (kinpinv_mem->d_redojac) --#define beta (kinpinv_mem->d_beta) --#define reg_param (kinpinv_mem->d_reg_param) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINSLUG -- * ----------------------------------------------------------------- -- * This routine initializes the memory record and sets various function -- * fields specific to the sparse regularization solver module. -- * KINSLUG sets the kin_linit, kin_lsetup, kin_lsolve, kin_lfree fields -- * in *kinmem to be kinSLUGInit, kinSLUGSetup, kinSLUGSolve, and -- * kinSLUGFree, respectively. -- * It allocates memory for a structure of type KINPinvMemRec and sets -- * the kin_lmem field in *kinmem to the address of this structure. -- * It sets setupNonNull in *kinmem to TRUE, and the djac field to the -- * default kinPinvDenseDQJac. -- * Finally, it allocates memory for J and pivots. -- * -- * NOTE: The sparse regularized linear solver assumes a serial implementation -- * of the NVECTOR package. Therefore, KINSLUG will first -- * test for compatible a compatible N_Vector internal -- * representation by checking that N_VGetArrayPointer and -- * N_VSetArrayPointer exist. -- * ----------------------------------------------------------------- -- */ -- --int KINSLUG(void *kinmem, int N) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- /* Check if kinmem is different from NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinv", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- /* Test if the NVECTOR package is present */ -- if (vec_tmpl->ops->nvgetarraypointer == NULL || -- vec_tmpl->ops->nvsetarraypointer == NULL) { -- KINProcessError(kin_mem, KINPINV_ILL_INPUT, "KINPINV", "KINPinv", MSGD_BAD_NVECTOR); -- return(KINPINV_ILL_INPUT); -- } -- -- if (lfree !=NULL) lfree(kin_mem); -- -- /* Set four main function fields in kin_mem */ -- linit = kinSLUGInit; -- lsetup = kinSLUGSetup; -- lsolve = kinSLUGSolve; -- lfree = kinSLUGFree; -- -- /* Get memory for KINDlsMemRec */ -- kinpinv_mem = NULL; -- kinpinv_mem = (KINPinvMem) malloc(sizeof(struct KINPinvMemRec)); -- if (kinpinv_mem == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- return(KINPINV_MEM_FAIL); -- } -- -- /* Set matrix type */ -- mtype = SUNDIALS_DENSE; -- -- /* Set default Jacobian routine and Jacobian data */ -- jacDQ = TRUE; -- spjac = NULL; -- J_data = NULL; -- last_flag = KINPINV_SUCCESS; -- -- setupNonNull = TRUE; -- -- /* Set problem dimension */ -- n = N; -- -- /* Allocate memory for J,JTJ and pivots */ -- -- spJ = NULL; -- spJ = calloc(1,sizeof(SuperMatrix)) ; -- /*spJ = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );*/ -- if (spJ == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- L = NULL; -- L = calloc(1,sizeof(SuperMatrix)) ; -- if (L == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- U = NULL; -- U = calloc(1,sizeof(SuperMatrix)) ; -- if (U == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- B = NULL; -- B = calloc(1,sizeof(SuperMatrix)) ; -- if (B == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- perm_r = NULL; -- perm_r = intMalloc(N); -- if (perm_r == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- perm_c = NULL; -- perm_c = calloc(N,sizeof(int)); -- if (perm_c == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- options = NULL; -- options = (superlu_options_t *) calloc(1,sizeof(superlu_options_t)); -- if (options == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(perm_c); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- stat = NULL; -- stat = (SuperLUStat_t *) calloc(1,sizeof(SuperLUStat_t)); -- if (stat == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(perm_c); -- free(options); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- spJTJ = NULL; -- spJTJ= (SuperMatrix *) calloc(1,sizeof(SuperMatrix)); -- if (spJTJ == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(perm_c); -- free(options); -- free(stat); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- if (printfl>0) printf("Sparse data structures created \n"); -- -- /* Set the default input options. */ -- set_default_options(options); -- options->ColPerm = NATURAL; -- if (printfl>0) printf("Sparse options set \n"); -- -- /* Initialize the statistics variables. */ -- StatInit(stat); -- if (printfl>0) printf("Sparse stat initialized \n"); -- -- /* This is a direct linear solver */ -- inexact_ls = FALSE; -- -- /* Attach linear solver memory to integrator memory */ -- lmem = kinpinv_mem; -- -- /* Set reg_param 'not set' */ -- reg_param = 0; -- -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ================================================================= -- * PRIVATE FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGInit -- * ----------------------------------------------------------------- -- * This routine does remaining initializations specific to the Sparse -- * linear solver. -- * ----------------------------------------------------------------- -- */ -- --static int kinSLUGInit(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- nje = 0; -- nfeDQ = 0; -- -- /* -- * Set where to find jacobian data -- * if jacDQ=True it will be calculated by finite differenses -- * otherwise a user-supplied jacobian will be used -- */ -- -- -- if (jacDQ) { -- djac = kinPinvDQJac; -- J_data = kin_mem; -- } else { -- J_data = kin_mem->kin_user_data; -- } -- -- /* Set regularization parameter */ -- if (reg_param == 0){ -- reg_param = 1; -- } -- last_flag = KINPINV_SUCCESS; -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGSetup -- * ----------------------------------------------------------------- -- * This routine does the setup operations for the linear solver and -- * prepares J transpose J + h^2 I if necessary for regularization. -- * ----------------------------------------------------------------- -- */ -- --static int kinSLUGSetup(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- SuperMatrix spJC ; /* Matrix postmultiplied by P_c */ -- int retval,relax,panel_size,*etree,permc_spec; -- double rp = 0; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- regularized = FALSE; -- -- /* Calculate value of jacobian */ -- nje++; -- -- retval = spjac(n, uu, fval, spJ, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- last_flag = -1; -- return(-1); -- } -- -- /* Set flag that the matrix needs to be factorized */ -- -- options->Fact = DOFACT; -- -- -- /* Get column permutation vector perm_c */ -- permc_spec = options->ColPerm; -- get_perm_c(permc_spec,spJ,perm_c); -- -- /* Permute Jacobian and calculate and calculate etree */ -- etree = calloc(n,sizeof(int)); -- -- sp_preorder(options, spJ, perm_c, etree, &spJC); -- -- /* Get sizes for supernodes and panels from SUperLU */ -- panel_size = sp_ienv(1); -- relax = sp_ienv(2); -- -- -- /* Try to do a LU factorization of spJ */ -- dgstrf(options, &spJC, relax, panel_size, etree, -- NULL, 0,perm_c, perm_r, L, U, stat, &retval); -- if (retval != 0){ -- /* Factorization failed, determine reason */ -- if (retval < 0) { -- printf("Argument %i has illegal value\n",(-retval)); -- } else if (retval > n) { -- printf("%i nb of bytes allocated when memory allocation failed \n",(retval+n)); -- } else { -- /* Jacobian is singular, perform regularization */ -- if (printfl > 0) printf("Jacobian is singular at %i, regularizing \n ",retval); -- -- -- /* Recalculate jacobian */ -- -- Destroy_CompCol_Matrix(spJ); -- nje++; -- retval = spjac(n, uu, fval, spJ, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- last_flag = -1; -- free(etree); -- return(-1); -- } -- -- /* Calculate regularization parameter */ -- /* Put rhs in B2 */ -- B = calloc(1,sizeof(SuperMatrix)); -- Vector2DenSLU(fval,B); -- -- rp = getRegParam(spJ,B); -- -- printf("rp: %f \n",rp*rp); -- if (rp > 1) { -- reg_param = 1; -- } else { -- reg_param = rp; -- } -- if (printfl > 0) printf("Regparam: %f, \n ",reg_param); -- /* Calculate regularized matrix */ -- spJTJ = regSparseMatrix(spJ,sqrt(reg_param)); -- -- /* Get column permutation vector perm_c */ -- permc_spec = options->ColPerm; -- get_perm_c(permc_spec,spJTJ,perm_c); -- -- /* Do a preorder on the column etree */ -- sp_preorder(options, spJTJ, perm_c, etree, &spJC); -- -- /* Get sizes for supernodes and panels from SUperLU */ -- panel_size = sp_ienv(1); -- relax = sp_ienv(2); -- -- -- /* Try to do a LU factorization of spJ */ -- dgstrf(options, &spJC, relax, panel_size, etree, -- NULL, 0,perm_c, perm_r, L, U, stat, &retval); -- -- regularized = TRUE; -- if (retval == 0) { -- if (printfl > 0) printf("Jacobian regularized \n"); -- /* set flag that the matrix is factorized */ -- options->Fact = FACTORED; -- } else { -- if (printfl >0) { -- if (retval > n) { -- printf("Error in Jacobian regularization, %d nb of bytes allocated when allocation failed. n: %d",retval,n); -- } else { -- printf("Error in Jacobian regularization, retval: %d, n: %d,reg_param: %f",retval,n,reg_param); -- } -- } -- } -- } -- -- -- } else { -- /* set flag that the matrix is factorized */ -- options->Fact = FACTORED; -- fflush(stdout); -- } -- -- redojac = FALSE; -- /* Free allocated memory and return flags */ -- SUPERLU_FREE(etree); -- Destroy_CompCol_Permuted(&spJC); -- -- last_flag = retval; -- if (retval > 0) return(-1); -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * Vector2SDenSLU -- * ----------------------------------------------------------------- -- * Routine used to get the information from a NVector into the structure -- * needed by SuperLU. -- * ----------------------------------------------------------------- -- */ --void Vector2DenSLU(N_Vector x, SuperMatrix *res) { -- realtype *data; -- int length = NV_LENGTH_S(x); -- Stype_t S = SLU_DN ; -- Dtype_t D = SLU_D ; -- Mtype_t M = SLU_GE; -- -- /* Get data from the N_Vector */ -- data = N_VGetArrayPointer(x); -- -- /* Create dense Supermatrix */ -- dCreate_Dense_Matrix(res,length,1,data,length,S,D,M); -- --} -- --/* -- * ----------------------------------------------------------------- -- * DenSLU2Vector -- * ----------------------------------------------------------------- -- * Routine used to get the information from a SuperMatrix into a -- * N_Vector. -- * ----------------------------------------------------------------- -- */ --void DenSLU2Vector(N_Vector res, SuperMatrix *x) { -- realtype *data; -- DNformat *storage; -- -- /* Get data pointer from SuperMatrix */ -- storage = x->Store; -- data = storage->nzval; -- -- /* Set data pointer in N_Vector */ -- N_VSetArrayPointer(data,res); -- --} -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGSolve -- * ----------------------------------------------------------------- -- * This routine handles the solve operation for the sparse linear solver -- * by calling the sparse backsolve routine. The returned value is 0. -- * ----------------------------------------------------------------- -- */ -- --static int kinSLUGSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm) --{ -- KINPinvMem kinpinv_mem; -- SuperMatrix *B2; -- int info; -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (redojac) { -- return 1; -- } -- /* Copy the right-hand side into x*/ -- N_VScale(ONE, b, x); -- -- -- if (regularized) { -- -- /* Put rhs in B2 */ -- B2 = calloc(1,sizeof(SuperMatrix)); -- Vector2DenSLU(x,B2); -- -- /* Get regularized rhs */ -- B = getRegRHS(spJ,B2); -- -- regularized = FALSE; -- redojac = TRUE; -- -- /* Take out the garbage */ -- free(B2); -- } else { -- /* put rhs in B */ -- Vector2DenSLU(x,B); -- -- } -- -- /* Back-solve, get solution in B and propagate it to x*/ -- dgstrs (NOTRANS, L, U, perm_c, perm_r, B, stat, &info); -- -- if (info < 0) { -- printf("Solving failed with flag %i \n",info); -- last_flag = info; -- -- } else { -- DenSLU2Vector(x,B); -- } -- -- -- -- -- /* Compute the terms Jpnorm and sfdotJp for use in the global strategy -- routines and in KINForcingTerm. Both of these terms are subsequently -- corrected if the step is reduced by constraints or the line search. -- -- sJpnorm is the norm of the scaled product (scaled by fscale) of -- the current Jacobian matrix J and the step vector p. -- -- sfdotJp is the dot product of the scaled f vector and the scaled -- vector J*p, where the scaling uses fscale. */ -- -- sJpnorm = N_VWL2Norm(b,fscale); -- N_VProd(b, fscale, b); -- N_VProd(b, fscale, b); -- sfdotJp = N_VDotProd(fval, b); -- -- -- last_flag = KINPINV_SUCCESS; -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGFree -- * ----------------------------------------------------------------- -- * This routine frees memory specific to the dense linear solver. -- * ----------------------------------------------------------------- -- */ -- --static void kinSLUGFree(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- -- kinpinv_mem = (KINPinvMem) lmem; -- /* De-allocate storage */ -- -- SUPERLU_FREE (perm_r); -- SUPERLU_FREE (perm_c); -- Destroy_CompCol_Matrix(spJ); -- Destroy_Dense_Matrix(B); -- Destroy_SuperNode_Matrix(L); -- Destroy_CompCol_Matrix(U); -- Destroy_CompCol_Matrix(spJTJ); -- StatFree(stat); -- free(stat); -- free(options); -- free(spJ); -- free(B); -- free(L); -- free(U); -- free(spJTJ); -- -- free(kinpinv_mem); kinpinv_mem = NULL; --} -- -Index: src/lib/kinsol_jmod.h -=================================================================== ---- assimulo/lib/kinsol_jmod.h (revision 838) -+++ assimulo/lib/kinsol_jmod.h (nonexistent) -@@ -1,241 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- -- --#ifndef _KINJMOD_H --#define _KINJMOD_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include <sundials/sundials_direct.h> --#include <sundials/sundials_nvector.h> -- --/* -- * ================================================================= -- * K I N P I N V C O N S T A N T S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPINV return values -- * ----------------------------------------------------------------- -- */ -- --#define KINPINV_SUCCESS 0 --#define KINPINV_MEM_NULL -1 --#define KINPINV_LMEM_NULL -2 --#define KINPINV_ILL_INPUT -3 --#define KINPINV_MEM_FAIL -4 -- --/* Additional last_flag values */ -- --#define KINPINV_JACFUNC_UNRECVR -5 --#define KINPINV_JACFUNC_RECVR -6 -- --/* -- * ================================================================= -- * F U N C T I O N T Y P E S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Type: KINPinvJacFn -- * ----------------------------------------------------------------- -- * -- * A dense Jacobian approximation function Jac must be of type -- * KINDlsDenseJacFn. Its parameters are: -- * -- * N - problem size. -- * -- * u - current iterate (unscaled) [input] -- * -- * fu - vector (type N_Vector) containing result of nonlinear -- * system function evaluated at current iterate: -- * fu = F(u) [input] -- * -- * J - dense matrix (of type DlsMat) that will be loaded -- * by a KINDlsDenseJacFn with an approximation to the -- * Jacobian matrix J = (dF_i/dy_j). -- * -- * user_data - pointer to user data - the same as the user_data -- * parameter passed to KINSetFdata. -- * -- * tmp1, tmp2 - available scratch vectors (volatile storage) -- * -- * A KINPinvJacFn should return 0 if successful, a positive -- * value if a recoverable error occurred, and a negative value if -- * an unrecoverable error occurred. -- * -- * ----------------------------------------------------------------- -- * -- * NOTE: The following are two efficient ways to load a dense Jac: -- * (1) (with macros - no explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = DENSE_COL(Jac,j); -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * (2) (without macros - explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = (Jac->data)[j]; -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * A third way, using the DENSE_ELEM(A,i,j) macro, is much less -- * efficient in general. It is only appropriate for use in small -- * problems in which efficiency of access is NOT a major concern. -- * -- * ----------------------------------------------------------------- -- */ -- -- --typedef int (*KINPinvJacFn)(int N, -- N_Vector u, N_Vector fu, -- DlsMat J, void *user_data, -- N_Vector tmp1, N_Vector tmp2); -- -- --/* -- * ================================================================= -- * E X P O R T E D F U N C T I O N S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Optional inputs to the KINPinv linear solver -- * ----------------------------------------------------------------- -- * -- * KINDlsSetDenseJacFn specifies the dense Jacobian approximation -- * routine to be used for a direct dense linear solver. -- * -- * By default, a difference quotient approximation, supplied with -- * the solver is used. -- * -- * The return value is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- --SUNDIALS_EXPORT int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac); -- --/* -- Set regularization parameter -- */ --SUNDIALS_EXPORT int KINPinvSetRegParam(void *kinmem, realtype reg_p); -- --/* -- * ----------------------------------------------------------------- -- * Optional outputs from a KINDLS linear solver -- * ----------------------------------------------------------------- -- * -- * KINPinvGetWorkSpace returns the real and integer workspace used -- * by the KINDLS linear solver. -- * KINPinvGetNumJacEvals returns the number of calls made to the -- * Jacobian evaluation routine. -- * KINPinvGetNumFuncEvals returns the number of calls to the user's F -- * routine due to finite difference Jacobian -- * evaluation. -- * KINPinvGetLastFlag returns the last error flag set by any of -- * the KINDLS interface functions. -- * KINPinvGetReturnFlagName returns the name of the constant -- * associated with a KINDLS return flag -- * -- * The return value of KINPinvGet* is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- -- --SUNDIALS_EXPORT int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB); --SUNDIALS_EXPORT int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB); --SUNDIALS_EXPORT int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB); --SUNDIALS_EXPORT int KINPinvGetLastFlag(void *kinmem, int *flag); --SUNDIALS_EXPORT char *KINPinvGetReturnFlagName(int flag); -- --#ifdef __cplusplus --} --#endif -- --#endif -Index: src/lib/kinpinv.c -=================================================================== ---- assimulo/lib/kinpinv.c (revision 838) -+++ assimulo/lib/kinpinv.c (nonexistent) -@@ -1,578 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --/* -- *This is a linear solver using Tikhonov regularization -- */ --#include <stdio.h> --#include <stdlib.h> --#include <math.h> -- --#include "kinsol_jmod_impl.h" --#include "kinpinv.h" --#include "kinsol/kinsol_impl.h" -- --#include <sundials/sundials_direct.h> --#include <nvector/nvector_serial.h> --#include <sundials/sundials_math.h> -- --/* Constants */ -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --/* -- * ================================================================= -- * PROTOTYPES FOR PRIVATE FUNCTIONS -- * ================================================================= -- */ -- -- --/* help functions */ --void regMatrix(realtype **JTJ_c, realtype **jac, realtype h,int size); -- --/* KINPinv linit, lsetup, lsolve, and lfree routines */ -- --static int kinPinvInit(KINMem kin_mem); --static int kinPinvSetup(KINMem kin_mem); --static int kinPinvSolve(KINMem kin_mem, N_Vector x, N_Vector b, -- realtype *res_norm); --static void kinPinvFree(KINMem kin_mem); -- -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define func (kin_mem->kin_func) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#if SUNDIALS_26 -- #define sFdotJp (kin_mem->kin_sFdotJp) --#else -- #define sfdotJp (kin_mem->kin_sfdotJp) --#endif --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define J (kinpinv_mem->d_J) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define JTJ (kinpinv_mem->d_JTJ) --#define regularized (kinpinv_mem->d_regularized) -- --#define redojac (kinpinv_mem->d_redojac) --#define beta (kinpinv_mem->d_beta) --#define reg_param (kinpinv_mem->d_reg_param) -- --#define ihfun (kin_mem->kin_ihfun) --#define ih_data (kin_mem->kin_ih_data) --#define ehfun (kin_mem->kin_ehfun) --#define eh_data (kin_mem->kin_eh_data) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPinv -- * ----------------------------------------------------------------- -- * This routine initializes the memory record and sets various function -- * fields specific to the dense pseude-inverse linear solver module. -- * KINPinv sets the kin_linit, kin_lsetup, kin_lsolve, kin_lfree fields -- * in *kinmem to be kinPinvInit, kinPinvSetup, kinPinvSolve, and -- * kinPinvFree, respectively. -- * It allocates memory for a structure of type KINDlsMemRec and sets -- * the kin_lmem field in *kinmem to the address of this structure. -- * It sets setupNonNull in *kinmem to TRUE, and the djac field to the -- * default kinDlsDenseDQJac. -- * Finally, it allocates memory for J and pivots. -- * -- * NOTE: The dense pseudo-inverse linear solver assumes a serial implementation -- * of the NVECTOR package. Therefore, KINPinv will first -- * test for compatible a compatible N_Vector internal -- * representation by checking that N_VGetArrayPointer and -- * N_VSetArrayPointer exist. -- * ----------------------------------------------------------------- -- */ -- --int KINPinv(void *kinmem, int N) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- /* Check if kinmem is different from NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinv", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- /* Test if the NVECTOR package is present */ -- if (vec_tmpl->ops->nvgetarraypointer == NULL || -- vec_tmpl->ops->nvsetarraypointer == NULL) { -- KINProcessError(kin_mem, KINPINV_ILL_INPUT, "KINPINV", "KINPinv", MSGD_BAD_NVECTOR); -- return(KINPINV_ILL_INPUT); -- } -- -- if (lfree !=NULL) lfree(kin_mem); -- -- /* Set four main function fields in kin_mem */ -- linit = kinPinvInit; -- lsetup = kinPinvSetup; -- lsolve = kinPinvSolve; -- lfree = kinPinvFree; -- -- /* Get memory for KINDlsMemRec */ -- kinpinv_mem = NULL; -- kinpinv_mem = (KINPinvMem) malloc(sizeof(struct KINPinvMemRec)); -- if (kinpinv_mem == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- return(KINPINV_MEM_FAIL); -- } -- -- /* Set matrix type */ -- mtype = SUNDIALS_DENSE; -- -- /* Set default Jacobian routine and Jacobian data */ -- jacDQ = TRUE; -- djac = NULL; -- J_data = NULL; -- last_flag = KINPINV_SUCCESS; -- -- setupNonNull = TRUE; -- -- /* Set problem dimension */ -- n = N; -- -- /* Allocate memory for J,RTR and pivots */ -- -- J = NULL; -- J = NewDenseMat(N, N); -- if (J == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- JTJ = NULL; -- JTJ = NewDenseMat(n, n); -- if (JTJ == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- DestroyMat(J); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- pivots = NULL; -- pivots = NewIntArray(N); -- if (pivots == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- DestroyMat(J); -- DestroyMat(JTJ); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- beta = NULL; -- beta = NewRealArray(N); -- if (beta == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- DestroyMat(J); -- DestroyMat(JTJ); -- DestroyArray(pivots); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- /* This is a direct linear solver */ -- inexact_ls = FALSE; -- -- /* Attach linear solver memory to integrator memory */ -- lmem = kinpinv_mem; -- /* Set reg_param 'not set' */ -- reg_param = 0; -- redojac=1 ; -- regularized = FALSE; -- nje = 0; -- nfeDQ = 0; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ================================================================= -- * PRIVATE FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * kinPinvInit -- * ----------------------------------------------------------------- -- * This routine does remaining initializations specific to the dense -- * linear solver. -- * ----------------------------------------------------------------- -- */ -- --static int kinPinvInit(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- nje = 0; -- nfeDQ = 0; -- -- /* -- * Set where to find jacobian data -- * if jacDQ=True it will be calculated by finite differenses -- * otherwise a user-supplied jacobian will be used -- */ -- -- if (jacDQ) { -- djac = kinPinvDQJac; -- J_data = kin_mem; -- } else { -- J_data = kin_mem->kin_user_data; -- } -- -- /* Set regularization parameter */ -- if (reg_param == 0){ -- reg_param = 1; -- } -- -- last_flag = KINPINV_SUCCESS; -- return(0); --} -- -- -- --/* -- * Function that calculates the regularized matrix for a given h -- * The result is returned in res while h is the regularization parameter -- * and J is the current prblem jacobian. -- */ -- --void regMatrix(realtype **JTJ_c, realtype **jac, realtype h, int size) --{ -- int i,j,k; -- -- -- i = 0; -- j = 0; -- k = 0; -- for (i=0;i<size;i++) { -- for (j=0;j<size;j++){ -- -- /*Calculate value at RTR(i,j) */ -- JTJ_c[j][i] = 0; -- for (k=0;k<size;k++) JTJ_c[j][i] += jac[j][k]*jac[i][k]; -- -- /* add the regularization parameter on the diagonal */ -- if (i==j)JTJ_c[j][i] += h*h; -- } -- } --} -- --/* -- * ----------------------------------------------------------------- -- * kinPinvSetup -- * ----------------------------------------------------------------- -- * This routine does the setup operations for the linear solver and -- * prepares J transpose J + h^2 I if necessary for regularization. -- * ----------------------------------------------------------------- -- */ -- --static int kinPinvSetup(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- long int ier; -- int retval; -- int i,j; -- -- realtype **JTJ_c ; -- realtype **jac ; -- realtype rp ; -- double data; -- realtype *b; -- -- char msg[200]; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- /* Calculate value of jacobian */ -- nje++; -- SetToZero(J); -- retval = djac(n, uu, fval, J, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- ehfun(KIN_LSETUP_FAIL,"KINPINV", "kinPinvSetup", "Jacobian evaluation function failed", eh_data); -- last_flag = -1; -- return(-1); -- } -- -- /* Try to do a LU factorization of J */ -- ier = DenseGETRF(J, (long int*)pivots); -- -- /* If the LU factorization failed, perform regularization */ -- if (ier > 0) { -- /* Calculate value of jacobian */ -- SetToZero(J); -- -- /* SetToZero(pivots);*/ -- for (i=0; i<n; i++) pivots[i] = 0; -- nje++; -- retval = djac(n, uu, fval, J, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- ehfun(KIN_LSETUP_FAIL, "KINPINV", "kinPinvSetup", "Jacobian evaluation function failed", eh_data); -- last_flag = -1; -- return(-1); -- } -- -- /* Calculate J tranpose J */ -- SetToZero(JTJ); -- -- /* Calculate the regularization parameter */ -- -- jac = J->cols; -- b = N_VGetArrayPointer(fval); -- rp = 0; -- for (i = 0;i<n;i++) { -- data = 0; -- for(j = 0;j<n;j++) data += jac[i][j]*b[j]; -- rp += data*data ; -- } -- /*printf("rp: %f \n",rp);*/ -- if (sqrt(rp)<1) { -- reg_param = sqrt(rp); -- } else { -- reg_param = 1; -- } -- if (printfl > 0) { -- sprintf(msg, "Singular jacobian detected, regparam: %f",reg_param); -- ihfun("KINPINV", "kinPinvSetup", msg, ih_data); -- } -- /* calculate a regularized matrix */ -- JTJ_c = JTJ->cols; -- -- regMatrix(JTJ_c,jac,reg_param,n); -- -- /* LU-factorize the regularized matrix*/ -- ier = DenseGETRF(JTJ,(long int*)pivots); -- regularized = TRUE; -- } -- else { -- regularized = FALSE; -- } -- -- redojac = FALSE; -- /* Return 0 if the LU was complete; otherwise return -1 */ -- last_flag = ier; -- if (ier > 0) return(-1); -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinPinvSolve -- * ----------------------------------------------------------------- -- * This routine handles the solve operation for the dense linear solver -- * by calling the dense backsolve routine. The returned value is 0. -- * ----------------------------------------------------------------- -- */ -- --static int kinPinvSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm) --{ -- KINPinvMem kinpinv_mem; -- realtype **jac; -- realtype *xd; -- realtype *bx; -- int i,j; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (redojac) { -- return(1); -- } -- -- if (regularized) { -- if (printfl > 0) -- ihfun("KINPINV", "kinPinvSetup", "Solving regularized problem", ih_data); -- /* Calculate new right hand side b = J transpose * b */ -- bx = N_VGetArrayPointer(b); -- xd = N_VGetArrayPointer(x); -- jac = J->cols; -- i=0; -- j=0; -- for (i=0;i<n;i++){ -- xd[i] = 0; -- for (j=0;j<n;j++) xd[i] += jac[i][j]*bx[j]; -- } -- /* Back-solve and get solution in x */ -- -- /*N_VScale(ONE,x,b);*/ -- -- -- DenseGETRS(JTJ, (long int*)pivots, xd); -- -- /* Calculate a fresh jacobian - not really needed */ -- redojac = TRUE; -- } else { -- -- /* Copy the right-hand side into x */ -- -- N_VScale(ONE, b, x); -- xd = N_VGetArrayPointer(x); -- -- /* Back-solve and get solution in x */ -- DenseGETRS(J, (long int*)pivots, xd); -- -- } -- -- /* Compute the terms Jpnorm and sfdotJp for use in the global strategy -- routines and in KINForcingTerm. Both of these terms are subsequently -- corrected if the step is reduced by constraints or the line search. -- -- sJpnorm is the norm of the scaled product (scaled by fscale) of -- the current Jacobian matrix J and the step vector p. -- -- sfdotJp is the dot product of the scaled f vector and the scaled -- vector J*p, where the scaling uses fscale. */ -- -- -- sJpnorm = N_VWL2Norm(b,fscale); -- N_VProd(b, fscale, b); -- N_VProd(b, fscale, b); --#if SUNDIALS_26 -- sFdotJp = N_VDotProd(fval, b); --#else -- sfdotJp = N_VDotProd(fval, b); --#endif -- -- last_flag = KINPINV_SUCCESS; -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinPinvFree -- * ----------------------------------------------------------------- -- * This routine frees memory specific to the dense linear solver. -- * ----------------------------------------------------------------- -- */ -- --static void kinPinvFree(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- kinpinv_mem = (KINPinvMem) lmem; -- -- DestroyMat(J); -- DestroyMat(JTJ); -- DestroyArray(pivots); -- DestroyArray(beta); -- free(kinpinv_mem); kinpinv_mem = NULL; --} -- -Index: src/lib/sundials_kinsol_core_wSLU.pyx -=================================================================== ---- assimulo/lib/sundials_kinsol_core_wSLU.pyx (revision 838) -+++ assimulo/lib/sundials_kinsol_core_wSLU.pyx (nonexistent) -@@ -1,413 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# 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 Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" --from __future__ import division -- --import scipy.sparse as ss -- --include "sundials_kinsol_core_wSLU.pxi" # Includes fcts from other header files --include "sundials_kinsol_core.pxi" # Includes the constants (textual include) -- -- --cdef int kin_res(N_Vector xv, N_Vector fval, void *problem_data): -- """ -- Residual fct called by KINSOL -- """ -- -- cdef: -- int i -- ProblemData pData = <ProblemData>problem_data -- ndarray[realtype, ndim = 1, mode = 'c'] rhs # Container holding return from user fct -- -- ndarray x = nv2arr(xv) -- realtype* resptr = (<N_VectorContent_Serial>fval.content).data -- -- try: -- rhs = (<object>pData.RHS)(x) -- for i in range(pData.dim): -- -- resptr[i] = rhs[i] -- -- return KIN_SUCCESS -- except: -- return KIN_REC_ERR -- --cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian, -- void *problem_data, N_Vector tmp1, N_Vector tmp2): -- """ -- This method is used to connect the Assimulo.Problem.jac to the Sundials -- Jacobian function. -- """ -- cdef: -- ProblemData pData = <ProblemData>problem_data -- realtype* col_i=DENSE_COL(Jacobian,0) -- ndarray x = nv2arr(xv) -- -- try: -- -- jac=(<object>pData.JAC)(x) -- #This needs further investigations: -- #memcpy(Jacobian.data,<realtype*>jac.data, pData.memSizeJac) -- -- for i in range(Neq): -- col_i = DENSE_COL(Jacobian, i) -- for j in range(Neq): -- col_i[j] = jac[j,i] -- -- return KINDLS_SUCCESS -- except: -- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description) -- -- --cdef int kin_sp_jac(int Neq, N_Vector xv, N_Vector fval, SuperMatrix *Jacobian, -- void *problem_data, N_Vector tmp1, N_Vector tmp2): -- """ -- This method is used to connect the Assimulo.Problem.sparsejac to the KINSOL -- Jacobian function. -- """ -- cdef: -- ProblemData pData = <ProblemData>problem_data -- ndarray x = nv2arr(xv) -- -- int m,n,nnz,i -- realtype* nzval = NULL -- int* rowind = NULL -- int* colptr = NULL -- Stype_t S = SLU_NC -- Dtype_t D = SLU_D -- Mtype_t M = SLU_GE -- -- try: -- # Get jacobian and extract data -- jac=(<object>pData.JAC)(x) -- nnz = jac.nnz -- m = jac.shape[0] -- n = jac.shape[1] -- nzval = <realtype*>calloc(len(jac.data),sizeof(realtype)) -- rowind = <int*>calloc(len(jac.indices),sizeof(int)) -- colptr = <int*>calloc(len(jac.indptr),sizeof(int)) -- arr2realtype(jac.data,nzval,len(jac.data)) -- arr2int(jac.indices,rowind,len(jac.indices)) -- arr2int(jac.indptr,colptr,len(jac.indptr)) -- -- # Create the necessary matrix -- dCreate_CompCol_Matrix(Jacobian, m, n, nnz, nzval, rowind, colptr, S, D, M) -- -- -- return KINDLS_SUCCESS -- except: -- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description) -- -- --class KINError(Exception): -- """ -- Kinsol exception -- """ -- -- msg = {KIN_MEM_FAIL : 'Memory allocation failed.', -- KIN_MEM_NULL : 'KINSOL not properly created.', -- KIN_ILL_INPUT: 'There was an input error.', -- KIN_NO_MALLOC: 'Memory not allocated, call KINSOL_init(...)', -- KIN_LINESEARCH_NONCONV: 'Linesearch could not find a step big enough.', -- KIN_MAXITER_REACHED: 'Max number of iterations reached', -- KIN_MXNEWT_5X_EXCEEDED: 'Max size of newton step exceeded 5 times.', -- KIN_LINESEARCH_BCFAIL: 'Linesearch beta-test failed, probably because of too poor progress.', -- KIN_LINSOLV_NO_RECOVERY: 'psolve encountered an error, but preconditioner is current', -- KIN_LINIT_FAIL: 'Init of linear solver failed.', -- KIN_LSETUP_FAIL: 'pset (preconditioner setup fct) encountered an error', -- KIN_LSOLVE_FAIL: 'Error in either psolve or in linear solver routine.', -- KIN_SYSFUNC_FAIL: 'Call to RHS failed.', -- KIN_FIRST_SYSFUNC_ERR: 'Call to RHS failed on first call', -- KIN_REPTD_SYSFUNC_ERR: 'Call to RHS failed multiple times.', -- KIN_STEP_LT_STPTOL: 'Scaled step length too small. Either an approximate solution or a local minimum is reached. Check value of residual.'} -- value = 0 -- def __init__(self, value): -- self.value = value -- -- def __str__(self): -- try: -- return repr(self.msg[self.value]) -- except KeyError: -- return repr('KINSOL failed with flag %s.'%(self.value)) -- --cdef class KINSOL_wrap: -- -- cdef: -- void* solver -- ProblemData pData # A struct containing problem data -- #ProblemData *ppData # Pointer to above mentioned struct -- N_Vector x_cur,x_init, x_scale, f_scale, con_nv -- booleantype noInitSetup, sparse -- realtype norm_of_res,reg_param -- ndarray con -- int print_level -- ndarray fscale -- -- def __cinit__(self): -- self.pData = ProblemData() #Create a new problem struct -- #self.ppData = &self.pData -- -- cpdef KINSOL_set_problem_info(self, RHS, dim, JAC= None): -- """ -- Sets the problem info to the desired struct -- """ -- -- # Sets residual function -- self.pData.RHS = <void*>RHS -- self.pData.dim = dim -- if JAC != None: -- self.pData.JAC = <void*>JAC -- else: -- self.pData.JAC = NULL -- -- -- def KINSOL_init(self,RHS,x0,dim, JAC = None, con = None,sparse = False, print_level = 0,norm_of_res = 0,reg_param=0,fscale=None): -- """ -- Initializes solver -- """ -- cdef int flag -- # Set problem info -- self.noInitSetup = False -- self.print_level = print_level -- self.norm_of_res = norm_of_res -- self.reg_param = reg_param -- self.KINSOL_set_problem_info(RHS,dim,JAC) -- self.con = con -- self.sparse = sparse -- self.fscale = fscale -- -- -- # Create initial guess from the supplied numpy array -- # print "x0 got in KINSOL wrapper: ", x0 -- self.x_cur = arr2nv(x0) -- self.x_init = arr2nv(x0) -- -- -- if self.solver == NULL: # solver runs for the first time -- -- # Create solver -- self.solver = KINCreate() -- if self.solver == NULL: -- raise KINError(KIN_MEM_FAIL) -- if self.print_level >0: -- print "KINSOL created" -- -- # Stop solver from performing precond setup -- flag = KINSetNoInitSetup(self.solver, self.noInitSetup) -- if flag < 0: -- raise KINError(flag) -- -- flag = KINSetPrintLevel(self.solver, self.print_level) -- if flag < 0: -- raise KINError(flag) -- -- if self.norm_of_res >0: -- flag = KINSetMaxNewtonStep(self.solver, self.norm_of_res) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetMaxNewtonStep(self.solver, 1e20) -- if flag < 0: -- raise KINError(flag) -- -- # If the user has specified constraints, connect them -- if self.con != None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(self.con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- -- # Allocate internal memory -- flag = KINInit(self.solver,kin_res,self.x_cur) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "KINSOL initialized" -- -- # Link to linear solver -- if self.sparse: -- flag = KINSLUG(self.solver,self.pData.dim) -- else: -- flag = KINPinv(self.solver,self.pData.dim) -- #flag = KINDense(self.solver,self.pData.dim) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Linear solver connected" -- -- # Set regularization parameter, if necessary -- if self.reg_param != 0.0: -- if self.sparse: -- flag = KINSLUGSetRegParam(self.solver, self.reg_param) -- else: -- flag = KINPinvSetRegParam(self.solver, self.reg_param) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Reg param ", self.reg_param, " set" -- else: -- # If the user has specified constraints, connect them -- if con != None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetConstraints(self.solver, NULL) -- if flag < 0: -- raise KINError(flag) -- -- # If the user supplied a Jacobien, link it to the solver -- if self.pData.JAC != NULL: -- if sparse: -- flag = KINSLUGSetJacFn(self.solver,kin_sp_jac) -- else: -- flag = KINPinvSetJacFn(self.solver,kin_jac) -- #flag = KINDlsSetDenseJacFn(self.solver,kin_jac) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Jacobian supplied by user connected" -- else: -- if sparse: -- flag = KINSLUGSetJacFn(self.solver,NULL) -- else: -- flag = KINPinvSetJacFn(self.solver,NULL) -- #flag = KINDlsSetDenseJacFn(self.solver,NULL) -- -- if flag < 0: -- raise KINError(flag) -- -- # Link the solver to the supplied problem data -- flag = KINSetUserData(self.solver,<void*>self.pData) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "User data set" -- -- # Create scaling vectors filled with ones -- # since the problem is assumed to be scaled -- self.x_scale = arr2nv(np.ones(self.pData.dim)) -- if self.fscale != None: -- self.f_scale = arr2nv(self.fscale) -- else: -- self.f_scale = arr2nv(np.ones(self.pData.dim)) -- -- -- def KINSOL_solve(self,wo_LineSearch): -- """ -- Function that should be called after a call to KINSOL_init -- solves the function supplied as RHS -- """ -- -- cdef: -- int flag, DLSflag, lsflag, count,ind -- long int nb_iters -- ndarray x0 -- if self.print_level >0: -- print "Calling solver..." -- -- if not wo_LineSearch: -- # Use LineSearch -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_LINESEARCH,self.x_scale,self.f_scale) -- count = 0 -- if flag > 1: -- # Probably stuck at a minimum -- self.x_cur = self.x_init -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- -- if flag > 1: -- # Stuck at minimum again -- raise KINError(flag) -- else: -- # Don't use linesearch -- flag = -7 -- -- while flag <0: -- # Get nbr of non-linear iterations -- lsflag = KINGetNumNonlinSolvIters(<void*>self.solver, &nb_iters) -- if lsflag != 0: -- raise KINError(lsflag) -- -- -- if (flag == -5 and nb_iters == 1): -- # Cannot take a step bigger than zero on the first iteration -- # Try with heuristic initial guess -- raise KINError(42) -- -- if (flag == -7 or flag == -6 or (flag == -5 and nb_iters > 1)): -- if count == 0: -- # First time -- print "Calling solver again as simple newton." -- self.x_cur = self.x_init -- DLSflag = KINSetMaxSetupCalls(self.solver, 1); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 1: -- # second time -- print "Turning of residual monitoring" -- DLSflag = KINSetNoResMon(self.solver, True); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 3: -- # fourth time, release the bound on the step -- -- print "Remove limit on stepsize" -- DLSflag = KINSetMaxNewtonStep(self.solver, 1.0e34) -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- count += 1 -- if count >= 10: -- if self.print_level >0: -- print "Tried with simple newton ten times" -- raise KINError(-7) -- -- else: -- -- print "Error in solve, flag: ", flag -- lsflag = KINPinvGetLastFlag(self.solver, &lsflag) -- #lsflag = KINDlsGetLastFlag(self.solver, &lsflag) -- if lsflag != 0: -- if lsflag <0: -- print "Last flag from Linear solver: ", lsflag -- else: -- print "Jacobian singular at column ", lsflag -- raise KINError(flag) -- -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- if self.print_level >0: -- print "Problem solved, returning result" -- -- return nv2arr(self.x_cur) -Index: src/lib/kinslug.h -=================================================================== ---- assimulo/lib/kinslug.h (revision 838) -+++ assimulo/lib/kinslug.h (nonexistent) -@@ -1,58 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB / Copyright (c) 2002, The Regents of the University of California. -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * 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 Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- * This file is a modification of the file kinsol_dense.h , revision 1.5, -- * from the SUNDIALS suite. The file is modified by: -- * -- * Johan Ylikiiskilä - johan.ylikiiskila@gmail.com -- * -- */ -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#ifndef _KINSLU_H --#define _KINSLU_H -- --#include "kinsol_jmod_wSLU.h" --#include <sundials/sundials_dense.h> -- --/* -- * ----------------------------------------------------------------- -- * Function : KINSLUG -- * ----------------------------------------------------------------- -- * A call to the KINSLUG function links the main solver with the -- * SLUG (SuperLU reGularization implementation) linear solver. -- * Its arguments are as follows: -- * -- * kinmem - pointer to an internal memory block allocated during a -- * prior call to KINCreate -- * -- * N - problem size -- * -- * The return value of KINSLUG is one of: -- * 0 if successful -- * int different from zero otherwise -- * ----------------------------------------------------------------- -- */ -- -- SUNDIALS_EXPORT int KINSLUG(void *kinmem, int N); -- --#endif -- --#ifdef __cplusplus --} --#endif -Index: src/lib/kinsol_jmod_wSLU.c -=================================================================== ---- assimulo/lib/kinsol_jmod_wSLU.c (revision 838) -+++ assimulo/lib/kinsol_jmod_wSLU.c (nonexistent) -@@ -1,530 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --/* -- * ================================================================= -- * IMPORTED HEADER FILES -- * ================================================================= -- */ -- -- --#include <stdio.h> --#include <stdlib.h> -- --#include "kinsol/kinsol_impl.h" --#include "kinsol_jmod_impl_wSLU.h" --#include <sundials/sundials_math.h> -- --/* -- * ================================================================= -- * FUNCTION SPECIFIC CONSTANTS -- * ================================================================= -- */ -- --/* Constant for DQ Jacobian approximation */ --#define MIN_INC_MULT RCONST(1000.0) -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define uround (kin_mem->kin_uround) --#define func (kin_mem->kin_func) --#define user_data (kin_mem->kin_user_data) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#define sfdotJp (kin_mem->kin_sfdotJp) --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define spjac (kinpinv_mem->d_spjac) --#define bjac (kinpinv_mem->d_bjac) --#define J (kinpinv_mem->d_J) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define RTR (kinpinv_mem->d_RTR) --#define reg_param (kinpinv_mem->d_reg_param) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINDlsSetJacFn -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (jac != NULL) { -- jacDQ = FALSE; -- djac = jac; -- } else { -- jacDQ = TRUE; -- } -- -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvSetRegParam -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetRegParam(void *kinmem, realtype reg_p) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- reg_param = reg_p; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINSLUGSetRegParam -- * ----------------------------------------------------------------- -- */ -- --int KINSLUGSetRegParam(void *kinmem, realtype reg_p) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- reg_param = reg_p; -- -- return(KINPINV_SUCCESS); --} -- --int KINSLUGSetJacFn(void *kinmem, KINSLUGJacFn jac) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINSLUG", "KINSLUGSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINSLUG", "KINSLUGSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (jac != NULL) { -- jacDQ = FALSE; -- spjac = jac; -- } else { -- jacDQ = TRUE; -- } -- -- return(KINPINV_SUCCESS); --} --/* -- * ----------------------------------------------------------------- -- * KINDlsGetWorkSpace -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetWorkSpace(void *kinmem, long int *lenrwLS, long int *leniwLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (mtype == SUNDIALS_DENSE) { -- *lenrwLS = n*n; -- *leniwLS = n; -- } else if (mtype == SUNDIALS_BAND) { -- *lenrwLS = n*(smu + mu + 2*ml + 2); -- *leniwLS = n; -- } -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumJacEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumJacEvals(void *kinmem, long int *njevals) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *njevals = nje; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumFuncEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumFuncEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumGuncEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *nfevalsLS = nfeDQ; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetLastFlag -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetLastFlag(void *kinmem, int *flag) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *flag = last_flag; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetReturnFlagName -- * ----------------------------------------------------------------- -- */ -- --char *KINPinvGetReturnFlagName(int flag) --{ -- char *name; -- -- name = (char *)malloc(30*sizeof(char)); -- -- switch(flag) { -- case KINPINV_SUCCESS: -- sprintf(name, "KINPINV_SUCCESS"); -- break; -- case KINPINV_MEM_NULL: -- sprintf(name, "KINPINV_MEM_NULL"); -- break; -- case KINPINV_LMEM_NULL: -- sprintf(name, "KINPINV_LMEM_NULL"); -- break; -- case KINPINV_ILL_INPUT: -- sprintf(name, "KINPINV_ILL_INPUT"); -- break; -- case KINPINV_MEM_FAIL: -- sprintf(name, "KINPINV_MEM_FAIL"); -- break; -- default: -- sprintf(name, "NONE"); -- } -- -- return(name); --} -- --/* -- * ================================================================= -- * DQ JACOBIAN APPROXIMATIONS -- * ================================================================= -- */ -- -- -- -- --/* -- * ----------------------------------------------------------------- -- * kinPinvDQJac -- * ----------------------------------------------------------------- -- * This routine generates a dense difference quotient approximation to -- * the Jacobian of F(u). It assumes that a dense matrix of type -- * DlsMat is stored column-wise, and that elements within each column -- * are contiguous. The address of the jth column of J is obtained via -- * the macro DENSE_COL and this pointer is associated with an N_Vector -- * using the N_VGetArrayPointer/N_VSetArrayPointer functions. -- * Finally, the actual computation of the jth column of the Jacobian is -- * done with a call to N_VLinearSum. -- * -- * The increment used in the finite-difference approximation -- * J_ij = ( F_i(u+sigma_j * e_j) - F_i(u) ) / sigma_j -- * is -- * sigma_j = max{|u_j|, |1/uscale_j|} * sqrt(uround) -- * -- * Note: uscale_j = 1/typ(u_j) -- * -- * NOTE: Any type of failure of the system function her leads to an -- * unrecoverable failure of the Jacobian function and thus -- * of the linear solver setup function, stopping KINSOL. -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2) --{ -- realtype inc, inc_inv, ujsaved, ujscale, sign; -- realtype *tmp2_data, *u_data, *uscale_data; -- N_Vector ftemp, jthCol; -- long int j; -- int retval = 0; -- -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* data points to kin_mem */ -- kin_mem = (KINMem) data; -- kinpinv_mem = (KINPinvMem) lmem; -- -- /* Save pointer to the array in tmp2 */ -- tmp2_data = N_VGetArrayPointer(tmp2); -- -- /* Rename work vectors for readibility */ -- ftemp = tmp1; -- jthCol = tmp2; -- -- /* Obtain pointers to the data for u and uscale */ -- u_data = N_VGetArrayPointer(u); -- uscale_data = N_VGetArrayPointer(uscale); -- -- /* This is the only for loop for 0..N-1 in KINSOL */ -- -- for (j = 0; j < N; j++) { -- -- /* Generate the jth col of Jac(u) */ -- -- N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol); -- -- ujsaved = u_data[j]; -- ujscale = ONE/uscale_data[j]; -- sign = (ujsaved >= ZERO) ? ONE : -ONE; --#if SUNDIALS_26 -- inc = sqrt_relfunc*SUNMAX(SUNRabs(ujsaved), ujscale)*sign; --#else -- inc = sqrt_relfunc*MAX(ABS(ujsaved), ujscale)*sign; --#endif -- u_data[j] += inc; -- -- retval = func(u, ftemp, user_data); -- nfeDQ++; -- if (retval != 0) break; -- -- u_data[j] = ujsaved; -- -- inc_inv = ONE/inc; -- N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol); -- -- } -- -- /* Restore original array pointer in tmp2 */ -- N_VSetArrayPointer(tmp2_data, tmp2); -- -- return(retval); --} -Index: src/lib/kinpinv.h -=================================================================== ---- assimulo/lib/kinpinv.h (revision 838) -+++ assimulo/lib/kinpinv.h (nonexistent) -@@ -1,107 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#ifndef _PINV_H --#define _PINV_H -- --#include "kinsol_jmod.h" --#include <sundials/sundials_dense.h> -- --/* -- * ----------------------------------------------------------------- -- * Function : KINpinv -- * ----------------------------------------------------------------- -- * A call to the KINpinv function links the main solver with the -- * pinv linear solver. Its arguments are as follows: -- * -- * kinmem - pointer to an internal memory block allocated during a -- * prior call to KINCreate -- * -- * N - problem size -- * -- * The return value of KINpinv is one of: -- * 0 if successful -- * int different from zero otherwise -- * ----------------------------------------------------------------- -- */ -- -- SUNDIALS_EXPORT int KINPinv(void *kinmem, int N); -- --#endif -- --#ifdef __cplusplus --} --#endif -Index: src/lib/kinsol_jmod_wSLU.h -=================================================================== ---- assimulo/lib/kinsol_jmod_wSLU.h (revision 838) -+++ assimulo/lib/kinsol_jmod_wSLU.h (nonexistent) -@@ -1,255 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifndef _KINJMOD_H --#define _KINJMOD_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include <sundials/sundials_direct.h> --#include <sundials/sundials_nvector.h> -- --#include "slu_ddefs.h" -- --/* -- * ================================================================= -- * K I N P I N V C O N S T A N T S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPINV return values -- * ----------------------------------------------------------------- -- */ -- --#define KINPINV_SUCCESS 0 --#define KINPINV_MEM_NULL -1 --#define KINPINV_LMEM_NULL -2 --#define KINPINV_ILL_INPUT -3 --#define KINPINV_MEM_FAIL -4 -- --/* Additional last_flag values */ -- --#define KINPINV_JACFUNC_UNRECVR -5 --#define KINPINV_JACFUNC_RECVR -6 -- --/* -- * ================================================================= -- * F U N C T I O N T Y P E S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Type: KINPinvJacFn -- * ----------------------------------------------------------------- -- * -- * A dense Jacobian approximation function Jac must be of type -- * KINDlsDenseJacFn. Its parameters are: -- * -- * N - problem size. -- * -- * u - current iterate (unscaled) [input] -- * -- * fu - vector (type N_Vector) containing result of nonlinear -- * system function evaluated at current iterate: -- * fu = F(u) [input] -- * -- * J - dense matrix (of type DlsMat) that will be loaded -- * by a KINDlsDenseJacFn with an approximation to the -- * Jacobian matrix J = (dF_i/dy_j). -- * -- * user_data - pointer to user data - the same as the user_data -- * parameter passed to KINSetFdata. -- * -- * tmp1, tmp2 - available scratch vectors (volatile storage) -- * -- * A KINPinvJacFn should return 0 if successful, a positive -- * value if a recoverable error occurred, and a negative value if -- * an unrecoverable error occurred. -- * -- * ----------------------------------------------------------------- -- * -- * NOTE: The following are two efficient ways to load a dense Jac: -- * (1) (with macros - no explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = DENSE_COL(Jac,j); -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * (2) (without macros - explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = (Jac->data)[j]; -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * A third way, using the DENSE_ELEM(A,i,j) macro, is much less -- * efficient in general. It is only appropriate for use in small -- * problems in which efficiency of access is NOT a major concern. -- * -- * ----------------------------------------------------------------- -- * ----------------------------------------------------------------- -- * Type: KINSLUGJacFn -- * ----------------------------------------------------------------- -- * -- * Basically the same as above except that instead of passing a -- * DlsMat to store the Jacobian SuperMatrix is passed. -- * -- * ----------------------------------------------------------------- -- */ -- -- --typedef int (*KINPinvJacFn)(int N, -- N_Vector u, N_Vector fu, -- DlsMat J, void *user_data, -- N_Vector tmp1, N_Vector tmp2); -- --typedef int (*KINSLUGJacFn)(int N, -- N_Vector u, N_Vector fu, -- SuperMatrix *J, void *user_data, -- N_Vector tmp1, N_Vector tmp2); -- -- --/* -- * ================================================================= -- * E X P O R T E D F U N C T I O N S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Optional inputs to the KINPinv linear solver -- * ----------------------------------------------------------------- -- * -- * KINDlsSetDenseJacFn specifies the dense Jacobian approximation -- * routine to be used for a direct dense linear solver. -- * -- * By default, a difference quotient approximation, supplied with -- * the solver is used. -- * -- * The return value is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- --SUNDIALS_EXPORT int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac); -- --SUNDIALS_EXPORT int KINSLUGSetJacFn(void *kinmem, KINSLUGJacFn jac); -- --SUNDIALS_EXPORT int KINSLUGSetRegParam(void *kinmem, realtype reg_p); -- --SUNDIALS_EXPORT int KINPinvSetRegParam(void *kinmem, realtype reg_p); -- --/* -- * ----------------------------------------------------------------- -- * Optional outputs from a KINDLS linear solver -- * ----------------------------------------------------------------- -- * -- * KINPinvGetWorkSpace returns the real and integer workspace used -- * by the KINDLS linear solver. -- * KINPinvGetNumJacEvals returns the number of calls made to the -- * Jacobian evaluation routine. -- * KINPinvGetNumFuncEvals returns the number of calls to the user's F -- * routine due to finite difference Jacobian -- * evaluation. -- * KINPinvGetLastFlag returns the last error flag set by any of -- * the KINDLS interface functions. -- * KINPinvGetReturnFlagName returns the name of the constant -- * associated with a KINDLS return flag -- * -- * The return value of KINPinvGet* is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- --SUNDIALS_EXPORT int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB); --SUNDIALS_EXPORT int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB); --SUNDIALS_EXPORT int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB); --SUNDIALS_EXPORT int KINPinvGetLastFlag(void *kinmem, int *flag); --SUNDIALS_EXPORT char *KINPinvGetReturnFlagName(int flag); -- --#ifdef __cplusplus --} --#endif -- --#endif -Index: src/lib/sundials_kinsol_core.pxi -=================================================================== ---- assimulo/lib/sundials_kinsol_core.pxi (revision 838) -+++ assimulo/lib/sundials_kinsol_core.pxi (nonexistent) -@@ -1,144 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# 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 Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --This file contains the constants from Sundials. --""" --# ----------------------------------------------------------------- --# KINSOL return flags --# ----------------------------------------------------------------- -- --DEF KIN_SUCCESS =0 --DEF KIN_INITIAL_GUESS_OK =1 --DEF KIN_STEP_LT_STPTOL =2 -- --DEF KIN_REC_ERR =5 -- --DEF KIN_WARNING =99 -- --DEF KIN_MEM_NULL =-1 --DEF KIN_ILL_INPUT =-2 --DEF KIN_NO_MALLOC =-3 --DEF KIN_MEM_FAIL =-4 --DEF KIN_LINESEARCH_NONCONV =-5 --DEF KIN_MAXITER_REACHED =-6 --DEF KIN_MXNEWT_5X_EXCEEDED =-7 --DEF KIN_LINESEARCH_BCFAIL =-8 --DEF KIN_LINSOLV_NO_RECOVERY =-9 --DEF KIN_LINIT_FAIL =-10 --DEF KIN_LSETUP_FAIL =-11 --DEF KIN_LSOLVE_FAIL =-12 -- --DEF KIN_SYSFUNC_FAIL =-13 --DEF KIN_FIRST_SYSFUNC_ERR =-14 --DEF KIN_REPTD_SYSFUNC_ERR =-15 -- -- --# ----------------------------------------------------------------- --# Enumeration for inputs to KINSetEtaForm (eta choice) --# ----------------------------------------------------------------- --# KIN_ETACONSTANT : use constant value for eta (default value is --# 0.1 but a different value can be specified via --# a call to KINSetEtaConstValue) --# --# KIN_ETACHOICE1 : use choice #1 as given in Eisenstat and Walker's --# paper of SIAM J.Sci.Comput.,17 (1996), pp 16-32, --# wherein eta is defined to be: --# --# eta(k+1) = ABS(||F(u_k+1)||_L2-||F(u_k)+J(u_k)*p_k||_L2) --# --------------------------------------------- --# ||F(u_k)||_L2 --# --# 1+sqrt(5) --# eta_safe = eta(k)^ealpha where ealpha = --------- --# 2 --# --# KIN_ETACHOICE2 : use choice #2 as given in Eisenstat and Walker's --# paper wherein eta is defined to be: --# --# [ ||F(u_k+1)||_L2 ]^ealpha --# eta(k+1) = egamma * [ --------------- ] --# [ ||F(u_k)||_L2 ] --# --# where egamma = [0,1] and ealpha = (1,2] --# --# eta_safe = egamma*(eta(k)^ealpha) --# --# Note: The default values of the scalar --# coefficients egamma and ealpha (both required) --# are egamma = 0.9 and ealpha = 2.0, but the --# routine KINSetEtaParams can be used to specify --# different values. --# --# When using either KIN_ETACHOICE1 or KIN_ETACHOICE2, if --# eta_safe > 0.1 then the following safeguard is applied: --# --# eta(k+1) = MAX {eta(k+1), eta_safe} --# --# The following safeguards are always applied when using either --# KIN_ETACHOICE1 or KIN_ETACHOICE2 so that eta_min <= eta <= eta_max: --# --# eta(k+1) = MAX {eta(k+1), eta_min} --# eta(k+1) = MIN {eta(k+1), eta_max} --# --# where eta_min = 1.0e-4 and eta_max = 0.9 (see KINForcingTerm). --# ----------------------------------------------------------------- -- --DEF KIN_ETACHOICE1 =1 --DEF KIN_ETACHOICE2 =2 --DEF KIN_ETACONSTANT =3 -- --# ----------------------------------------------------------------- --# Enumeration for global strategy --# ----------------------------------------------------------------- --# Choices are KIN_NONE and KIN_LINESEARCH. --# ----------------------------------------------------------------- -- --DEF KIN_NONE =0 --DEF KIN_LINESEARCH =1 -- --# ----------------------------------------------------------------- --# KINDirect constants --# ----------------------------------------------------------------- -- --DEF KINDLS_SUCCESS =0 -- --DEF KINDLS_MEM_NULL =-1 --DEF KINDLS_LMEM_NULL =-2 --DEF KINDLS_ILL_INPUT =-3 --DEF KINDLS_MEM_FAIL =-4 --DEF KINDLS_JACFUNC_UNRECVR =-5 --DEF KINDLS_JACFUNC_RECVR =-6 -- -- --# ----------------------------------------------------------------- --# KINPinv constants --# ----------------------------------------------------------------- -- --DEF KINPINV_SUCCESS = 0 --DEF KINPINV_MEM_NULL = -1 --DEF KINPINV_LMEM_NULL = -2 --DEF KINPINV_ILL_INPUT = -3 --DEF KINPINV_MEM_FAIL = -4 --DEF KINPINV_JACFUNC_UNRECVR = -5 --DEF KINPINV_JACFUNC_RECVR = -6 -Index: src/lib/kinsol_jmod_impl.h -=================================================================== ---- assimulo/lib/kinsol_jmod_impl.h (revision 838) -+++ assimulo/lib/kinsol_jmod_impl.h (nonexistent) -@@ -1,153 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifndef _KINJMOD_IMPL_H --#define _KINJMOD_IMPL_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include "kinsol_jmod.h" -- -- --/* -- * ----------------------------------------------------------------- -- * Types: KINPinvMemRec, KINPinvMem -- * ----------------------------------------------------------------- -- * The type KINPinvMem is pointer to a KINPinvMemRec. -- * This structure contains KINPinv solver-specific data. -- * ----------------------------------------------------------------- -- */ -- --typedef struct KINPinvMemRec { -- -- int d_type; /* SUNDIALS_DENSE or SUNDIALS_BAND */ -- -- int d_n; /* problem dimension */ -- -- int d_ml; /* lower bandwidth of Jacobian */ -- int d_mu; /* upper bandwidth of Jacobian */ -- int d_smu; /* upper bandwith of M = MIN(N-1,d_mu+d_ml) */ -- -- booleantype d_jacDQ; /* TRUE if using internal DQ Jacobian approx. */ -- KINPinvJacFn d_djac; /* dense Jacobian routine to be called */ -- -- void *d_J_data; /* J_data is passed to djac or bjac */ -- -- DlsMat d_J; /* problem Jacobian */ -- -- int *d_pivots; /* pivot array for PM = LU */ -- realtype *d_beta; -- realtype d_reg_param; /* Regularization parameter */ -- long int d_nje; /* no. of calls to jac */ -- -- long int d_nfeDQ; /* no. of calls to F due to DQ Jacobian approx. */ -- -- int d_last_flag; /* last error return flag */ -- -- DlsMat d_JTJ; -- booleantype d_regularized; /* Flag if the regularized problem should be solved*/ -- booleantype d_redojac; -- --} *KINPinvMem; -- -- --/* -- * ----------------------------------------------------------------- -- * Prototypes of internal functions -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2); -- --/* -- * ----------------------------------------------------------------- -- * Error Messages -- * ----------------------------------------------------------------- -- */ -- --#define MSGD_KINMEM_NULL "KINSOL memory is NULL." --#define MSGD_BAD_NVECTOR "A required vector operation is not implemented." --#define MSGD_MEM_FAIL "A memory request failed." --#define MSGD_LMEM_NULL "Linear solver memory is NULL." --#define MSGD_BAD_SIZES "Illegal bandwidth parameter(s). Must have 0 <= ml, mu <= N-1." --#define MSGD_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner." -- --#ifdef __cplusplus --} --#endif -- --#endif |