summarylogtreecommitdiffstats
path: root/r845.patch
diff options
context:
space:
mode:
Diffstat (limited to 'r845.patch')
-rw-r--r--r845.patch1689
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)
-