diff options
Diffstat (limited to 'r845.patch')
-rw-r--r-- | r845.patch | 1689 |
1 files changed, 0 insertions, 1689 deletions
diff --git a/r845.patch b/r845.patch deleted file mode 100644 index dec23fc92e53..000000000000 --- a/r845.patch +++ /dev/null @@ -1,1689 +0,0 @@ -Index: setup.py -=================================================================== ---- setup.py (revision 844) -+++ setup.py (revision 845) -@@ -309,15 +309,30 @@ - self.with_SUNDIALS=True - L.debug('SUNDIALS found.') - sundials_version = None -+ sundials_vector_type_size = None - - try: - if os.path.exists(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')): - with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f: - for line in f: -- if "SUNDIALS_PACKAGE_VERSION" in line: -+ if "SUNDIALS_PACKAGE_VERSION" in line or "SUNDIALS_VERSION" in line: - sundials_version = tuple([int(f) for f in line.split()[-1][1:-1].split(".")]) - L.debug('SUNDIALS %d.%d found.'%(sundials_version[0], sundials_version[1])) -- except Exception: -+ break -+ with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f: -+ for line in f: -+ if "SUNDIALS_INT32_T" in line and line.startswith("#define"): -+ sundials_vector_type_size = "32" -+ L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size)) -+ break -+ if "SUNDIALS_INT64_T" in line and line.startswith("#define"): -+ sundials_vector_type_size = "64" -+ L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size)) -+ if self.with_SLU: -+ L.warning("It is recommended to set the SUNDIALS_INDEX_TYPE to an 32bit integer when using SUNDIALS together with SuperLU.") -+ L.warning("SuperLU may not function properly.") -+ break -+ except Exception as e: - if os.path.exists(os.path.join(os.path.join(self.incdirs,'arkode'), 'arkode.h')): #This was added in 2.6 - sundials_version = (2,6,0) - L.debug('SUNDIALS 2.6 found.') -@@ -326,6 +341,7 @@ - L.debug('SUNDIALS 2.5 found.') - - self.SUNDIALS_version = sundials_version -+ self.SUNDIALS_vector_size = sundials_vector_type_size - - else: - L.warning(("Could not find Sundials, check the provided path (--sundials-home={}) "+ -@@ -375,7 +391,8 @@ - # SUNDIALS - if self.with_SUNDIALS: - compile_time_env = {'SUNDIALS_VERSION': self.SUNDIALS_version, -- 'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu} -+ 'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu, -+ 'SUNDIALS_VECTOR_SIZE': self.SUNDIALS_vector_size} - #CVode and IDA - ext_list += cythonize(["assimulo" + os.path.sep + "solvers" + os.path.sep + "sundials.pyx"], - include_path=[".","assimulo","assimulo" + os.sep + "lib"], -@@ -382,11 +399,19 @@ - compile_time_env=compile_time_env, force=True) - ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs] - ext_list[-1].library_dirs = [self.libdirs] -- ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"] -+ -+ if self.SUNDIALS_version >= (3,0,0): -+ ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas", "sundials_sunlinsoldense", "sundials_sunlinsolspgmr", "sundials_sunmatrixdense", "sundials_sunmatrixsparse"] -+ else: -+ ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"] - if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU -+ if self.SUNDIALS_version >= (3,0,0): -+ ext_list[-1].libraries.extend(["sundials_sunlinsolsuperlumt"]) -+ - ext_list[-1].include_dirs.append(self.SLUincdir) - ext_list[-1].library_dirs.append(self.SLUlibdir) - ext_list[-1].libraries.extend(self.superLUFiles) -+ - - #Kinsol - ext_list += cythonize(["assimulo"+os.path.sep+"solvers"+os.path.sep+"kinsol.pyx"], -@@ -395,7 +420,11 @@ - ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs] - ext_list[-1].library_dirs = [self.libdirs] - ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial"] -- -+ -+ if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU -+ ext_list[-1].include_dirs.append(self.SLUincdir) -+ ext_list[-1].library_dirs.append(self.SLUlibdir) -+ ext_list[-1].libraries.extend(self.superLUFiles) - - for el in ext_list: - #Debug -Index: assimulo/solvers/__init__.py -=================================================================== ---- assimulo/solvers/__init__.py (revision 844) -+++ assimulo/solvers/__init__.py (revision 845) -@@ -25,7 +25,7 @@ - from .euler import ExplicitEuler
- from .euler import ImplicitEuler
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .radau5 import Radau5ODE
- from .radau5 import Radau5DAE
-@@ -32,43 +32,43 @@ - from .radau5 import _Radau5ODE
- from .radau5 import _Radau5DAE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .sundials import IDA
- from .sundials import CVode
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .kinsol import KINSOL
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .runge_kutta import RungeKutta34
- from .runge_kutta import RungeKutta4
- from .runge_kutta import Dopri5
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .rosenbrock import RodasODE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .odassl import ODASSL
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .odepack import LSODAR
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .radar5 import Radar5ODE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .dasp3 import DASP3ODE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .glimda import GLIMDA
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
-Index: src/solvers/kinsol.pyx -=================================================================== ---- assimulo/solvers/kinsol.pyx (revision 844) -+++ assimulo/solvers/kinsol.pyx (revision 845) -@@ -30,8 +30,8 @@ - cimport sundials_includes as SUNDIALS
-
- #Various C includes transfered to namespace
--from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL
--from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat
-+from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype
-+from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse
- from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial
- from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial
-
-@@ -57,6 +57,8 @@ -
- cdef object pt_fcn, pt_jac, pt_jacv, pt_prec_setup, pt_prec_solve
- cdef object _added_linear_solver
-+ cdef SUNDIALS.SUNMatrix sun_matrix
-+ cdef SUNDIALS.SUNLinearSolver sun_linearsolver
-
- def __init__(self, problem):
- Algebraic.__init__(self, problem) #Calls the base class
-@@ -86,6 +88,7 @@ - self.options["no_min_epsilon"] = False #Specifies wheter the scaled linear residual is bounded from below
- self.options["max_beta_fails"] = 10
- self.options["max_krylov"] = 0
-+ self.options["precond"] = PREC_NONE
-
- #Statistics
- self.statistics["nfevals"] = 0 #Function evaluations
-@@ -109,6 +112,13 @@ - if self.kinsol_mem != NULL:
- #Free Memory
- SUNDIALS.KINFree(&self.kinsol_mem)
-+
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ if self.sun_matrix != NULL:
-+ SUNDIALS.SUNMatDestroy(self.sun_matrix)
-+
-+ if self.sun_linearsolver != NULL:
-+ SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
-
- def update_variable_scaling(self, value="Automatic"):
- """
-@@ -170,7 +180,7 @@ -
- cdef initialize_kinsol(self):
- cdef int flag #Used for return
--
-+
- self.y_temp = arr2nv(self.y)
- self.y_scale = arr2nv([1.0]*self.problem_info["dim"])
- self.f_scale = arr2nv([1.0]*self.problem_info["dim"])
-@@ -208,16 +218,34 @@ -
- cpdef add_linear_solver(self):
- if self.options["linear_solver"] == "DENSE":
-- flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ #Create a dense Sundials matrix
-+ self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
-+ #Create a dense Sundials linear solver
-+ self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.y_temp, self.sun_matrix)
-+ #Attach it to Kinsol
-+ flag = SUNDIALS.KINDlsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver, self.sun_matrix)
-+ ELSE:
-+ flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
- if flag < 0:
- raise KINSOLError(flag)
-
- if self.problem_info["jac_fcn"]:
-- flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ flag = SUNDIALS.KINDlsSetJacFn(self.kinsol_mem, kin_jac);
-+ ELSE:
-+ flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
- if flag < 0:
- raise KINSOLError(flag)
- elif self.options["linear_solver"] == "SPGMR":
-- flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ #Create the linear solver
-+ self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.y_temp, self.options["precond"], self.options["max_krylov"])
-+ #Attach it to Kinsol
-+ flag = SUNDIALS.KINSpilsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver)
-+ ELSE:
-+ #Specify the use of KINSpgmr linear solver.
-+ flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
- if flag < 0:
- raise KINSOLError(flag)
-
-Index: src/solvers/sundials.pyx -=================================================================== ---- assimulo/solvers/sundials.pyx (revision 844) -+++ assimulo/solvers/sundials.pyx (revision 845) -@@ -34,8 +34,8 @@ - cimport sundials_includes as SUNDIALS - - #Various C includes transfered to namespace --from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL --from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat -+from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype -+from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse - from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial - from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial - -@@ -70,6 +70,8 @@ - cdef public N.ndarray yS0 - #cdef N.ndarray _event_info - cdef public N.ndarray g_old -+ cdef SUNDIALS.SUNMatrix sun_matrix -+ cdef SUNDIALS.SUNLinearSolver sun_linearsolver - - def __init__(self, problem): - Implicit_ODE.__init__(self, problem) #Calls the base class -@@ -100,7 +102,9 @@ - self.options["dqrhomax"] = 0.0 - self.options["pbar"] = [1]*self.problem_info["dimSens"] - self.options["external_event_detection"] = False #Sundials rootfinding is used for event location as default -+ self.options["precond"] = PREC_NONE - -+ - #Solver support - self.supports["report_continuously"] = True - self.supports["interpolated_output"] = True -@@ -173,6 +177,13 @@ - if self.ida_mem != NULL: - #Free Memory - SUNDIALS.IDAFree(&self.ida_mem) -+ -+ IF SUNDIALS_VERSION >= (3,0,0): -+ if self.sun_matrix != NULL: -+ SUNDIALS.SUNMatDestroy(self.sun_matrix) -+ -+ if self.sun_linearsolver != NULL: -+ SUNDIALS.SUNLinSolFree(self.sun_linearsolver) - - cpdef state_event_info(self): - """ -@@ -255,21 +266,30 @@ - if flag < 0: - raise IDAError(flag, self.t) - -- #Specify the use of the internal dense linear algebra functions. -- flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim) -- if flag < 0: -- raise IDAError(flag, self.t) -- -- #Choose a linear solver if and only if NEWTON is choosen -+ #Choose a linear solver if and only if NEWTON is choosen - if self.options["linear_solver"] == 'DENSE': -- #Specify the use of the internal dense linear algebra functions. -- flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create a dense Sundials matrix -+ self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim) -+ #Create a dense Sundials linear solver -+ self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix) -+ #Attach it to IDA -+ flag = SUNDIALS.IDADlsSetLinearSolver(self.ida_mem, self.sun_linearsolver, self.sun_matrix); -+ ELSE: -+ #Specify the use of the internal dense linear algebra functions. -+ flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim) - if flag < 0: - raise IDAError(flag, self.t) - - elif self.options["linear_solver"] == 'SPGMR': -- #Specify the use of SPGMR linear solver. -- flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create the linear solver -+ self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], 0) -+ #Attach it to IDAS -+ flag = SUNDIALS.IDASpilsSetLinearSolver(self.ida_mem, self.sun_linearsolver) -+ ELSE: -+ #Specify the use of SPGMR linear solver. -+ flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations - if flag < 0: - raise IDAError(flag, self.t) - -@@ -310,12 +330,17 @@ - if self.options["linear_solver"] == 'DENSE': - #Specify the jacobian to the solver - if self.pData.JAC != NULL and self.options["usejac"]: -- -- flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, ida_jac) -+ ELSE: -+ flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac) - if flag < 0: - raise IDAError(flag,self.t) - else: -- flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, NULL) -+ ELSE: -+ flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL) - if flag < 0: - raise IDAError(flag,self.t) - -@@ -322,11 +347,17 @@ - elif self.options["linear_solver"] == 'SPGMR': - #Specify the jacobian times vector function - if self.pData.JACV != NULL and self.options["usejac"]: -- flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv); -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, SUNDIALS.ida_spils_jtsetup_dummy, ida_jacv); -+ ELSE: -+ flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv); - if flag < 0: - raise IDAError(flag, self.t) - else: -- flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL); -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, NULL, NULL); -+ ELSE: -+ flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL); - if flag < 0: - raise IDAError(flag, self.t) - else: -@@ -1415,6 +1446,8 @@ - cdef public N.ndarray yS0 - #cdef N.ndarray _event_info - cdef public N.ndarray g_old -+ cdef SUNDIALS.SUNMatrix sun_matrix -+ cdef SUNDIALS.SUNLinearSolver sun_linearsolver - - def __init__(self, problem): - Explicit_ODE.__init__(self, problem) #Calls the base class -@@ -1479,6 +1512,13 @@ - if self.cvode_mem != NULL: - #Free Memory - SUNDIALS.CVodeFree(&self.cvode_mem) -+ -+ IF SUNDIALS_VERSION >= (3,0,0): -+ if self.sun_matrix != NULL: -+ SUNDIALS.SUNMatDestroy(self.sun_matrix) -+ -+ if self.sun_linearsolver != NULL: -+ SUNDIALS.SUNLinSolFree(self.sun_linearsolver) - - cpdef get_local_errors(self): - """ -@@ -2019,29 +2059,51 @@ - Updates the simulation options. - """ - cdef flag -- -+ - #Choose a linear solver if and only if NEWTON is choosen - if self.options["linear_solver"] == 'DENSE' and self.options["iter"] == "Newton": -- #Specify the use of the internal dense linear algebra functions. -- flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim) -- if flag < 0: -- raise CVodeError(flag) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create a dense Sundials matrix -+ self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim) -+ #Create a dense Sundials linear solver -+ self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix) -+ #Attach it to CVode -+ flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix); -+ if flag < 0: -+ raise CVodeError(flag) -+ ELSE: -+ #Specify the use of the internal dense linear algebra functions. -+ flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim) -+ if flag < 0: -+ raise CVodeError(flag) - - #Specify the jacobian to the solver - if self.pData.JAC != NULL and self.options["usejac"]: -- flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac); -+ ELSE: -+ flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac) - if flag < 0: - raise CVodeError(flag) - else: -- flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, NULL); -+ ELSE: -+ flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL) - if flag < 0: - raise CVodeError(flag) - - elif self.options["linear_solver"] == 'SPGMR' and self.options["iter"] == "Newton": -- #Specify the use of CVSPGMR linear solver. -- flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"]) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create the linear solver -+ self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], self.options["maxkrylov"]) -+ #Attach it to CVode -+ flag = SUNDIALS.CVSpilsSetLinearSolver(self.cvode_mem, self.sun_linearsolver) -+ ELSE: -+ #Specify the use of CVSPGMR linear solver. -+ flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"]) - if flag < 0: -- raise CVodeError(flag) -+ raise CVodeError(flag) - - if self.pData.PREC_SOLVE != NULL: - if self.pData.PREC_SETUP != NULL: -@@ -2055,11 +2117,17 @@ - - #Specify the jacobian times vector function - if self.pData.JACV != NULL and self.options["usejac"]: -- flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv) -- if flag < 0: -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, SUNDIALS.cv_spils_jtsetup_dummy, cv_jacv); -+ ELSE: -+ flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv) -+ if flag < 0: - raise CVodeError(flag) - else: -- flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, NULL, NULL); -+ ELSE: -+ flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL) - if flag < 0: - raise CVodeError(flag) - elif self.options["linear_solver"] == 'SPARSE' and self.options["iter"] == "Newton": -@@ -2066,17 +2134,28 @@ - - if SUNDIALS.version() < (2,6,0): - raise AssimuloException("Not supported with this SUNDIALS version.") -- -+ if SUNDIALS.with_superlu() == 0: -+ raise AssimuloException("No support for SuperLU was detected, please verify that SuperLU and SUNDIALS has been installed correctly.") -+ - #Specify the use of CVSPGMR linear solver. - if self.problem_info["jac_fcn_nnz"] == -1: - raise AssimuloException("Need to specify the number of non zero elements in the Jacobian via the option 'jac_nnz'") -- flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"]) -+ -+ IF SUNDIALS_VERSION >= (3,0,0): -+ self.sun_matrix = SUNDIALS.SUNSparseMatrix(self.pData.dim, self.pData.dim, self.problem_info["jac_fcn_nnz"], CSC_MAT) -+ self.sun_linearsolver = SUNDIALS.SUNSuperLUMT(self.yTemp, self.sun_matrix, self.options["num_threads"]) -+ flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix) -+ ELSE: -+ flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"]) - if flag < 0: - raise CVodeError(flag) - - #Specify the jacobian to the solver - if self.pData.JAC != NULL and self.options["usejac"]: -- flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac_sparse) -+ ELSE: -+ flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse) - if flag < 0: - raise CVodeError(flag) - else: -@@ -2899,7 +2978,10 @@ - flag = SUNDIALS.CVSpilsGetNumRhsEvals(self.cvode_mem, &nfevalsLS) #Number of rhs due to jac*vector - self.statistics["njacvecs"] += njvevals - elif self.options["linear_solver"] == "SPARSE": -- flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals) -+ ELSE: -+ flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals) - self.statistics["njacs"] += njevals - else: - flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals) #Number of jac evals -Index: assimulo/lib/sundials_includes.pxd -=================================================================== ---- assimulo/lib/sundials_includes.pxd (revision 844) -+++ assimulo/lib/sundials_includes.pxd (revision 845) -@@ -37,19 +37,6 @@ - ctypedef double realtype
- ctypedef bint booleantype # should be bool instead of bint, but there is a bug in Cython
-
--#==============================================
--# C headers
--#==============================================
--cdef extern from "string.h":
-- void *memcpy(void *s1, void *s2, int n)
--cdef extern from "stdlib.h":
-- void *malloc(int size)
-- void free(void *ptr)
--
--#==============================================
--#External definitions from Sundials headers
--#==============================================
--
- cdef extern from "sundials/sundials_nvector.h":
- ctypedef _generic_N_Vector* N_Vector
-
-@@ -77,6 +64,107 @@ - void N_VDestroy_Serial(N_Vector v)
- void N_VPrint_Serial(N_Vector v)
-
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "sundials/sundials_types.h":
-+ IF SUNDIALS_VECTOR_SIZE == "64":
-+ ctypedef long int sunindextype
-+ ELSE:
-+ ctypedef int sunindextype
-+ cdef extern from "sundials/sundials_matrix.h":
-+ ctypedef _generic_SUNMatrix *SUNMatrix;
-+ void SUNMatDestroy(SUNMatrix A);
-+
-+ cdef struct _generic_SUNMatrix_Ops:
-+ SUNMatrix_ID (*getid)(SUNMatrix);
-+ SUNMatrix (*clone)(SUNMatrix);
-+ void (*destroy)(SUNMatrix);
-+ int (*zero)(SUNMatrix);
-+ int (*copy)(SUNMatrix, SUNMatrix);
-+ int (*scaleadd)(realtype, SUNMatrix, SUNMatrix);
-+ int (*scaleaddi)(realtype, SUNMatrix);
-+ int (*matvec)(SUNMatrix, N_Vector, N_Vector);
-+ int (*space)(SUNMatrix, long int*, long int*);
-+
-+ cdef struct _generic_SUNMatrix:
-+ void *content;
-+ _generic_SUNMatrix_Ops *ops;
-+
-+ cdef enum SUNMatrix_ID:
-+ SUNMATRIX_DENSE,
-+ SUNMATRIX_BAND,
-+ SUNMATRIX_SPARSE,
-+ SUNMATRIX_CUSTOM
-+
-+ cdef extern from "sundials/sundials_linearsolver.h":
-+ ctypedef _generic_SUNLinearSolver *SUNLinearSolver;
-+ int SUNLinSolFree(SUNLinearSolver S);
-+
-+ cdef struct _generic_SUNLinearSolver_Ops:
-+ SUNLinearSolver_Type (*gettype)(SUNLinearSolver);
-+ int (*setatimes)(SUNLinearSolver, void*, ATimesFn);
-+ int (*setpreconditioner)(SUNLinearSolver, void*,
-+ PSetupFn, PSolveFn);
-+ int (*setscalingvectors)(SUNLinearSolver,
-+ N_Vector, N_Vector);
-+ int (*initialize)(SUNLinearSolver);
-+ int (*setup)(SUNLinearSolver, SUNMatrix);
-+ int (*solve)(SUNLinearSolver, SUNMatrix, N_Vector,
-+ N_Vector, realtype);
-+ int (*numiters)(SUNLinearSolver);
-+ realtype (*resnorm)(SUNLinearSolver);
-+ long int (*lastflag)(SUNLinearSolver);
-+ int (*space)(SUNLinearSolver, long int*, long int*);
-+ N_Vector (*resid)(SUNLinearSolver);
-+ int (*free)(SUNLinearSolver);
-+
-+ cdef struct _generic_SUNLinearSolver:
-+ void *content;
-+ _generic_SUNLinearSolver_Ops *ops;
-+
-+ cdef enum SUNLinearSolver_Type:
-+ SUNLINEARSOLVER_DIRECT,
-+ SUNLINEARSOLVER_ITERATIVE,
-+ SUNLINEARSOLVER_CUSTOM
-+
-+ cdef extern from "sunmatrix/sunmatrix_dense.h":
-+ ctypedef _SUNMatrixContent_Dense *SUNMatrixContent_Dense;
-+ cdef struct _SUNMatrixContent_Dense:
-+ sunindextype M;
-+ sunindextype N;
-+ realtype *data;
-+ sunindextype ldata;
-+ realtype **cols;
-+ SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N);
-+ cdef extern from "sunmatrix/sunmatrix_sparse.h":
-+ ctypedef _SUNMatrixContent_Sparse *SUNMatrixContent_Sparse;
-+ cdef struct _SUNMatrixContent_Sparse:
-+ sunindextype M;
-+ sunindextype N;
-+ sunindextype NNZ;
-+ sunindextype NP;
-+ realtype *data;
-+ int sparsetype;
-+ sunindextype *indexvals;
-+ sunindextype *indexptrs;
-+ sunindextype **rowvals;
-+ sunindextype **colptrs;
-+ sunindextype **colvals;
-+ sunindextype **rowptrs;
-+ SUNMatrix SUNSparseMatrix(sunindextype M, sunindextype N, sunindextype NNZ, int sparsetype);
-+ cdef extern from "sunlinsol/sunlinsol_dense.h":
-+ SUNLinearSolver SUNDenseLinearSolver(N_Vector y, SUNMatrix A);
-+ cdef extern from "sunlinsol/sunlinsol_spgmr.h":
-+ SUNLinearSolver SUNSPGMR(N_Vector y, int pretype, int maxl);
-+
-+ELSE:
-+ #Dummy defines
-+ ctypedef void *SUNLinearSolver
-+ ctypedef void *SUNMatrix
-+ ctypedef void *SUNMatrixContent_Dense
-+ ctypedef void *SUNMatrixContent_Sparse
-+ ctypedef int sunindextype
-+
-+
- #Struct for handling the Jacobian data
- cdef extern from "sundials/sundials_direct.h":
- cdef struct _DlsMat:
-@@ -92,7 +180,7 @@ - realtype **cols
- ctypedef _DlsMat* DlsMat
- cdef realtype* DENSE_COL(DlsMat A, int j)
--
-+
- IF SUNDIALS_VERSION >= (2,6,3):
- cdef extern from "sundials/sundials_sparse.h":
- cdef struct _SlsMat:
-@@ -128,7 +216,26 @@ - int *rowvals
- int *colptrs
- ctypedef _SlsMat* SlsMat
-+
-+#==============================================
-+# C headers
-+#==============================================
-+cdef extern from "string.h":
-+ void *memcpy(void *s1, void *s2, int n)
-+cdef extern from "stdlib.h":
-+ void *malloc(int size)
-+ void free(void *ptr)
-+
-+#==============================================
-+#External definitions from Sundials headers
-+#==============================================
-
-+IF SUNDIALS_WITH_SUPERLU:
-+ cdef inline int with_superlu(): return 1
-+ELSE:
-+ cdef inline int with_superlu(): return 0
-+
-+
- cdef extern from "cvodes/cvodes.h":
- void* CVodeCreate(int lmm, int iter)
- ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void *f_data)
-@@ -225,45 +332,67 @@ - int CVodeGetSensNonlinSolvStats(void *cvode_mem, long int *nSniters, long int *nSncfails)
- int CVodeGetStgrSensNumNonlinSolvIters(void *cvode_mem, long int *nSTGR1niters)
- int CVodeGetStgrSensNumNonlinSolvConvFails(void *cvode_mem, long int *nSTGR1ncfails)
-+
-+cdef extern from "cvodes/cvodes_spils.h":
-+ ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
-+ N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "cvodes/cvodes_direct.h":
-+ ctypedef int (*CVDlsDenseJacFn)(realtype t, N_Vector y, N_Vector fy,
-+ SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int CVDlsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS, SUNMatrix A);
-+ int CVDlsSetJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
-+ cdef extern from "cvodes/cvodes_spils.h":
-+ int CVSpilsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS);
-+ ctypedef int (*CVSpilsJacTimesSetupFn)(realtype t, N_Vector y, N_Vector fy, void *user_data);
-+ int CVSpilsSetJacTimes(void *cvode_mem, CVSpilsJacTimesSetupFn jtsetup, CVSpilsJacTimesVecFn jtimes);
-
-
--cdef extern from "cvodes/cvodes_dense.h":
-- int CVDense(void *cvode_mem, long int n)
-- ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy,
-- DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-- int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
-+ IF SUNDIALS_WITH_SUPERLU:
-+ cdef extern from "sunlinsol/sunlinsol_superlumt.h":
-+ SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads)
-+ ELSE:
-+ cdef inline SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads): return NULL
-+
-+ cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0
-+ cdef inline tuple version(): return (3,0,0)
-+ELSE:
-+ cdef extern from "cvodes/cvodes_dense.h":
-+ int CVDense(void *cvode_mem, long int n)
-+ ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy,
-+ DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
-
--cdef extern from "cvodes/cvodes_spgmr.h":
-- int CVSpgmr(void *cvode_mem, int pretype, int max1)
-+ cdef extern from "cvodes/cvodes_spgmr.h":
-+ int CVSpgmr(void *cvode_mem, int pretype, int max1)
-
--IF SUNDIALS_VERSION >= (2,6,0):
-- cdef extern from "cvodes/cvodes_sparse.h":
-+ cdef extern from "cvodes/cvodes_spils.h":
-+ int CVSpilsSetJacTimesVecFn(void *cvode_mem, CVSpilsJacTimesVecFn jtv)
-+
-+ IF SUNDIALS_VERSION >= (2,6,0):
-+ cdef extern from "cvodes/cvodes_sparse.h":
-+ ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
-+ SlsMat Jac, void *user_data, N_Vector tmp1,
-+ N_Vector tmp2, N_Vector tmp3)
-+ int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
-+ int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
-+ cdef inline tuple version(): return (2,6,0)
-+ IF SUNDIALS_WITH_SUPERLU:
-+ cdef extern from "cvodes/cvodes_sparse.h":
-+ int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
-+ ELSE:
-+ cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
-+ ELSE:
-+ cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
- ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
- SlsMat Jac, void *user_data, N_Vector tmp1,
- N_Vector tmp2, N_Vector tmp3)
-- int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
-- int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
-- #cdef inline char* version(): return "2.6.0"
-- cdef inline tuple version(): return (2,6,0)
-- IF SUNDIALS_WITH_SUPERLU:
-- cdef extern from "cvodes/cvodes_sparse.h":
-- int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
-- ELSE:
-- cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
--ELSE:
-- cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
-- ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
-- SlsMat Jac, void *user_data, N_Vector tmp1,
-- N_Vector tmp2, N_Vector tmp3)
-- cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
-- cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
-- #cdef inline char* version(): return "2.5.0"
-- cdef inline tuple version(): return (2,5,0)
-+ cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
-+ cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
-+ cdef inline tuple version(): return (2,5,0)
-
- cdef extern from "cvodes/cvodes_spils.h":
-- ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
-- N_Vector y, N_Vector fy,
-- void *user_data, N_Vector tmp)
- ctypedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy,
- booleantype jok, booleantype *jcurPtr,
- realtype gamma, void *user_data,
-@@ -273,13 +402,12 @@ - N_Vector r, N_Vector z,
- realtype gamma, realtype delta,
- int lr, void *user_data, N_Vector tmp)
-- int CVSpilsSetJacTimesVecFn(void *cvode_mem, CVSpilsJacTimesVecFn jtv)
-+
- int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
- int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
- int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
- int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals)
- int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves)
--
-
- cdef extern from "idas/idas.h":
- ctypedef int (*IDAResFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
-@@ -380,19 +508,40 @@ -
- #End Sensitivities
- #=================
--
--cdef extern from "idas/idas_dense.h":
-- int IDADense(void *ida_mem, long int n)
-- ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy,
-- N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data,
-- N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-- int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
-
--cdef extern from "idas/idas_dense.h":
-- int IDASpgmr(void *ida_mem, int max1)
-- ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
-- int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
-+cdef extern from "idas/idas_spils.h":
-+ ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr,
-+ N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
-
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "idas/idas_direct.h":
-+ ctypedef int (*IDADlsDenseJacFn)(realtype tt, realtype cj, N_Vector yy,
-+ N_Vector yp, N_Vector rr, SUNMatrix Jac, void *user_data,
-+ N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int IDADlsSetJacFn(void *ida_mem, IDADlsDenseJacFn djac)
-+ int IDADlsSetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix A);
-+
-+ cdef extern from "idas/idas_spils.h":
-+ int IDASpilsSetLinearSolver(void *ida_mem, SUNLinearSolver LS);
-+ ctypedef int (*IDASpilsJacTimesSetupFn)(realtype tt, N_Vector yy,
-+ N_Vector yp, N_Vector rr, realtype c_j, void *user_data);
-+ int IDASpilsSetJacTimes(void *ida_mem,
-+ IDASpilsJacTimesSetupFn jtsetup, IDASpilsJacTimesVecFn jtimes);
-+
-+ cdef inline int ida_spils_jtsetup_dummy(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, realtype c_j, void *user_data): return 0
-+ELSE:
-+ cdef extern from "idas/idas_dense.h":
-+ int IDADense(void *ida_mem, long int n)
-+ ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy,
-+ N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data,
-+ N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
-+
-+ int IDASpgmr(void *ida_mem, int max1)
-+
-+ cdef extern from "idas/idas_spils.h":
-+ int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
-+
- cdef extern from "idas/idas_spils.h":
- int IDASpilsGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
- int IDASpilsGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
-@@ -455,14 +604,31 @@ - # fuction used to deallocate memory used by KINSOL
- void KINFree(void **kinmem)
-
--# functions used for supplying jacobian, and receiving info from linear solver
--cdef extern from "kinsol/kinsol_direct.h":
-- # user functions
-- ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "kinsol/kinsol_direct.h":
-+ ctypedef int (*KINDlsDenseJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
-+ int KINDlsSetLinearSolver(void *kinmem, SUNLinearSolver LS, SUNMatrix A);
-+ int KINDlsSetJacFn(void *kinmem, KINDlsDenseJacFn djac)
-
-- # function used to link user functions to KINSOL
-- int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
-+ cdef extern from "kinsol/kinsol_spils.h":
-+ int KINSpilsSetLinearSolver(void *kinsol_mem, SUNLinearSolver LS);
-+ELSE:
-+ # functions used for supplying jacobian, and receiving info from linear solver
-+ cdef extern from "kinsol/kinsol_direct.h":
-+ # user functions
-+ ctypedef int (*KINDlsDenseJacFn)(int dim, 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)
-+
-+ cdef extern from "kinsol/kinsol_dense.h":
-+ int KINDense(void *kinmem, int dim)
-+
-+ cdef extern from "kinsol/kinsol_spgmr.h":
-+ int KINSpgmr(void *kinmem, int maxl)
-
-+cdef extern from "kinsol/kinsol_direct.h":
- # optional output fcts for linear direct solver
- int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB)
- int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB)
-@@ -470,12 +636,6 @@ - int KINDlsGetLastFlag(void *kinmem, int *flag)
- char *KINDlsGetReturnFlagName(int flag)
-
--cdef extern from "kinsol/kinsol_dense.h":
-- int KINDense(void *kinmem, int dim)
--
--cdef extern from "kinsol/kinsol_spgmr.h":
-- int KINSpgmr(void *kinmem, int maxl)
--
- cdef extern from "kinsol/kinsol_spils.h":
- ctypedef int (*KINSpilsJacTimesVecFn)(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
- void *problem_data)
-Index: src/lib/sundials_constants.pxi -=================================================================== ---- assimulo/lib/sundials_constants.pxi (revision 844) -+++ assimulo/lib/sundials_constants.pxi (revision 845) -@@ -95,6 +95,9 @@ - DEF CV_BAD_IS = -45
-
-
-+DEF CSC_MAT = 0
-+DEF CSR_MAT = 1
-+
- #==========
- # IDA
- #==========
-Index: src/lib/sundials_callbacks.pxi -=================================================================== ---- assimulo/lib/sundials_callbacks.pxi (revision 844) -+++ assimulo/lib/sundials_callbacks.pxi (revision 845) -@@ -94,98 +94,55 @@ - traceback.print_exc()
- return CV_UNREC_RHSFUNC_ERR
-
--@cython.boundscheck(False)
--@cython.wraparound(False)
--cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
-- void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-- """
-- This method is used to connect the Assimulo.Problem.jac to the Sundials
-- Sparse Jacobian function.
-- """
-- cdef ProblemData pData = <ProblemData>problem_data
-- #cdef N.ndarray y = nv2arr(yv)
-- cdef N.ndarray y = pData.work_y
-- cdef int i
-- cdef int nnz = Jacobian.NNZ
-- cdef int ret_nnz
-- cdef int dim = Jacobian.N
-- cdef realtype* data = Jacobian.data
--
-- IF SUNDIALS_VERSION >= (2,6,3):
-- cdef int* rowvals = Jacobian.rowvals[0]
-- cdef int* colptrs = Jacobian.colptrs[0]
-- ELSE:
-- cdef int* rowvals = Jacobian.rowvals
-- cdef int* colptrs = Jacobian.colptrs
--
-- nv2arr_inplace(yv, y)
-- """
-- realtype *data;
-- int *rowvals;
-- int *colptrs;
-- """
-- try:
-- if pData.dimSens > 0: #Sensitivity activated
-- p = realtype2arr(pData.p,pData.dimSens)
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-- else:
-- jac=(<object>pData.JAC)(t,y,p=p)
-- else:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-- else:
-- jac=(<object>pData.JAC)(t,y)
--
-- if not isinstance(jac, sparse.csc.csc_matrix):
-- jac = sparse.csc.csc_matrix(jac)
-- raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-- ret_nnz = jac.nnz
-- if ret_nnz > nnz:
-- raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")
-
-- for i in range(min(ret_nnz,nnz)):
-- data[i] = jac.data[i]
-- rowvals[i] = jac.indices[i]
-- for i in range(dim+1):
-- colptrs[i] = jac.indptr[i]
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ @cython.boundscheck(False)
-+ @cython.wraparound(False)
-+ cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Sparse Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef SUNMatrixContent_Sparse Jacobian = <SUNMatrixContent_Sparse>Jac.content
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i
-+ cdef sunindextype nnz = Jacobian.NNZ
-+ cdef int ret_nnz
-+ cdef sunindextype dim = Jacobian.N
-+ cdef realtype* data = Jacobian.data
-+ cdef sunindextype* rowvals = Jacobian.rowvals[0]
-+ cdef sunindextype* colptrs = Jacobian.colptrs[0]
-
-- return CVDLS_SUCCESS
-- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-- return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-- except:
-- traceback.print_exc()
-- return CVDLS_JACFUNC_UNRECVR
-+ nv2arr_inplace(yv, y)
-
--cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian,
-- void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-- """
-- This method is used to connect the Assimulo.Problem.jac to the Sundials
-- Jacobian function.
-- """
-- cdef ProblemData pData = <ProblemData>problem_data
-- #cdef ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-- cdef realtype* col_i=DENSE_COL(Jacobian,0)
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #cdef N.ndarray y = nv2arr(yv)
-- cdef N.ndarray y = pData.work_y
-- cdef int i,j
--
-- nv2arr_inplace(yv, y)
--
-- if pData.dimSens>0: #Sensitivity activated
-- p = realtype2arr(pData.p,pData.dimSens)
- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
-+ if pData.dimSens > 0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p=p)
- else:
-- jac=(<object>pData.JAC)(t,y,p)
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-+ if not isinstance(jac, sparse.csc.csc_matrix):
-+ jac = sparse.csc.csc_matrix(jac)
-+ raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-+ ret_nnz = jac.nnz
-+ if ret_nnz > nnz:
-+ raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")
-
-+ for i in range(min(ret_nnz,nnz)):
-+ data[i] = jac.data[i]
-+ rowvals[i] = jac.indices[i]
-+ for i in range(dim+1):
-+ colptrs[i] = jac.indptr[i]
-+
- return CVDLS_SUCCESS
- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
- return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-@@ -192,18 +149,62 @@ - except:
- traceback.print_exc()
- return CVDLS_JACFUNC_UNRECVR
-- else:
-+ELSE:
-+ @cython.boundscheck(False)
-+ @cython.wraparound(False)
-+ cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Sparse Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i
-+ cdef int nnz = Jacobian.NNZ
-+ cdef int ret_nnz
-+ cdef int dim = Jacobian.N
-+ cdef realtype* data = Jacobian.data
-+
-+ IF SUNDIALS_VERSION >= (2,6,3):
-+ cdef int* rowvals = Jacobian.rowvals[0]
-+ cdef int* colptrs = Jacobian.colptrs[0]
-+ ELSE:
-+ cdef int* rowvals = Jacobian.rowvals
-+ cdef int* colptrs = Jacobian.colptrs
-+
-+ nv2arr_inplace(yv, y)
-+ """
-+ realtype *data;
-+ int *rowvals;
-+ int *colptrs;
-+ """
- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ if pData.dimSens > 0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p=p)
- else:
-- jac=(<object>pData.JAC)(t,y)
--
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-+
-+ if not isinstance(jac, sparse.csc.csc_matrix):
-+ jac = sparse.csc.csc_matrix(jac)
-+ raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-+ ret_nnz = jac.nnz
-+ if ret_nnz > nnz:
-+ raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")
-
-+ for i in range(min(ret_nnz,nnz)):
-+ data[i] = jac.data[i]
-+ rowvals[i] = jac.indices[i]
-+ for i in range(dim+1):
-+ colptrs[i] = jac.indptr[i]
-+
- return CVDLS_SUCCESS
- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
- return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-@@ -210,7 +211,112 @@ - except:
- traceback.print_exc()
- return CVDLS_JACFUNC_UNRECVR
-+
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef int cv_jac(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef realtype* col_i=Jacobian.cols[0]
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i,j, Neq = pData.dim
-
-+ nv2arr_inplace(yv, y)
-+
-+ if pData.dimSens>0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+ELSE:
-+ cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef realtype* col_i=DENSE_COL(Jacobian,0)
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i,j
-+
-+ nv2arr_inplace(yv, y)
-+
-+ if pData.dimSens>0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+
- cdef int cv_jacv(N_Vector vv, N_Vector Jv, realtype t, N_Vector yv, N_Vector fyv,
- void *problem_data, N_Vector tmp):
- """
-@@ -337,9 +443,6 @@ - Root-finding function.
- """
- cdef ProblemData pData = <ProblemData>problem_data
-- #cdef ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #cdef N.ndarray y = nv2arr(yv)
- cdef N.ndarray y = pData.work_y
- cdef int i
-
-@@ -366,12 +469,8 @@ - """
- cdef ProblemData pData = <ProblemData>problem_data
- cdef N.ndarray[realtype, ndim=1, mode='c'] res #Used for return from the user function
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #(<ndarray>pData.yd).data = <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
- cdef N.ndarray y = pData.work_y
- cdef N.ndarray yd = pData.work_yd
-- # cdef N.ndarray y = nv2arr(yv)
-- # cdef N.ndarray yd = nv2arr(yvdot)
- cdef realtype* resptr=(<N_VectorContent_Serial>residual.content).data
- cdef int i
-
-@@ -414,63 +513,113 @@ - except:
- traceback.print_exc()
- return IDA_RES_FAIL
--
--cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef int ida_jac(realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-+ cdef realtype* col_i=Jacobian.cols[0]
-+ cdef N.ndarray y = pData.work_y
-+ cdef N.ndarray yd = pData.work_yd
-+ cdef int i,j, Neq = pData.dim
-+
-+ nv2arr_inplace(yv, y)
-+ nv2arr_inplace(yvdot, yd)
-+
-+ if pData.dimSens!=0: #SENSITIVITY
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd,p=p)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-+ELSE:
-+ cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
- void* problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-- """
-- This method is used to connect the Assimulo.Problem.jac to the Sundials
-- Jacobian function.
-- """
-- cdef ProblemData pData = <ProblemData>problem_data
-- cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-- cdef realtype* col_i=DENSE_COL(Jacobian,0)
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #(<ndarray>pData.yd).data = <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
-- cdef N.ndarray y = pData.work_y
-- cdef N.ndarray yd = pData.work_yd
-- #cdef N.ndarray y = nv2arr(yv)
-- #cdef N.ndarray yd = nv2arr(yvdot)
-- cdef int i,j
--
-- nv2arr_inplace(yv, y)
-- nv2arr_inplace(yvdot, yd)
--
-- if pData.dimSens!=0: #SENSITIVITY
-- p = realtype2arr(pData.p,pData.dimSens)
-- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p) # call to the python residual function
-- else:
-- jac=(<object>pData.JAC)(c,t,y,yd,p=p)
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-+ cdef realtype* col_i=DENSE_COL(Jacobian,0)
-+ cdef N.ndarray y = pData.work_y
-+ cdef N.ndarray yd = pData.work_yd
-+ cdef int i,j
-+
-+ nv2arr_inplace(yv, y)
-+ nv2arr_inplace(yvdot, yd)
-+
-+ if pData.dimSens!=0: #SENSITIVITY
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd,p=p)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-- return IDADLS_SUCCESS
-- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-- return IDADLS_JACFUNC_RECVR #Recoverable Error
-- except:
-- traceback.print_exc()
-- return IDADLS_JACFUNC_UNRECVR
-- else:
-- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw) # call to the python residual function
-- else:
-- jac=(<object>pData.JAC)(c,t,y,yd)
--
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-- return IDADLS_SUCCESS
-- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-- return IDADLS_JACFUNC_RECVR #Recoverable Error
-- except:
-- traceback.print_exc()
-- return IDADLS_JACFUNC_UNRECVR
--
-
- cdef int ida_root(realtype t, N_Vector yv, N_Vector yvdot, realtype *gout, void* problem_data):
- """
-@@ -479,12 +628,8 @@ - """
- cdef ProblemData pData = <ProblemData>problem_data
- cdef N.ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #(<ndarray>pData.yd).data = <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
- cdef N.ndarray y = pData.work_y
- cdef N.ndarray yd = pData.work_yd
-- #cdef N.ndarray y = nv2arr(yv)
-- #cdef N.ndarray yd = nv2arr(yvdot)
- cdef int i
-
- nv2arr_inplace(yv, y)
-@@ -553,30 +698,54 @@ - traceback.print_exc()
- return SPGMR_PSOLVE_FAIL_UNREC
-
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef int kin_jac(N_Vector xv, N_Vector fval, SUNMatrix Jac,
-+ 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 SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
-+ cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-+ cdef realtype* col_i=Jacobian.cols[0]
-+ cdef N.ndarray x = nv2arr(xv)
-+ cdef int i,j, Neq = pData.dim
-+
-+ try:
-+ jac=(<object>pData.JAC)(x)
-
--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 ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-- cdef realtype* col_i=DENSE_COL(Jacobian,0)
-- cdef N.ndarray x = nv2arr(xv)
-- cdef int i,j
--
-- try:
-- jac=(<object>pData.JAC)(x)
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-
-- 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)
-+ELSE:
-+ 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 ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-+ cdef realtype* col_i=DENSE_COL(Jacobian,0)
-+ cdef N.ndarray x = nv2arr(xv)
-+ cdef int i,j
-+
-+ try:
-+ jac=(<object>pData.JAC)(x)
-
-- return KINDLS_SUCCESS
-- except:
-- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
--
-+ 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_jacv(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
- void *problem_data):
- cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-Index: examples/cvode_with_jac.py -=================================================================== ---- assimulo/examples/cvode_with_jac.py (revision 844) -+++ assimulo/examples/cvode_with_jac.py (revision 845) -@@ -46,7 +46,6 @@ - def f(t,y): - yd_0 = y[1] - yd_1 = -9.82 -- #print y, yd_0, yd_1 - return N.array([yd_0,yd_1]) - - #Defines the Jacobian -Index: examples/__init__.py -=================================================================== ---- assimulo/examples/__init__.py (revision 844) -+++ assimulo/examples/__init__.py (revision 845) -@@ -24,6 +24,6 @@ - "mech_system_pendulum", "euler_vanderpol", "cvode_with_parameters_modified", - "cvode_basic_backward","ida_basic_backward","dasp3_basic", "cvode_with_preconditioning", - "kinsol_basic","kinsol_with_jac", "radau5dae_time_events", "kinsol_ors", "lsodar_bouncing_ball", -- "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result"] -+ "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result", "cvode_with_jac_sparse"] - - -Index: examples/cvode_with_jac_sparse.py -=================================================================== ---- assimulo/examples/cvode_with_jac_sparse.py (nonexistent) -+++ assimulo/examples/cvode_with_jac_sparse.py (revision 845) -@@ -0,0 +1,102 @@ -+#!/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/>. -+ -+import numpy as N -+import scipy.sparse as SP -+import pylab as P -+import nose -+from assimulo.solvers import CVode -+from assimulo.problem import Explicit_Problem -+ -+ -+def run_example(with_plots=True): -+ r""" -+ Example for demonstrating the use of a user supplied Jacobian (sparse). -+ Note that this will only work if Assimulo has been configured with -+ Sundials + SuperLU. Based on the SUNDIALS example cvRoberts_sps.c -+ -+ ODE: -+ -+ .. math:: -+ -+ \dot y_1 &= -0.04y_1 + 1e4 y_2 y_3 \\ -+ \dot y_2 &= - \dot y_1 - \dot y_3 \\ -+ \dot y_3 &= 3e7 y_2^2 -+ -+ -+ on return: -+ -+ - :dfn:`exp_mod` problem instance -+ -+ - :dfn:`exp_sim` solver instance -+ -+ """ -+ -+ #Defines the rhs -+ def f(t,y): -+ yd_0 = -0.04*y[0] + 1e4*y[1]*y[2] -+ yd_2 = 3e7*y[1]*y[1] -+ yd_1 = -yd_0 - yd_2 -+ return N.array([yd_0,yd_1,yd_2]) -+ -+ #Defines the Jacobian -+ def jac(t,y): -+ -+ colptrs = [0,3,6,9] -+ rowvals = [0, 1, 2, 0, 1, 2, 0, 1, 2] -+ data = [-0.04, 0.04, 0.0, 1e4*y[2], -1e4*y[2]-6e7*y[1], 6e7*y[1], 1e4*y[1], -1e4*y[1], 0.0] -+ -+ J = SP.csc_matrix((data, rowvals, colptrs)) -+ return J -+ -+ #Defines an Assimulo explicit problem -+ y0 = [1.0,0.0,0.0] #Initial conditions -+ -+ exp_mod = Explicit_Problem(f,y0, name = 'Example using analytic (sparse) Jacobian') -+ -+ exp_mod.jac = jac #Sets the Jacobian -+ exp_mod.jac_nnz = 9 -+ -+ -+ exp_sim = CVode(exp_mod) #Create a CVode solver -+ -+ #Set the parameters -+ exp_sim.iter = 'Newton' #Default 'FixedPoint' -+ exp_sim.discr = 'BDF' #Default 'Adams' -+ exp_sim.atol = [1e-8,1e-14,1e-6] #Default 1e-6 -+ exp_sim.rtol = 1e-4 #Default 1e-6 -+ exp_sim.linear_solver = "sparse" -+ -+ #Simulate -+ t, y = exp_sim.simulate(0.4) #Simulate 0.4 seconds -+ -+ #Basic tests -+ nose.tools.assert_almost_equal(y[-1][0],0.9851,3) -+ -+ #Plot -+ if with_plots: -+ P.plot(t,y[:,1],linestyle="dashed",marker="o") #Plot the solution -+ P.xlabel('Time') -+ P.ylabel('State') -+ P.title(exp_mod.name) -+ P.show() -+ -+ return exp_mod, exp_sim -+ -+ -+if __name__=='__main__': -+ mod,sim = run_example() -Index: tests/test_examples.py -=================================================================== ---- assimulo/tests/test_examples.py (revision 844) -+++ assimulo/tests/test_examples.py (revision 845) -@@ -22,7 +22,15 @@ -
- class Test_Examples:
-
-+
- @testattr(stddist = True)
-+ def test_cvode_with_jac_sparse(self):
-+ try:
-+ cvode_with_jac_sparse.run_example(with_plots=False)
-+ except AssimuloException:
-+ pass #Handle the case when SuperLU is not installed
-+
-+ @testattr(stddist = True)
- def test_ida_with_user_defined_handle_result(self):
- ida_with_user_defined_handle_result.run_example(with_plots=False)
-
|