diff options
author | Michel Zou | 2018-05-16 15:15:21 +0200 |
---|---|---|
committer | Michel Zou | 2018-05-16 15:15:21 +0200 |
commit | eda92d5048e66f469b2ec3a30de46c3196cc0be0 (patch) | |
tree | e220568aff2a0e96ce6d7164cc4df7bf093d0003 | |
parent | 791c0310530f3096903ece1025fbdfe454fc9c17 (diff) | |
download | aur-eda92d5048e66f469b2ec3a30de46c3196cc0be0.tar.gz |
fix lf
-rw-r--r-- | .SRCINFO | 2 | ||||
-rw-r--r-- | r831.patch | 118 | ||||
-rw-r--r-- | r833.patch | 148 | ||||
-rw-r--r-- | r836.patch | 164 | ||||
-rw-r--r-- | r837.patch | 274 | ||||
-rw-r--r-- | r838.patch | 166 | ||||
-rw-r--r-- | r840.patch | 114 | ||||
-rw-r--r-- | r845.patch | 2276 |
8 files changed, 1631 insertions, 1631 deletions
@@ -1,5 +1,5 @@ # Generated by mksrcinfo v8 -# Wed May 16 13:13:51 UTC 2018 +# Wed May 16 13:15:21 UTC 2018 pkgbase = python-assimulo pkgdesc = A package for solving ordinary differential equations and differential algebraic equations pkgver = 2.9 diff --git a/r831.patch b/r831.patch index 82b7f7ef5d13..1215207005ad 100644 --- a/r831.patch +++ b/r831.patch @@ -3,62 +3,62 @@ Index: src/solvers/__init__.py --- assimulo/solvers/__init__.py (revision 830) +++ assimulo/solvers/__init__.py (revision 831) @@ -19,16 +19,50 @@ - "glimda","odepack","radar5","dasp3","odassl"] - - #Import all the solvers from the different modules --from .euler import ExplicitEuler, ImplicitEuler --from .radau5 import Radau5ODE, Radau5DAE, _Radau5ODE, _Radau5DAE --from .sundials import IDA, CVode --from .kinsol import KINSOL --from .runge_kutta import RungeKutta34, RungeKutta4, Dopri5 --from .rosenbrock import RodasODE --from .odassl import ODASSL --from .odepack import LSODAR --from .radar5 import Radar5ODE - try: -+ from .euler import ExplicitEuler -+ from .euler import ImplicitEuler -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .radau5 import Radau5ODE -+ from .radau5 import Radau5DAE -+ from .radau5 import _Radau5ODE -+ from .radau5 import _Radau5DAE -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .sundials import IDA -+ from .sundials import CVode -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .kinsol import KINSOL -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .runge_kutta import RungeKutta34 -+ from .runge_kutta import RungeKutta4 -+ from .runge_kutta import Dopri5 -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .rosenbrock import RodasODE -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .odassl import ODASSL -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .odepack import LSODAR -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: -+ from .radar5 import Radar5ODE -+except ImportError as ie: -+ print("Could not find {}".format(ie.args[0][16:])) -+try: - from .dasp3 import DASP3ODE - except ImportError as ie: - print("Could not find {}".format(ie.args[0][16:])) + "glimda","odepack","radar5","dasp3","odassl"]
+
+ #Import all the solvers from the different modules
+-from .euler import ExplicitEuler, ImplicitEuler
+-from .radau5 import Radau5ODE, Radau5DAE, _Radau5ODE, _Radau5DAE
+-from .sundials import IDA, CVode
+-from .kinsol import KINSOL
+-from .runge_kutta import RungeKutta34, RungeKutta4, Dopri5
+-from .rosenbrock import RodasODE
+-from .odassl import ODASSL
+-from .odepack import LSODAR
+-from .radar5 import Radar5ODE
+ try:
++ from .euler import ExplicitEuler
++ from .euler import ImplicitEuler
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .radau5 import Radau5ODE
++ from .radau5 import Radau5DAE
++ from .radau5 import _Radau5ODE
++ from .radau5 import _Radau5DAE
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .sundials import IDA
++ from .sundials import CVode
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .kinsol import KINSOL
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .runge_kutta import RungeKutta34
++ from .runge_kutta import RungeKutta4
++ from .runge_kutta import Dopri5
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .rosenbrock import RodasODE
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .odassl import ODASSL
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .odepack import LSODAR
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
++ from .radar5 import Radar5ODE
++except ImportError as ie:
++ print("Could not find {}".format(ie.args[0][16:]))
++try:
+ from .dasp3 import DASP3ODE
+ except ImportError as ie:
+ print("Could not find {}".format(ie.args[0][16:]))
diff --git a/r833.patch b/r833.patch index 090d2adcc4dc..7385443600c1 100644 --- a/r833.patch +++ b/r833.patch @@ -3,83 +3,83 @@ Index: src/solvers/radar5.py --- assimulo/solvers/radar5.py (revision 832) +++ assimulo/solvers/radar5.py (revision 833) @@ -26,9 +26,11 @@ - from assimulo.explicit_ode import Explicit_ODE - from assimulo.implicit_ode import Implicit_ODE - --from assimulo.lib import radar5 -+try: -+ from assimulo.lib import radar5 -+except ImportError: -+ print("Could not find RADAR5") - -- - class Radar_Exception(Exception): - pass - + from assimulo.explicit_ode import Explicit_ODE
+ from assimulo.implicit_ode import Implicit_ODE
+
+-from assimulo.lib import radar5
++try:
++ from assimulo.lib import radar5
++except ImportError:
++ print("Could not find RADAR5")
+
+-
+ class Radar_Exception(Exception):
+ pass
+
Index: src/solvers/__init__.py =================================================================== --- assimulo/solvers/__init__.py (revision 832) +++ assimulo/solvers/__init__.py (revision 833) @@ -23,7 +23,7 @@ - from .euler import ExplicitEuler - from .euler import ImplicitEuler - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .radau5 import Radau5ODE - from .radau5 import Radau5DAE + from .euler import ExplicitEuler
+ from .euler import ImplicitEuler
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .radau5 import Radau5ODE
+ from .radau5 import Radau5DAE
@@ -30,43 +30,43 @@ - from .radau5 import _Radau5ODE - from .radau5 import _Radau5DAE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .sundials import IDA - from .sundials import CVode - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .kinsol import KINSOL - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .runge_kutta import RungeKutta34 - from .runge_kutta import RungeKutta4 - from .runge_kutta import Dopri5 - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .rosenbrock import RodasODE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .odassl import ODASSL - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .odepack import LSODAR - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .radar5 import Radar5ODE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .dasp3 import DASP3ODE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) - try: - from .glimda import GLIMDA - except ImportError as ie: -- print("Could not find {}".format(ie.args[0][16:])) -+ print("Could not find {}".format(ie.args[0].split("'")[1])) + from .radau5 import _Radau5ODE
+ from .radau5 import _Radau5DAE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .sundials import IDA
+ from .sundials import CVode
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .kinsol import KINSOL
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .runge_kutta import RungeKutta34
+ from .runge_kutta import RungeKutta4
+ from .runge_kutta import Dopri5
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .rosenbrock import RodasODE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .odassl import ODASSL
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .odepack import LSODAR
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .radar5 import Radar5ODE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .dasp3 import DASP3ODE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
+ try:
+ from .glimda import GLIMDA
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0][16:]))
++ print("Could not find {}".format(ie.args[0].split("'")[1]))
diff --git a/r836.patch b/r836.patch index 50012932459b..b3cc0a96ba5d 100644 --- a/r836.patch +++ b/r836.patch @@ -3,90 +3,90 @@ Index: src/lib/sundials_callbacks.pxi --- assimulo/lib/sundials_callbacks.pxi (revision 835) +++ assimulo/lib/sundials_callbacks.pxi (revision 836) @@ -103,7 +103,8 @@ - Sparse Jacobian function. - """ - cdef ProblemData pData = <ProblemData>problem_data -- cdef N.ndarray y = nv2arr(yv) -+ #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 + Sparse Jacobian function.
+ """
+ cdef ProblemData pData = <ProblemData>problem_data
+- cdef N.ndarray y = nv2arr(yv)
++ #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
@@ -111,7 +112,8 @@ - cdef realtype* data = Jacobian.data - cdef int* rowvals = Jacobian.rowvals - cdef int* colptrs = Jacobian.colptrs -- -+ -+ nv2arr_inplace(yv, y) - """ - realtype *data; - int *rowvals; + cdef realtype* data = Jacobian.data
+ cdef int* rowvals = Jacobian.rowvals
+ cdef int* colptrs = Jacobian.colptrs
+-
++
++ nv2arr_inplace(yv, y)
+ """
+ realtype *data;
+ int *rowvals;
@@ -160,8 +162,11 @@ - #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 = 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) + #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 = 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)
@@ -358,11 +363,16 @@ - 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 = nv2arr(yv) -- cdef N.ndarray yd = nv2arr(yvdot) -+ 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 - -+ nv2arr_inplace(yv, y) -+ nv2arr_inplace(yvdot, yd) -+ - if pData.dimSens!=0: #SENSITIVITY - p = realtype2arr(pData.p,pData.dimSens) - try: + 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 = nv2arr(yv)
+- cdef N.ndarray yd = nv2arr(yvdot)
++ 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
+
++ nv2arr_inplace(yv, y)
++ nv2arr_inplace(yvdot, yd)
++
+ if pData.dimSens!=0: #SENSITIVITY
+ p = realtype2arr(pData.p,pData.dimSens)
+ try:
@@ -411,10 +421,15 @@ - 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 = nv2arr(yv) -- cdef N.ndarray yd = nv2arr(yvdot) -+ 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: + 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 = nv2arr(yv)
+- cdef N.ndarray yd = nv2arr(yvdot)
++ 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:
@@ -461,10 +476,15 @@ - 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 = nv2arr(yv) -- cdef N.ndarray yd = nv2arr(yvdot) -+ 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) -+ nv2arr_inplace(yvdot, yd) -+ - try: - if pData.sw != NULL: - root=(<object>pData.ROOT)(t,y,yd,<list>pData.sw) #Call to the Python root function + 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 = nv2arr(yv)
+- cdef N.ndarray yd = nv2arr(yvdot)
++ 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)
++ nv2arr_inplace(yvdot, yd)
++
+ try:
+ if pData.sw != NULL:
+ root=(<object>pData.ROOT)(t,y,yd,<list>pData.sw) #Call to the Python root function
diff --git a/r837.patch b/r837.patch index 6123aaacdee6..e5fdd761d3a3 100644 --- a/r837.patch +++ b/r837.patch @@ -3,82 +3,82 @@ Index: src/implicit_ode.pyx --- assimulo/implicit_ode.pyx (revision 836) +++ assimulo/implicit_ode.pyx (revision 837) @@ -26,7 +26,7 @@ - cimport numpy as N - - from exception import * --from time import clock, time -+from timeit import default_timer as timer - import warnings - - realtype = N.float + cimport numpy as N
+
+ from exception import *
+-from time import clock, time
++from timeit import default_timer as timer
+ import warnings
+
+ realtype = N.float
@@ -187,7 +187,7 @@ - output_index = 0 - - self.time_limit_activated = 1 if self.time_limit > 0 else 0 -- self.time_integration_start = time() -+ self.time_integration_start = timer() - - while (flag == ID_COMPLETE and tevent == tfinal) is False and (self.t-eps > tfinal) if backward else (self.t+eps < tfinal): - + output_index = 0
+
+ self.time_limit_activated = 1 if self.time_limit > 0 else 0
+- self.time_integration_start = time()
++ self.time_integration_start = timer()
+
+ while (flag == ID_COMPLETE and tevent == tfinal) is False and (self.t-eps > tfinal) if backward else (self.t+eps < tfinal):
+
@@ -203,7 +203,7 @@ - - #Initialize the clock, enabling storing elapsed time for each step - if REPORT_CONTINUOUSLY and self.options["clock_step"]: -- self.clock_start = clock() -+ self.clock_start = timer() - - [flag, tlist, ylist, ydlist] = self.integrate(self.t, self.y, self.yd, tevent, opts) - +
+ #Initialize the clock, enabling storing elapsed time for each step
+ if REPORT_CONTINUOUSLY and self.options["clock_step"]:
+- self.clock_start = clock()
++ self.clock_start = timer()
+
+ [flag, tlist, ylist, ydlist] = self.integrate(self.t, self.y, self.yd, tevent, opts)
+
@@ -274,16 +274,16 @@ - - #Store the elapsed time for a single step - if self.options["clock_step"]: -- self.elapsed_step_time = clock() - self.clock_start -- self.clock_start = clock() -+ self.elapsed_step_time = timer() - self.clock_start -+ self.clock_start = timer() - - #Check elapsed timed - if self.time_limit_activated: -- if self.time_limit-(time()-self.time_integration_start) < 0.0: -+ if self.time_limit-(timer()-self.time_integration_start) < 0.0: - raise TimeLimitExceeded("The time limit was exceeded at integration time %.8E."%self.t) - - if self.display_progress: -- if (time() - self.time_integration_start) > self.display_counter*10: -+ if (timer() - self.time_integration_start) > self.display_counter*10: - self.display_counter += 1 - - sys.stdout.write(" Integrator time: %e" % self.t) +
+ #Store the elapsed time for a single step
+ if self.options["clock_step"]:
+- self.elapsed_step_time = clock() - self.clock_start
+- self.clock_start = clock()
++ self.elapsed_step_time = timer() - self.clock_start
++ self.clock_start = timer()
+
+ #Check elapsed timed
+ if self.time_limit_activated:
+- if self.time_limit-(time()-self.time_integration_start) < 0.0:
++ if self.time_limit-(timer()-self.time_integration_start) < 0.0:
+ raise TimeLimitExceeded("The time limit was exceeded at integration time %.8E."%self.t)
+
+ if self.display_progress:
+- if (time() - self.time_integration_start) > self.display_counter*10:
++ if (timer() - self.time_integration_start) > self.display_counter*10:
+ self.display_counter += 1
+
+ sys.stdout.write(" Integrator time: %e" % self.t)
Index: src/ode.pyx =================================================================== --- assimulo/ode.pyx (revision 836) +++ assimulo/ode.pyx (revision 837) @@ -17,7 +17,7 @@ - - import numpy as N - cimport numpy as N --import time -+from timeit import default_timer as timer - import pylab as P - - import itertools +
+ import numpy as N
+ cimport numpy as N
+-import time
++from timeit import default_timer as timer
+ import pylab as P
+
+ import itertools
@@ -283,13 +283,13 @@ - self.initialize() - - #Start of simulation, start the clock -- time_start = time.clock() -+ time_start = timer() - - #Start the simulation - self._simulate(t0, tfinal, output_list, REPORT_CONTINUOUSLY, INTERPOLATE_OUTPUT, TIME_EVENT) - - #End of simulation, stop the clock -- time_stop = time.clock() -+ time_stop = timer() - - #Simulation complete, call finalize - self.finalize() + self.initialize()
+
+ #Start of simulation, start the clock
+- time_start = time.clock()
++ time_start = timer()
+
+ #Start the simulation
+ self._simulate(t0, tfinal, output_list, REPORT_CONTINUOUSLY, INTERPOLATE_OUTPUT, TIME_EVENT)
+
+ #End of simulation, stop the clock
+- time_stop = time.clock()
++ time_stop = timer()
+
+ #Simulation complete, call finalize
+ self.finalize()
Index: src/kinsol.py =================================================================== --- assimulo/kinsol.py (revision 836) @@ -109,83 +109,83 @@ Index: src/explicit_ode.pyx --- assimulo/explicit_ode.pyx (revision 836) +++ assimulo/explicit_ode.pyx (revision 837) @@ -25,7 +25,7 @@ - cimport numpy as N - - from exception import * --from time import clock, time -+from timeit import default_timer as timer - - include "constants.pxi" #Includes the constants (textual include) - + cimport numpy as N
+
+ from exception import *
+-from time import clock, time
++from timeit import default_timer as timer
+
+ include "constants.pxi" #Includes the constants (textual include)
+
@@ -170,7 +170,7 @@ - output_index = 0 - - self.time_limit_activated = 1 if self.time_limit > 0 else 0 -- self.time_integration_start = time() -+ self.time_integration_start = timer() - - while (flag == ID_COMPLETE and tevent == tfinal) is False and (self.t-eps > tfinal) if backward else (self.t+eps < tfinal): - + output_index = 0
+
+ self.time_limit_activated = 1 if self.time_limit > 0 else 0
+- self.time_integration_start = time()
++ self.time_integration_start = timer()
+
+ while (flag == ID_COMPLETE and tevent == tfinal) is False and (self.t-eps > tfinal) if backward else (self.t+eps < tfinal):
+
@@ -183,7 +183,7 @@ - - #Initialize the clock, enabling storing elapsed time for each step - if REPORT_CONTINUOUSLY and self.options["clock_step"]: -- self.clock_start = clock() -+ self.clock_start = timer() - - flag, tlist, ylist = self.integrate(self.t, self.y, tevent, opts) - +
+ #Initialize the clock, enabling storing elapsed time for each step
+ if REPORT_CONTINUOUSLY and self.options["clock_step"]:
+- self.clock_start = clock()
++ self.clock_start = timer()
+
+ flag, tlist, ylist = self.integrate(self.t, self.y, tevent, opts)
+
@@ -248,12 +248,12 @@ - - #Store the elapsed time for a single step - if self.options["clock_step"]: -- self.elapsed_step_time = clock() - self.clock_start -- self.clock_start = clock() -+ self.elapsed_step_time = timer() - self.clock_start -+ self.clock_start = timer() - - #Check elapsed timed - if self.time_limit_activated: -- if self.time_limit-(time()-self.time_integration_start) < 0.0: -+ if self.time_limit-(timer()-self.time_integration_start) < 0.0: - raise TimeLimitExceeded("The time limit was exceeded at integration time %.8E."%self.t) - - self.chattering_clear_counter += 1 +
+ #Store the elapsed time for a single step
+ if self.options["clock_step"]:
+- self.elapsed_step_time = clock() - self.clock_start
+- self.clock_start = clock()
++ self.elapsed_step_time = timer() - self.clock_start
++ self.clock_start = timer()
+
+ #Check elapsed timed
+ if self.time_limit_activated:
+- if self.time_limit-(time()-self.time_integration_start) < 0.0:
++ if self.time_limit-(timer()-self.time_integration_start) < 0.0:
+ raise TimeLimitExceeded("The time limit was exceeded at integration time %.8E."%self.t)
+
+ self.chattering_clear_counter += 1
@@ -262,7 +262,7 @@ - self.chattering_ok_print = 1 - - if self.display_progress: -- if (time() - self.time_integration_start) > self.display_counter*10: -+ if (timer() - self.time_integration_start) > self.display_counter*10: - self.display_counter += 1 - - sys.stdout.write(" Integrator time: %e" % self.t) + self.chattering_ok_print = 1
+
+ if self.display_progress:
+- if (time() - self.time_integration_start) > self.display_counter*10:
++ if (timer() - self.time_integration_start) > self.display_counter*10:
+ self.display_counter += 1
+
+ sys.stdout.write(" Integrator time: %e" % self.t)
Index: src/algebraic.pyx =================================================================== --- assimulo/algebraic.pyx (revision 836) +++ assimulo/algebraic.pyx (revision 837) @@ -17,7 +17,7 @@ - - import numpy as N - cimport numpy as N --import time -+from timeit import default_timer as timer - import pylab as P - - import itertools +
+ import numpy as N
+ cimport numpy as N
+-import time
++from timeit import default_timer as timer
+ import pylab as P
+
+ import itertools
@@ -102,13 +102,13 @@ - self.initialize() - - #Start of simulation, start the clock -- time_start = time.clock() -+ time_start = timer() - - #Start the simulation - self.y = self._solve(y) - - #End of simulation, stop the clock -- time_stop = time.clock() -+ time_stop = timer() - - #Simulation complete, call finalize - self.finalize() + self.initialize()
+
+ #Start of simulation, start the clock
+- time_start = time.clock()
++ time_start = timer()
+
+ #Start the simulation
+ self.y = self._solve(y)
+
+ #End of simulation, stop the clock
+- time_stop = time.clock()
++ time_stop = timer()
+
+ #Simulation complete, call finalize
+ self.finalize()
diff --git a/r838.patch b/r838.patch index 167f6f8b2c80..89104386f156 100644 --- a/r838.patch +++ b/r838.patch @@ -24,93 +24,93 @@ Index: src/solvers/__init__.py --- assimulo/solvers/__init__.py (revision 837) +++ assimulo/solvers/__init__.py (revision 838) @@ -18,12 +18,14 @@ - __all__ = ["euler","radau5","sundials","runge_kutta","rosenbrock", - "glimda","odepack","radar5","dasp3","odassl"] - -+import sys -+ - #Import all the solvers from the different modules - try: - from .euler import ExplicitEuler - from .euler import ImplicitEuler - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .radau5 import Radau5ODE - from .radau5 import Radau5DAE + __all__ = ["euler","radau5","sundials","runge_kutta","rosenbrock",
+ "glimda","odepack","radar5","dasp3","odassl"]
+
++import sys
++
+ #Import all the solvers from the different modules
+ try:
+ from .euler import ExplicitEuler
+ from .euler import ImplicitEuler
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .radau5 import Radau5ODE
+ from .radau5 import Radau5DAE
@@ -30,43 +32,43 @@ - from .radau5 import _Radau5ODE - from .radau5 import _Radau5DAE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .sundials import IDA - from .sundials import CVode - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .kinsol import KINSOL - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .runge_kutta import RungeKutta34 - from .runge_kutta import RungeKutta4 - from .runge_kutta import Dopri5 - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .rosenbrock import RodasODE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .odassl import ODASSL - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .odepack import LSODAR - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .radar5 import Radar5ODE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .dasp3 import DASP3ODE - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) - try: - from .glimda import GLIMDA - except ImportError as ie: -- print("Could not find {}".format(ie.args[0].split("'")[1])) -+ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1])) + from .radau5 import _Radau5ODE
+ from .radau5 import _Radau5DAE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .sundials import IDA
+ from .sundials import CVode
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .kinsol import KINSOL
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .runge_kutta import RungeKutta34
+ from .runge_kutta import RungeKutta4
+ from .runge_kutta import Dopri5
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .rosenbrock import RodasODE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .odassl import ODASSL
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .odepack import LSODAR
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .radar5 import Radar5ODE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .dasp3 import DASP3ODE
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
+ try:
+ from .glimda import GLIMDA
+ except ImportError as ie:
+- print("Could not find {}".format(ie.args[0].split("'")[1]))
++ sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
Index: src/solvers/odepack.py =================================================================== --- assimulo/solvers/odepack.py (revision 837) +++ assimulo/solvers/odepack.py (revision 838) @@ -17,6 +17,7 @@ - - import numpy as N - import scipy.linalg as Sc -+import sys - from assimulo.exception import * - from assimulo.ode import * - import logging +
+ import numpy as N
+ import scipy.linalg as Sc
++import sys
+ from assimulo.exception import *
+ from assimulo.ode import *
+ import logging
@@ -27,7 +28,7 @@ - from assimulo.lib.odepack import dlsodar, dcfode, dintdy - from assimulo.lib.odepack import set_lsod_common, get_lsod_common - except ImportError: -- print("Could not find ODEPACK functions") -+ sys.stderr.write("Could not find ODEPACK functions.\n") - - # Real work array - class common_like(object): + from assimulo.lib.odepack import dlsodar, dcfode, dintdy
+ from assimulo.lib.odepack import set_lsod_common, get_lsod_common
+ except ImportError:
+- print("Could not find ODEPACK functions")
++ sys.stderr.write("Could not find ODEPACK functions.\n")
+
+ # Real work array
+ class common_like(object):
diff --git a/r840.patch b/r840.patch index 0c8019ea95c1..eede3cd2b80e 100644 --- a/r840.patch +++ b/r840.patch @@ -141,66 +141,66 @@ Index: src/lib/sundials_callbacks.pxi --- assimulo/lib/sundials_callbacks.pxi (revision 839) +++ assimulo/lib/sundials_callbacks.pxi (revision 840) @@ -110,9 +110,14 @@ - cdef int ret_nnz - cdef int dim = Jacobian.N - cdef realtype* data = Jacobian.data -- cdef int* rowvals = Jacobian.rowvals -- cdef int* colptrs = Jacobian.colptrs - -+ 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; + cdef int ret_nnz
+ cdef int dim = Jacobian.N
+ cdef realtype* data = Jacobian.data
+- cdef int* rowvals = Jacobian.rowvals
+- cdef int* colptrs = Jacobian.colptrs
+
++ 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;
@@ -136,9 +141,9 @@ - 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: -+ 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] + 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:
++ 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]
Index: src/lib/sundials_includes.pxd =================================================================== --- assimulo/lib/sundials_includes.pxd (revision 839) +++ assimulo/lib/sundials_includes.pxd (revision 840) @@ -93,13 +93,29 @@ - ctypedef _DlsMat* DlsMat - cdef realtype* DENSE_COL(DlsMat A, int j) - --IF SUNDIALS_VERSION >= (2,6,0): -+IF SUNDIALS_VERSION >= (2,6,3): - cdef extern from "sundials/sundials_sparse.h": - cdef struct _SlsMat: - int M - int N - int NNZ -+ int NP - realtype *data -+ int sparsetype -+ int *indexvals -+ int *indexptrs -+ int **rowvals -+ int **colptrs -+ int **colvals -+ int **rowptrs -+ ctypedef _SlsMat* SlsMat -+ELIF SUNDIALS_VERSION >= (2,6,0): -+ cdef extern from "sundials/sundials_sparse.h": -+ cdef struct _SlsMat: -+ int M -+ int N -+ int NNZ -+ realtype *data - int *rowvals - int *colptrs - ctypedef _SlsMat* SlsMat + ctypedef _DlsMat* DlsMat
+ cdef realtype* DENSE_COL(DlsMat A, int j)
+
+-IF SUNDIALS_VERSION >= (2,6,0):
++IF SUNDIALS_VERSION >= (2,6,3):
+ cdef extern from "sundials/sundials_sparse.h":
+ cdef struct _SlsMat:
+ int M
+ int N
+ int NNZ
++ int NP
+ realtype *data
++ int sparsetype
++ int *indexvals
++ int *indexptrs
++ int **rowvals
++ int **colptrs
++ int **colvals
++ int **rowptrs
++ ctypedef _SlsMat* SlsMat
++ELIF SUNDIALS_VERSION >= (2,6,0):
++ cdef extern from "sundials/sundials_sparse.h":
++ cdef struct _SlsMat:
++ int M
++ int N
++ int NNZ
++ realtype *data
+ int *rowvals
+ int *colptrs
+ ctypedef _SlsMat* SlsMat
diff --git a/r845.patch b/r845.patch index 4069af48895c..dec23fc92e53 100644 --- a/r845.patch +++ b/r845.patch @@ -92,161 +92,161 @@ 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 + 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") + 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 - + 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 +
+ 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 + 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"): - """ + 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"]) +
+ 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) - +
+ 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) @@ -532,1010 +532,1010 @@ 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 - + 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: + 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: + 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) + 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, + 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) + 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 +
+ #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) + # 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) + 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 - #========== + 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) + 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) + 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): - """ + 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 - + 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 - + """
+ 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): - """ + 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) + """
+ 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 + 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) @@ -1672,18 +1672,18 @@ 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) - +
+ 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)
+
|