diff options
author | Michel Zou | 2019-01-08 19:34:11 +0100 |
---|---|---|
committer | Michel Zou | 2019-01-08 19:34:11 +0100 |
commit | 4a75591f31b4d04c24e2d23ff276ef1bdd3338b2 (patch) | |
tree | 85fcfb585086c88dd1ca123175bddd3c394450bd | |
parent | 4e155d7d2c6215da49a3b8a36542c04ec6b07d15 (diff) | |
download | aur-4a75591f31b4d04c24e2d23ff276ef1bdd3338b2.tar.gz |
3.0
-rw-r--r-- | .SRCINFO | 27 | ||||
-rw-r--r-- | PKGBUILD | 24 | ||||
-rw-r--r-- | r831.patch | 64 | ||||
-rw-r--r-- | r833.patch | 85 | ||||
-rw-r--r-- | r836.patch | 92 | ||||
-rw-r--r-- | r837.patch | 191 | ||||
-rw-r--r-- | r838.patch | 116 | ||||
-rw-r--r-- | r839.patch | 5467 | ||||
-rw-r--r-- | r840.patch | 206 | ||||
-rw-r--r-- | r841.patch | 13 | ||||
-rw-r--r-- | r845.patch | 1689 | ||||
-rw-r--r-- | sundials4.patch | 26 |
12 files changed, 40 insertions, 7960 deletions
@@ -1,9 +1,7 @@ -# Generated by mksrcinfo v8 -# Wed May 16 13:43:51 UTC 2018 pkgbase = python-assimulo pkgdesc = A package for solving ordinary differential equations and differential algebraic equations - pkgver = 2.9 - pkgrel = 4 + pkgver = 3.0 + pkgrel = 1 url = http://www.jmodelica.org/assimulo arch = i686 arch = x86_64 @@ -15,24 +13,9 @@ pkgbase = python-assimulo makedepends = gcc-fortran makedepends = sundials makedepends = lapack - source = https://pypi.python.org/packages/4c/c0/19a54949817204313efff9f83f1e4a247edebed0a1cc5a317a95d3f374ae/Assimulo-2.9.zip - source = r831.patch - source = r833.patch - source = r836.patch - source = r837.patch - source = r838.patch - source = r839.patch - source = r840.patch - source = r841.patch - source = r845.patch - md5sums = 3f28fd98011d2ec7a01703a1ef1dff45 - md5sums = SKIP - md5sums = SKIP - md5sums = SKIP - md5sums = SKIP - md5sums = SKIP - md5sums = SKIP - md5sums = SKIP + makedepends = subversion + source = svn+https://svn.jmodelica.org/assimulo/tags/Assimulo-3.0 + source = sundials4.patch md5sums = SKIP md5sums = SKIP @@ -1,28 +1,22 @@ pkgbase=python-assimulo pkgname=('python-assimulo' 'python2-assimulo') -pkgver=2.9 -pkgrel=4 +pkgver=3.0 +pkgrel=1 pkgdesc="A package for solving ordinary differential equations and differential algebraic equations" url="http://www.jmodelica.org/assimulo" arch=('i686' 'x86_64') license=('LGPL') -makedepends=('python-setuptools' 'python2-setuptools' 'cython' 'cython2' 'gcc-fortran' 'sundials' 'lapack') -source=("https://pypi.python.org/packages/4c/c0/19a54949817204313efff9f83f1e4a247edebed0a1cc5a317a95d3f374ae/Assimulo-2.9.zip" r831.patch r833.patch r836.patch r837.patch r838.patch r839.patch r840.patch r841.patch r845.patch) -md5sums=('3f28fd98011d2ec7a01703a1ef1dff45' SKIP SKIP SKIP SKIP SKIP SKIP SKIP SKIP SKIP) +makedepends=('python-setuptools' 'python2-setuptools' 'cython' 'cython2' 'gcc-fortran' 'sundials' 'lapack' 'subversion') +# PyPI archive is missing .pyf files +# source=("https://pypi.io/packages/source/A/Assimulo/Assimulo-${pkgver}.tar.gz") +source=("svn+https://svn.jmodelica.org/assimulo/tags/Assimulo-${pkgver}" sundials4.patch) +md5sums=('SKIP' SKIP ) prepare() { cd "${srcdir}"/Assimulo-$pkgver - # https://trac.jmodelica.org/assimulo/changeset/845 - patch -p0 -i "${srcdir}"/r831.patch - patch -p0 -i "${srcdir}"/r833.patch - patch -p0 -i "${srcdir}"/r836.patch - patch -p0 -i "${srcdir}"/r837.patch - patch -p0 -i "${srcdir}"/r838.patch - patch -p0 -i "${srcdir}"/r839.patch - patch -p0 -i "${srcdir}"/r840.patch - patch -p0 -i "${srcdir}"/r841.patch - patch -p0 -i "${srcdir}"/r845.patch + # try to build with sundials 4.0 + patch -p0 -i "${srcdir}"/sundials4.patch } build() { diff --git a/r831.patch b/r831.patch deleted file mode 100644 index 1215207005ad..000000000000 --- a/r831.patch +++ /dev/null @@ -1,64 +0,0 @@ -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:]))
diff --git a/r833.patch b/r833.patch deleted file mode 100644 index 7385443600c1..000000000000 --- a/r833.patch +++ /dev/null @@ -1,85 +0,0 @@ -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
-
-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
-@@ -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]))
diff --git a/r836.patch b/r836.patch deleted file mode 100644 index b3cc0a96ba5d..000000000000 --- a/r836.patch +++ /dev/null @@ -1,92 +0,0 @@ -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
-@@ -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;
-@@ -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)
-@@ -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:
-@@ -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:
-@@ -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
diff --git a/r837.patch b/r837.patch deleted file mode 100644 index e5fdd761d3a3..000000000000 --- a/r837.patch +++ /dev/null @@ -1,191 +0,0 @@ -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
-@@ -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):
-
-@@ -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)
-
-@@ -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)
-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
-@@ -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()
-Index: src/kinsol.py -=================================================================== ---- assimulo/kinsol.py (revision 836) -+++ assimulo/kinsol.py (revision 837) -@@ -26,7 +26,7 @@ - import pylab as P - import operator as O - import re --import time -+from timeit import default_timer as timer - from assimulo.problem_algebraic import * - - class KINSOL_Exception(Exception): -@@ -290,9 +290,9 @@ - self.solver.KINSOL_init(self.func,self.x0,self.dim,jac,self.constraints,self.use_sparse,self.verbosity,self.norm_of_res,self.reg_param,self.fscale) - else: - self.solver.KINSOL_init(self.func,self.x0,self.dim,jac,self.constraints,self.use_sparse,self.verbosity,self.norm_of_res,self.reg_param,None) -- start = time.clock() -+ start = timer() - res = self.solver.KINSOL_solve(not self._use_ls) -- stop = time.clock() -+ stop = timer() - self.exec_time += (stop - start) - solved = True - except KINError as error: -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)
-
-@@ -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):
-
-@@ -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)
-
-@@ -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
-@@ -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)
-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
-@@ -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()
diff --git a/r838.patch b/r838.patch deleted file mode 100644 index 89104386f156..000000000000 --- a/r838.patch +++ /dev/null @@ -1,116 +0,0 @@ -Index: src/solvers/glimda.py -=================================================================== ---- assimulo/solvers/glimda.py (revision 837) -+++ assimulo/solvers/glimda.py (revision 838) -@@ -16,6 +16,7 @@ - # along with this program. If not, see <http://www.gnu.org/licenses/>. - - import numpy as N -+import sys - - from assimulo.exception import * - from assimulo.ode import * -@@ -25,7 +26,7 @@ - try: - from assimulo.lib.glimda import glimda - except ImportError: -- print("Could not find GLIMDA") -+ sys.stderr.write("Could not find GLIMDA.\n") - - class GLIMDA(Implicit_ODE): - """ -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
-@@ -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]))
-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
-@@ -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):
diff --git a/r839.patch b/r839.patch deleted file mode 100644 index 9fde401d3862..000000000000 --- a/r839.patch +++ /dev/null @@ -1,5467 +0,0 @@ -Index: setup.py -=================================================================== ---- setup.py (revision 838) -+++ setup.py (revision 839) -@@ -387,39 +387,7 @@ - if self.platform == "mac": - el.extra_compile_args += ["-Wno-error=return-type"] - el.extra_compile_args += self.flag_32bit + self.extra_c_flags -- -- if self.with_SUNDIALS: -- cordir_KINSOL_wSLU = os.path.join(self.assimulo_lib,'sundials_kinsol_core_wSLU.pyx') -- cordir_KINSOL = os.path.join(self.assimulo_lib,'sundials_kinsol_core.pyx') -- -- cordir_KINSOL_jmod_wSLU = os.path.join(self.assimulo_lib,'kinsol_jmod_wSLU.c') -- cordir_KINSOL_jmod = os.path.join(self.assimulo_lib,'kinsol_jmod.c') -- -- cordir_kinpinv = os.path.join(self.assimulo_lib,'kinpinv.c') -- cordir_kinslug = os.path.join(self.assimulo_lib,'kinslug.c') -- cordir_reg_routines = os.path.join(self.assimulo_lib,'reg_routines.c') -- if self.with_SLU: -- ext_list = ext_list + cythonize([cordir_KINSOL_wSLU], include_path=[".","assimulo",self.assimulo_lib]) -- ext_list[-1].sources += [cordir_KINSOL_jmod_wSLU,cordir_kinpinv,cordir_kinslug,cordir_reg_routines] -- ext_list[-1].include_dirs = [np.get_include(), self.SLUincdir, self.incdirs] -- ext_list[-1].library_dirs = [self.libdirs,self.SLUlibdir,self.BLASdir] -- ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial", "superlu_4.1",self.BLASname,'gfortran'] -- else: -- ext_list = ext_list + cythonize([cordir_KINSOL])#, include_path=[".","assimulo",self.assimulo_lib]) -- ext_list[-1].sources += [cordir_KINSOL_jmod,cordir_kinpinv] -- ext_list[-1].include_dirs = [np.get_include(), self.incdirs] -- ext_list[-1].library_dirs = [self.libdirs] -- ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial"] -- if self.SUNDIALS_version > (2,5,0): -- ext_list[-1].define_macros.append(("SUNDIALS_26", 1)) -- if self.debug_flag: -- ext_list[-1].extra_compile_args = ["-g", "-fno-strict-aliasing"] -- else: -- ext_list[-1].extra_compile_args = ["-O2", "-fno-strict-aliasing"] -- if self.platform == "mac": -- ext_list[-1].extra_compile_args += ["-Wno-error=return-type"] -- ext_list[-1].extra_compile_args += self.flag_32bit + self.extra_c_flags -- -+ - for el in ext_list: - if self.is_python3: - el.cython_directives = {"language_level": 3} -Index: src/kinsol.py -=================================================================== ---- assimulo/kinsol.py (revision 838) -+++ assimulo/kinsol.py (nonexistent) -@@ -1,350 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# This program is distributed in the hope that it will be useful, --# but WITHOUT ANY WARRANTY; without even the implied warranty of --# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --# GNU Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. --try: -- from .lib.sundials_kinsol_core_wSLU import KINSOL_wrap, KINError --except: -- from .lib.sundials_kinsol_core import KINSOL_wrap, KINError -- --from scipy.linalg import pinv2 --from scipy.optimize import fminbound --from numpy.linalg import solve --import numpy as N --import pylab as P --import operator as O --import re --from timeit import default_timer as timer --from assimulo.problem_algebraic import * -- --class KINSOL_Exception(Exception): -- pass -- -- --class KINSOL: -- -- def __init__(self,problem): -- """ -- Create the solver -- -- Parameters:: -- problem-- -- instance of ProblemAlgebraic found in problem_algebraic.py -- """ -- -- self.solver = KINSOL_wrap() -- -- # extract info from problem -- self.problem = problem -- if hasattr(self.problem,'_x0'): -- try: -- _x0 = self.problem.get_x0() -- except ProblemAlg_Exception: -- raise KINSOL_Exception("Problem has not implemented method 'get_x0'") -- else: -- raise KINSOL_Exception("Problem has no instance '_x0'") -- -- # calculate dimension -- try: -- if isinstance(_x0, int) or isinstance(_x0, float): -- self.x0 = [_x0] -- else: -- self.x0 = _x0 -- self.dim = len([N.array(self.x0, dtype=float)][0]) -- except ValueError: -- raise KINSOL_Exception("Initial guess must be a Numpy.array with either ints or floats.") -- -- # check for functions and test them -- try: -- tmp = self.problem.f(self.x0) -- self.norm_of_res = 10000*N.linalg.norm(self.x0) -- #self.norm_of_res = 1e20 -- self.func = self.problem.f -- except ProblemAlg_Exception: -- raise KINSOL_Exception("Problem has not implemented method 'f'") -- except IndexError: -- raise KINSOL_Exception("Problem has mismatching dimensions of f and initial guess") -- -- # Calculate 'stupid' scaling -- res = self.problem.f(self.x0) -- self.fscale = N.ones(self.dim) -- -- for i,val in zip(N.arange(self.dim),res): -- if val > 10: -- self.fscale[i] = 1.0/val -- -- # check for constraints and test them -- broken_constraints = [] -- if hasattr(self.problem, 'get_constraints'): -- self.constraints = self.problem.get_constraints() -- -- if self.constraints is not None: -- # test if constraints are of correct type -- if type(self.constraints).__name__ != 'ndarray': -- raise KINSOL_Exception("Constraints must be of type numpy.ndarray") -- -- if len(self.constraints) != len(self.x0): -- raise KINSOL_Exception("Constraints must have same length as x0") -- # Test if initial guess x0 is consistant with constraints -- for c,xi,i in zip(self.constraints,self.x0,N.arange(0,self.x0.__len__())): -- if re.search('float',type(c).__name__) is None: -- print("Type problem with: ", c, type(c).__name__) -- raise KINSOL_Exception("Constraints must contain floats.") -- if abs(c) > 2: -- raise KINSOL_Exception("Entries in constraint vector must be between -2 and 2, see documentation.") -- -- if c != 0.0: -- if c == 1.0: -- if xi < 0.0: -- broken_constraints.append(i) -- elif c == 2.0: -- if xi <= 0.0: -- broken_constraints.append(i) -- elif c == -1.0: -- if xi > 0.0: -- broken_constraints.append(i) -- elif c == -2.0: -- if xi >= 0.0: -- broken_constraints.append(i) -- else: -- raise KINSOL_Exception("Constraint vector contains illegal elements.") -- -- else: -- self.constraints = None -- -- if broken_constraints != []: -- print("Variables breaking initial constraint: ") -- for i in broken_constraints: -- self.problem.print_var_info(i) -- -- raise KINSOL_Exception("Initial guess does not fulfill applied constraints.") -- -- -- if hasattr(self.problem, 'check_constraints'): -- self.check_with_model = True -- else: -- self.check_with_model = False -- -- self._use_jac = True -- self.reg_count = 0 -- self.lin_count = 0 -- self.verbosity = 0 -- self.max_reg = 2.0 -- self.use_sparse = False -- self.reg_param = 0.0 -- self.exec_time = 0 -- self._use_ls = True -- self._use_fscale = False -- -- def set_jac_usage(self,use_jac): -- """ -- Set whether to use the jacobian supplied by the model -- or if we are to calculate it numericaly -- -- Parameters:: -- -- use_jac -- -- Boolean set to True if the jacobian is to be -- supplied by the model -- -- """ -- if type(use_jac).__name__ == 'bool': -- self._use_jac = use_jac -- else: -- raise KINSOL_Exception("The variable sent to 'set_jac_usage' must be a boolean.") -- -- def use_LineSearch(self,use_ls): -- """ -- Set whether the solver starts using the Linesearch Algorithm -- or not -- -- Parameters:: -- -- use_ls -- -- Boolean set to True if LineSearch is to -- be used. -- -- """ -- if type(use_ls).__name__ == 'bool': -- self._use_ls = use_ls -- else: -- raise KINSOL_Exception("The variable sent to 'use_LineSearch' must be a boolean.") -- -- def use_fscale(self,use_fscale): -- """ -- Set whether the solver starts using the Linesearch Algorithm -- or not -- -- Parameters:: -- -- use_ls -- -- Boolean set to True if LineSearch is to -- be used. -- -- """ -- if type(use_fscale).__name__ == 'bool': -- self._use_fscale = use_fscale -- else: -- raise KINSOL_Exception("The variable sent to 'use_fscale' must be a boolean.") -- -- -- def set_verbosity(self,verbosity): -- """ -- Method used to set the verbosity of KINSOL -- -- Parameters:: -- -- verbosity -- -- Integer set to one of 0,1,2,3 where 0 is no info and 3 the -- most info. For more information please see the documentation for KINSOL. -- """ -- type_name = type(verbosity).__name__ -- if re.search('int',type_name) is not None: -- -- # It is an integer, tes bounds -- if verbosity < 4 and verbosity > -1: -- self.verbosity = verbosity -- else: -- raise KINSOL_Exception("The variable sent to 'set_verbosity' must be either 0, 1, 2 or 3.") -- else: -- raise KINSOL_Exception("The variable sent to 'set_verbosity' must be an integer.") -- -- def set_sparsity(self,use_sparse): -- """ -- Method used to set if the problem should be treated as sparse by -- the linear solver in KINSOL. If the problem supplied has not implemented -- a method sparse_jac an exception will be thrown -- -- Parameters:: -- -- use_sparse -- -- Boolean set to True if the problem is to be treated as -- sparse and False otherwise. -- -- """ -- -- if hasattr(self.problem,'sparse_jac'): -- self.use_sparse = use_sparse -- else: -- raise KINSOL_Exception("The problem must have implemented a method 'sparse_jac' for sparsity to by used.") -- -- def set_reg_param(self,reg_param): -- """ -- Method used to set the regularization parameter. -- If set to zero the parameter will be set by the solver. -- -- Parameters:: -- -- reg_param -- -- Float larger or equal to zero. -- -- """ -- -- type_name = type(reg_param).__name__ -- if re.search('float',type_name) is not None: -- if reg_param < 0: -- raise KINSOL_Exception("Value sent to set_reg_param must be equal to or larger than zero.") -- else: -- self.reg_param = reg_param -- else: -- raise KINSOL_Exception("The variable sent to 'set_reg_param' must be a float.") -- -- def solve(self): -- """ -- Function called when solving function rhs_fct -- -- """ -- # check for jacobian and set it if present and to be used -- if self.use_sparse: -- if self._use_jac and hasattr(self.problem,'sparse_jac'): -- jac = self.problem.sparse_jac -- else: -- jac = None -- else: -- if self._use_jac and hasattr(self.problem,'jac'): -- jac = self.problem.jac -- else: -- jac = None -- -- # Initialize solver and solve -- -- solved = False -- local_min = False -- -- res = N.zeros(self.x0.__len__()) -- while (not solved) and self.reg_count < 2: -- try: -- if self._use_fscale: -- self.solver.KINSOL_init(self.func,self.x0,self.dim,jac,self.constraints,self.use_sparse,self.verbosity,self.norm_of_res,self.reg_param,self.fscale) -- else: -- self.solver.KINSOL_init(self.func,self.x0,self.dim,jac,self.constraints,self.use_sparse,self.verbosity,self.norm_of_res,self.reg_param,None) -- start = timer() -- res = self.solver.KINSOL_solve(not self._use_ls) -- stop = timer() -- self.exec_time += (stop - start) -- solved = True -- except KINError as error: -- if error.value == 42: -- # Try the heuristic -- if hasattr(self.problem, 'get_heuristic_x0'): -- print("----------------------------------------------------") -- print(" Solver stuck with zero step-length.") -- print("----------------------------------------------------") -- print("The following variables have start value zero") -- print("and min set to zero causing the zero step-lenght.") -- print("These settings are either set by default or by user.") -- print("") -- -- self.x0 = self.problem.get_heuristic_x0() -- self.reg_count += 1 -- -- print("") -- print("This setting (start and min to zero) can often") -- print("cause problem when initializing the system. ") -- print("") -- print("To avoid this the above variables have") -- print("their start attributes reset to one.") -- print("") -- print("Trying to solve the system again...") -- else: -- raise KINSOL_Exception("Regularization failed due to constraints, tried getting heuristic initial guess but failed.") -- -- -- elif (error.value == 2): -- print("---------------------------------------------------------") -- print("") -- print(" !!! WARNING !!!") -- print("") -- print(" KINSOL has returned a result but the algorithm has converged") -- print(" to a local minima, the initial values are NOT consistant!") -- print("") -- print("---------------------------------------------------------") -- solved = True -- local_min = True -- else: -- # Other error, send onward as exception -- self.problem.check_constraints(res) -- raise KINSOL_Exception(error.msg[error.value]) -- -- if not solved: -- self.solver.Free_KINSOL() -- raise KINSOL_Exception("Algorithm exited solution loop without finding a solution, please contact Assimulo support.") -- -- if self.check_with_model: -- self.problem.check_constraints(res) -- if not local_min: -- print("Problem sent to KINSOL solved.") -- -- return res -Index: src/lib/sundials_kinsol_core_wSLU.pxi -=================================================================== ---- assimulo/lib/sundials_kinsol_core_wSLU.pxi (revision 838) -+++ assimulo/lib/sundials_kinsol_core_wSLU.pxi (nonexistent) -@@ -1,293 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# This program is distributed in the hope that it will be useful, --# but WITHOUT ANY WARRANTY; without even the implied warranty of --# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --# GNU Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" -- --import numpy as np --from numpy cimport NPY_DOUBLE, npy_intp, NPY_INT -- --#============================================== --#External definitions from numpy headers --#============================================== --cdef extern from "numpy/arrayobject.h": -- -- ctypedef int npy_intp -- -- ctypedef extern class numpy.dtype [object PyArray_Descr]: -- cdef int type_num, elsize, alignment -- cdef char type, kind, byteorder, hasobject -- cdef object fields, typeobj -- -- ctypedef int intp -- ctypedef extern class numpy.ndarray [object PyArrayObject]: -- cdef char *data -- cdef int nd -- cdef intp *dimensions -- cdef intp *strides -- cdef int flags -- -- cdef object PyArray_SimpleNew(int nd, npy_intp* dims, int typenum) -- cdef object PyArray_SimpleNewFromData(int nd, npy_intp *dims, -- int typenum, void *data) -- void import_array() -- void* PyArray_GetPtr(ndarray aobj, npy_intp* ind) -- void *PyArray_DATA(ndarray aobj) -- --import_array() -- --#============================================== --#External definitions from SuperLU headers --#============================================== --cdef extern from "slu_ddefs.h": -- ctypedef int int_t -- -- ctypedef enum Stype_t: -- SLU_NC -- ctypedef enum Dtype_t: -- SLU_D -- ctypedef enum Mtype_t: -- SLU_GE -- ctypedef struct SuperMatrix: -- Stype_t Stype -- Dtype_t Dtype -- Mtype_t Mtype -- int_t nrow -- int_t ncol -- void *Store -- void dCreate_CompCol_Matrix(SuperMatrix *, int, int, int, double *,int *, int *, Stype_t, Dtype_t, Mtype_t) --#============================================== --# C headers --#============================================== --cdef extern from "string.h": -- void *memcpy(void *s1, void *s2, int n) --cdef extern from "stdio.h": -- int printf (char *__format, ...) --cdef extern from "stdlib.h": -- ctypedef unsigned long size_t -- void *calloc(size_t nelem, size_t elsize) -- --#============================================== --#External definitions from Sundials headers --#============================================== -- --# import data types --# basic data types --cdef extern from "sundials/sundials_types.h": -- ctypedef double realtype -- ctypedef bint booleantype -- --# sundials handling of vectors, so far only serial, --# parallel to be implemented later? --cdef extern from "sundials/sundials_nvector.h": -- cdef struct _generic_N_Vector: -- void* content -- ctypedef _generic_N_Vector* N_Vector -- --cdef extern from "nvector/nvector_serial.h": -- cdef struct _N_VectorContent_Serial: -- long int length -- realtype* data -- ctypedef _N_VectorContent_Serial* N_VectorContent_Serial -- N_Vector N_VNew_Serial(long int vec_length) -- void N_VDestroy_Serial(N_Vector v) -- void N_VPrint_Serial(N_Vector v) -- cdef N_Vector N_VMake_Serial(long int vec_length, realtype *v_data) -- N_Vector *N_VCloneVectorArray_Serial(int count, N_Vector w) -- N_Vector *N_VCloneVectorArrayEmpty_Serial(int count, N_Vector w) -- void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v) -- void N_VConst_Serial(realtype c, N_Vector z) -- --#Struct for handling the Jacobian data --cdef extern from "sundials/sundials_direct.h": -- cdef struct _DlsMat: -- int type -- int M -- int N -- int ldim -- int mu -- int ml -- int s_mu -- realtype *data -- int ldata -- realtype **cols -- ctypedef _DlsMat* DlsMat -- cdef realtype* DENSE_COL(DlsMat A, int j) -- --# KINSOL functions and routines --cdef extern from "kinsol/kinsol.h": -- # user defined functions -- ctypedef int (*KINSysFn)(N_Vector uu, N_Vector fval, void *user_data ) -- ctypedef void (*KINErrHandlerFn)(int error_code, char *module, char *function, char *msg, void *user_data) -- ctypedef void (*KINInfoHandlerFn)(char *module, char *function, char *msg, void *user_data) -- # initialization routines -- void *KINCreate() -- int KINInit(void *kinmem, KINSysFn func, N_Vector tmpl) -- -- # optional input spec. functions, -- # for specificationsdocumentation cf. kinsol.h line 218-449 -- int KINSetErrHandlerFn(void *kinmem, KINErrHandlerFn ehfun, void *eh_data) -- int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_data) -- int KINSetUserData(void *kinmem, void *user_data) -- int KINSetPrintLevel(void *kinmemm, int printfl) -- int KINSetNumMaxIters(void *kinmem, long int mxiter) -- int KINSetNoInitSetup(void *kinmem, booleantype noInitSetup) -- int KINSetNoResMon(void *kinmem, booleantype noNNIResMon) -- int KINSetMaxSetupCalls(void *kinmem, long int msbset) -- int KINSetMaxSubSetupCalls(void *kinmem, long int msbsetsub) -- int KINSetEtaForm(void *kinmem, int etachoice) -- int KINSetEtaConstValue(void *kinmem, realtype eta) -- int KINSetEtaParams(void *kinmem, realtype egamma, realtype ealpha) -- int KINSetResMonParams(void *kinmem, realtype omegamin, realtype omegamax) -- int KINSetResMonConstValue(void *kinmem, realtype omegaconst) -- int KINSetNoMinEps(void *kinmem, booleantype noMinEps) -- int KINSetMaxNewtonStep(void *kinmem, realtype mxnewtstep) -- int KINSetMaxBetaFails(void *kinmem, long int mxnbcf) -- int KINSetRelErrFunc(void *kinmem, realtype relfunc) -- int KINSetFuncNormTol(void *kinmem, realtype fnormtol) -- int KINSetScaledStepTol(void *kinmem, realtype scsteptol) -- int KINSetConstraints(void *kinmem, N_Vector constraints) -- int KINSetSysFunc(void *kinmem, KINSysFn func) -- -- # solver routine -- int KINSol(void *kinmem, N_Vector uu, int strategy, N_Vector u_scale, N_Vector f_scale) -- -- # optional output routines. -- # Documentation see kinsol.h line 670-735 -- int KINGetWorkSpace(void *kinmem, long int *lenrw, long int *leniw) -- int KINGetNumNonlinSolvIters(void *kinmem, long int *nniters) -- int KINGetNumFuncEvals(void *kinmem, long int *nfevals) -- int KINGetNumBetaCondFails(void *kinmem, long int *nbcfails) -- int KINGetNumBacktrackOps(void *kinmem, long int *nbacktr) -- int KINGetFuncNorm(void *kinmem, realtype *fnorm) -- int KINGetStepLength(void *kinmem, realtype *steplength) -- char *KINGetReturnFlagName(int flag) -- -- # fuction used to deallocate memory used by KINSOL -- void KINFree(void **kinmem) -- --# linear solver modules --cdef extern from "kinsol/kinsol_dense.h": -- int KINDense(void *kinmem, int N) -- --cdef extern from "kinpinv.h": -- int KINPinv(void *kinmem, int N) -- --cdef extern from "kinslug.h": -- int KINSLUG(void *kinmem, int N) -- --# functions used for supplying jacobian, and receiving info from linear solver --cdef extern from "kinsol/kinsol_direct.h": -- # user functions -- ctypedef int (*KINDlsDenseJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user functions to KINSOL -- int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac) -- -- # optional output fcts for linear direct solver -- int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINDlsGetLastFlag(void *kinmem, int *flag) -- char *KINDlsGetReturnFlagName(int flag) -- --cdef extern from "kinsol_jmod_wSLU.h": -- # user functions -- ctypedef int (*KINPinvJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- ctypedef int (*KINSLUGJacFn)(int N, N_Vector u, N_Vector fu, SuperMatrix *J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user jacobian to KINSOL -- int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) -- int KINSLUGSetJacFn(void *kinmem, KINSLUGJacFn jac) -- -- # functions used to set regularizationparameter -- int KINSLUGSetRegParam(void *kinmem, realtype reg_p) -- int KINPinvSetRegParam(void *kinmem, realtype reg_p) -- -- # optional output fcts for linear direct solver -- int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINPinvGetLastFlag(void *kinmem, int *flag) -- char *KINPinvGetReturnFlagName(int flag) -- -- --#========================= --# END SUNDIALS DEFINITIONS --#========================= -- --#=========================== --# Problem information struct --#=========================== --cdef class ProblemData: -- cdef: -- void *RHS # Residual or right-hand-side -- void *JAC -- int dim # Dimension of the problem -- --#================= --# Module functions --#================= -- --cdef inline N_Vector arr2nv(x): -- x=np.array(x) -- cdef long int n = len(x) -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- cdef N_Vector v=N_VNew_Serial(n) -- memcpy((<N_VectorContent_Serial>v.content).data, data_ptr, n*sizeof(double)) -- return v -- --cdef inline nv2arr(N_Vector v): -- cdef long int n = (<N_VectorContent_Serial>v.content).length -- cdef realtype* v_data = (<N_VectorContent_Serial>v.content).data -- cdef long int i -- import_array() -- cdef ndarray x=np.empty(n) -- #cdef ndarray x = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE,v_data) -- memcpy(x.data, v_data, n*sizeof(double)) -- return x -- --cdef inline realtype2arr(realtype *data, int n): -- """ Create new numpy array from realtype*""" -- import_array() -- cdef ndarray[realtype, ndim=1, mode='c'] x=np.empty(n) -- memcpy(x.data, data, n*sizeof(double)) -- return x -- --cdef inline arr2realtype(ndarray x, realtype* out_data, int n): -- """ Create new realtype* from numpy """ -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(double)) -- -- --cdef inline arr2int(ndarray x, int* out_data, int n): -- """ Create new int* from numpy array """ -- cdef ndarray[int, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(int)) -- -Index: src/lib/sundials_kinsol_includes.pxi -=================================================================== ---- assimulo/lib/sundials_kinsol_includes.pxi (revision 838) -+++ assimulo/lib/sundials_kinsol_includes.pxi (nonexistent) -@@ -1,267 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# This program is distributed in the hope that it will be useful, --# but WITHOUT ANY WARRANTY; without even the implied warranty of --# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --# GNU Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" -- --import numpy as np --from numpy cimport NPY_DOUBLE, npy_intp, NPY_INT -- --#============================================== --#External definitions from numpy headers --#============================================== --cdef extern from "numpy/arrayobject.h": -- -- ctypedef int npy_intp -- -- ctypedef extern class numpy.dtype [object PyArray_Descr]: -- cdef int type_num, elsize, alignment -- cdef char type, kind, byteorder, hasobject -- cdef object fields, typeobj -- -- ctypedef int intp -- ctypedef extern class numpy.ndarray [object PyArrayObject]: -- cdef char *data -- cdef int nd -- cdef intp *dimensions -- cdef intp *strides -- cdef int flags -- -- cdef object PyArray_SimpleNew(int nd, npy_intp* dims, int typenum) -- cdef object PyArray_SimpleNewFromData(int nd, npy_intp *dims, -- int typenum, void *data) -- void import_array() -- void* PyArray_GetPtr(ndarray aobj, npy_intp* ind) -- void *PyArray_DATA(ndarray aobj) -- --import_array() -- --#============================================== --# C headers --#============================================== --cdef extern from "string.h": -- void *memcpy(void *s1, void *s2, int n) --cdef extern from "stdio.h": -- int printf (char *__format, ...) --cdef extern from "stdlib.h": -- ctypedef unsigned long size_t -- void *calloc(size_t nelem, size_t elsize) -- --#============================================== --#External definitions from Sundials headers --#============================================== -- --# import data types --# basic data types --cdef extern from "sundials/sundials_types.h": -- ctypedef double realtype -- ctypedef bint booleantype -- --# sundials handling of vectors, so far only serial, --# parallel to be implemented later? --cdef extern from "sundials/sundials_nvector.h": -- cdef struct _generic_N_Vector: -- void* content -- ctypedef _generic_N_Vector* N_Vector -- --cdef extern from "nvector/nvector_serial.h": -- cdef struct _N_VectorContent_Serial: -- long int length -- realtype* data -- ctypedef _N_VectorContent_Serial* N_VectorContent_Serial -- N_Vector N_VNew_Serial(long int vec_length) -- void N_VDestroy_Serial(N_Vector v) -- void N_VPrint_Serial(N_Vector v) -- cdef N_Vector N_VMake_Serial(long int vec_length, realtype *v_data) -- N_Vector *N_VCloneVectorArray_Serial(int count, N_Vector w) -- N_Vector *N_VCloneVectorArrayEmpty_Serial(int count, N_Vector w) -- void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v) -- void N_VConst_Serial(realtype c, N_Vector z) -- --#Struct for handling the Jacobian data --cdef extern from "sundials/sundials_direct.h": -- cdef struct _DlsMat: -- int type -- int M -- int N -- int ldim -- int mu -- int ml -- int s_mu -- realtype *data -- int ldata -- realtype **cols -- ctypedef _DlsMat* DlsMat -- cdef realtype* DENSE_COL(DlsMat A, int j) -- --# KINSOL functions and routines --cdef extern from "kinsol/kinsol.h": -- # user defined functions -- ctypedef int (*KINSysFn)(N_Vector uu, N_Vector fval, void *user_data ) -- ctypedef void (*KINErrHandlerFn)(int error_code, char *module, char *function, char *msg, void *user_data) -- ctypedef void (*KINInfoHandlerFn)(char *module, char *function, char *msg, void *user_data) -- # initialization routines -- void *KINCreate() -- int KINInit(void *kinmem, KINSysFn func, N_Vector tmpl) -- -- # optional input spec. functions, -- # for specificationsdocumentation cf. kinsol.h line 218-449 -- int KINSetErrHandlerFn(void *kinmem, KINErrHandlerFn ehfun, void *eh_data) -- int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_data) -- int KINSetUserData(void *kinmem, void *user_data) -- int KINSetPrintLevel(void *kinmemm, int printfl) -- int KINSetNumMaxIters(void *kinmem, long int mxiter) -- int KINSetNoInitSetup(void *kinmem, booleantype noInitSetup) -- int KINSetNoResMon(void *kinmem, booleantype noNNIResMon) -- int KINSetMaxSetupCalls(void *kinmem, long int msbset) -- int KINSetMaxSubSetupCalls(void *kinmem, long int msbsetsub) -- int KINSetEtaForm(void *kinmem, int etachoice) -- int KINSetEtaConstValue(void *kinmem, realtype eta) -- int KINSetEtaParams(void *kinmem, realtype egamma, realtype ealpha) -- int KINSetResMonParams(void *kinmem, realtype omegamin, realtype omegamax) -- int KINSetResMonConstValue(void *kinmem, realtype omegaconst) -- int KINSetNoMinEps(void *kinmem, booleantype noMinEps) -- int KINSetMaxNewtonStep(void *kinmem, realtype mxnewtstep) -- int KINSetMaxBetaFails(void *kinmem, long int mxnbcf) -- int KINSetRelErrFunc(void *kinmem, realtype relfunc) -- int KINSetFuncNormTol(void *kinmem, realtype fnormtol) -- int KINSetScaledStepTol(void *kinmem, realtype scsteptol) -- int KINSetConstraints(void *kinmem, N_Vector constraints) -- int KINSetSysFunc(void *kinmem, KINSysFn func) -- -- # solver routine -- int KINSol(void *kinmem, N_Vector uu, int strategy, N_Vector u_scale, N_Vector f_scale) -- -- # optional output routines. -- # Documentation see kinsol.h line 670-735 -- int KINGetWorkSpace(void *kinmem, long int *lenrw, long int *leniw) -- int KINGetNumNonlinSolvIters(void *kinmem, long int *nniters) -- int KINGetNumFuncEvals(void *kinmem, long int *nfevals) -- int KINGetNumBetaCondFails(void *kinmem, long int *nbcfails) -- int KINGetNumBacktrackOps(void *kinmem, long int *nbacktr) -- int KINGetFuncNorm(void *kinmem, realtype *fnorm) -- int KINGetStepLength(void *kinmem, realtype *steplength) -- char *KINGetReturnFlagName(int flag) -- -- # fuction used to deallocate memory used by KINSOL -- void KINFree(void **kinmem) -- --# linear solver modules --cdef extern from "kinsol/kinsol_dense.h": -- int KINDense(void *kinmem, int N) -- --cdef extern from "kinpinv.h": -- int KINPinv(void *kinmem, int N) -- --# functions used for supplying jacobian, and receiving info from linear solver --cdef extern from "kinsol/kinsol_direct.h": -- # user functions -- ctypedef int (*KINDlsDenseJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user functions to KINSOL -- int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac) -- -- # optional output fcts for linear direct solver -- int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINDlsGetLastFlag(void *kinmem, int *flag) -- char *KINDlsGetReturnFlagName(int flag) -- --cdef extern from "kinsol_jmod.h": -- # user functions -- ctypedef int (*KINPinvJacFn)(int N, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2) -- -- # function used to link user jacobian to KINSOL -- int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) -- -- # functions used to set regularizationparameter -- int KINPinvSetRegParam(void *kinmem, realtype reg_p) -- -- # optional output fcts for linear direct solver -- int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB) -- int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB) -- int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB) -- int KINPinvGetLastFlag(void *kinmem, int *flag) -- char *KINPinvGetReturnFlagName(int flag) -- -- --#========================= --# END SUNDIALS DEFINITIONS --#========================= -- --#=========================== --# Problem information struct --#=========================== --cdef class ProblemData: -- cdef: -- void *RHS # Residual or right-hand-side -- void *JAC -- int dim # Dimension of the problem -- --#================= --# Module functions --#================= -- --cdef inline N_Vector arr2nv(x): -- x=np.array(x) -- cdef long int n = len(x) -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- cdef N_Vector v=N_VNew_Serial(n) -- memcpy((<N_VectorContent_Serial>v.content).data, data_ptr, n*sizeof(double)) -- return v -- --cdef inline nv2arr(N_Vector v): -- cdef long int n = (<N_VectorContent_Serial>v.content).length -- cdef realtype* v_data = (<N_VectorContent_Serial>v.content).data -- cdef long int i -- import_array() -- cdef ndarray x=np.empty(n) -- #cdef ndarray x = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE,v_data) -- memcpy(x.data, v_data, n*sizeof(double)) -- return x -- --cdef inline realtype2arr(realtype *data, int n): -- """ Create new numpy array from realtype*""" -- import_array() -- cdef ndarray[realtype, ndim=1, mode='c'] x=np.empty(n) -- memcpy(x.data, data, n*sizeof(double)) -- return x -- --cdef inline arr2realtype(ndarray x, realtype* out_data, int n): -- """ Create new realtype* from numpy """ -- cdef ndarray[double, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(double)) -- -- --cdef inline arr2int(ndarray x, int* out_data, int n): -- """ Create new int* from numpy array """ -- cdef ndarray[int, ndim=1,mode='c'] ndx=x -- import_array() -- cdef void* data_ptr=PyArray_DATA(ndx) -- memcpy(out_data, data_ptr, n*sizeof(int)) -- -Index: src/lib/reg_routines.c -=================================================================== ---- assimulo/lib/reg_routines.c (revision 838) -+++ assimulo/lib/reg_routines.c (nonexistent) -@@ -1,202 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * This program is distributed in the hope that it will be useful, -- * but WITHOUT ANY WARRANTY; without even the implied warranty of -- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -- * GNU Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- * Johan Ylikiiskilä - johan.ylikiiskila@gmail.com -- * -- */ -- --#include "reg_routines.h" --#include <math.h> -- --SuperMatrix* regSparseMatrix(SuperMatrix *jac, double h) --{ -- SuperMatrix *res = NULL; -- int i=0,j=0,k=0,eli=0,elj=0,rowi=0,rowj=0,nnz=0,jnnz=0,jn=0,is_zero=0,max_nz=0; -- NCformat *Jstore = NULL; -- double *nzt = NULL, *jnzv=NULL,*nzval=NULL; -- double reg_param = h*h; -- int *rit = NULL, *colptr=NULL, *jri=NULL, *jcp=NULL,*rowind=NULL; -- Stype_t S = SLU_NC; -- Dtype_t D = SLU_D; -- Mtype_t M = SLU_GE; -- -- double *tmp_nzv = NULL; -- int *tmp_ri = NULL; -- -- -- /* Extract data */ -- Jstore = jac->Store; -- jn = jac->nrow; -- jnnz = Jstore->nnz; -- jnzv = Jstore->nzval; -- jri = Jstore->rowind; -- jcp = Jstore->colptr; -- -- /* Initialize data and allocate memory for result */ -- max_nz = 2*jnnz; -- res = calloc(1,sizeof(SuperMatrix)); -- nnz = 0; -- nzt = calloc(max_nz,sizeof(double)); -- rit = calloc(max_nz,sizeof(int)); -- colptr = calloc(jn+1,sizeof(int)); -- -- /* Calculate nnz,nzval,rowind and colptr */ -- colptr[0] = 0; -- -- /* Iterate over the columns of both J^T and J */ -- for (i = 0;i < jn;i++) { -- for (j = 0;j < jn;j++) { -- /* Get first index of first elt in both columns */ -- rowi = jcp[i]; -- rowj = jcp[j]; -- -- if (i==j) { -- /* Diagonal elt, perform regularization */ -- is_zero = 0; -- rit[nnz] = j; -- nzt[nnz++] += reg_param; -- } else { -- is_zero = 1; -- } -- -- /* Chack for matching and perform multiplication */ -- while ((rowi < jcp[i+1])&&(rowj < jcp[j+1])) { -- /* Get row numbers of the two current elts */ -- eli = jri[rowi]; -- elj = jri[rowj]; -- if (eli == elj) { -- /* Match! Non-zero elt to be inserted! */ -- /* If it is the first time, add to the non_zero structure */ -- if (is_zero) { -- is_zero = 0; -- rit[nnz++] = j; -- } -- nzt[nnz-1] += jnzv[rowi]*jnzv[rowj]; -- rowi++; -- rowj++; -- } else if (eli < elj) { -- rowi++; -- } else { -- rowj++; -- } -- } -- } -- /* New column, add info to colptr */ -- colptr[i+1] = nnz; -- -- /* -- * Check if we need to expand the allocated memory -- * and expand if necessary -- */ -- -- if (nnz > (max_nz-jn)) { -- tmp_nzv = calloc(2*max_nz,sizeof(double)); -- tmp_ri = calloc(2*max_nz,sizeof(int)); -- -- for (k = 0; k < nnz; k++) { -- tmp_nzv[k] = nzt[k]; -- tmp_ri[k] = rit[k]; -- } -- -- free(nzt); -- free(rit); -- -- nzt = tmp_nzv; -- rit = tmp_ri; -- -- tmp_nzv = NULL; -- tmp_ri = NULL; -- -- max_nz = 2*max_nz; -- } -- } -- -- /* Reallocate to the needed size of storage */ -- nzval = realloc(nzt,nnz*sizeof(double)); -- rowind = realloc(rit,nnz*sizeof(int)); -- -- /* Now create the matrix */ -- dCreate_CompCol_Matrix(res, jn, jn, nnz, nzval, rowind, colptr, S, D, M); -- return res ; --} -- --SuperMatrix* getRegRHS( SuperMatrix *jac, SuperMatrix *B) { -- SuperMatrix *res; -- NCformat *Jstore; -- DNformat *Bstore; -- double *data, *jac_data, *rhs_data; -- int *jac_ri, *jac_cp; -- int n,i,j,k ; -- -- /* Extract data from Jacobian */ -- Jstore = jac->Store; -- jac_data = Jstore->nzval; -- jac_ri = Jstore->rowind; -- jac_cp = Jstore->colptr; -- -- /* Extract data from right hand side */ -- n = B->nrow; -- Bstore = B->Store; -- rhs_data = Bstore->nzval; -- -- /* Allocate data for result */ -- data = calloc(n,sizeof(double)); -- -- /* Calculate data for result */ -- for (i = 0; i< n; i++) { -- for (j = jac_cp[i]; j < jac_cp[i+1]; j++) { -- k = jac_ri[j] ; -- data[i] += jac_data[j]*rhs_data[k]; -- } -- } -- -- /* Create the Matrix in SuperMatrix format */ -- res = calloc(1,sizeof(SuperMatrix)); -- dCreate_Dense_Matrix(res,n,1,data,n,SLU_DN,SLU_D,SLU_GE); -- -- return res; --} -- --double getRegParam(SuperMatrix *jac, SuperMatrix *B) { -- double res = 0; -- NCformat *Jstore; -- DNformat *Bstore; -- double data, *jac_data, *rhs_data; -- int *jac_ri, *jac_cp; -- int n,i,j,k ; -- -- /* Extract data from Jacobian */ -- Jstore = jac->Store; -- jac_data = Jstore->nzval; -- jac_ri = Jstore->rowind; -- jac_cp = Jstore->colptr; -- -- /* Extract data from right hand side */ -- n = B->nrow; -- Bstore = B->Store; -- rhs_data = Bstore->nzval; -- -- /* Calculate data for result */ -- for (i = 0; i< n; i++) { -- data = 0; -- for (j = jac_cp[i]; j < jac_cp[i+1]; j++) { -- k = jac_ri[j] ; -- data += jac_data[j]*rhs_data[k]; -- } -- res += data*data; -- } -- -- return sqrt(res); --} -Index: src/lib/sundials_kinsol_core.pyx -=================================================================== ---- assimulo/lib/sundials_kinsol_core.pyx (revision 838) -+++ assimulo/lib/sundials_kinsol_core.pyx (nonexistent) -@@ -1,371 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# This program is distributed in the hope that it will be useful, --# but WITHOUT ANY WARRANTY; without even the implied warranty of --# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --# GNU Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" -- --from __future__ import division -- --include "sundials_kinsol_includes.pxi" # Includes fcts from other header files --include "sundials_kinsol_core.pxi" # Includes the constants (textual include) -- --cdef int kin_res(N_Vector xv, N_Vector fval, void *problem_data): -- """ -- Residual fct called by KINSOL -- """ -- -- cdef: -- int i -- ProblemData pData = <ProblemData>problem_data -- ndarray[realtype, ndim = 1, mode = 'c'] rhs # Container holding return from user fct -- -- ndarray x = nv2arr(xv) -- realtype* resptr = (<N_VectorContent_Serial>fval.content).data -- -- try: -- rhs = (<object>pData.RHS)(x) -- for i in range(pData.dim): -- -- resptr[i] = rhs[i] -- -- return KIN_SUCCESS -- except: -- return KIN_REC_ERR -- --cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian, -- void *problem_data, N_Vector tmp1, N_Vector tmp2): -- """ -- This method is used to connect the Assimulo.Problem.jac to the Sundials -- Jacobian function. -- """ -- cdef: -- ProblemData pData = <ProblemData>problem_data -- realtype* col_i=DENSE_COL(Jacobian,0) -- ndarray x = nv2arr(xv) -- -- try: -- -- jac=(<object>pData.JAC)(x) -- #This needs further investigations: -- #memcpy(Jacobian.data,<realtype*>jac.data, pData.memSizeJac) -- -- for i in range(Neq): -- col_i = DENSE_COL(Jacobian, i) -- for j in range(Neq): -- col_i[j] = jac[j,i] -- -- return KINDLS_SUCCESS -- except: -- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description) -- -- --class KINError(Exception): -- """ -- Kinsol exception -- """ -- -- msg = {KIN_MEM_FAIL : 'Memory allocation failed.', -- KIN_MEM_NULL : 'KINSOL not properly created.', -- KIN_ILL_INPUT: 'There was an input error.', -- KIN_NO_MALLOC: 'Memory not allocated, call KINSOL_init(...)', -- KIN_LINESEARCH_NONCONV: 'Linesearch could not find a step big enough.', -- KIN_MAXITER_REACHED: 'Max number of iterations reached', -- KIN_MXNEWT_5X_EXCEEDED: 'Max size of newton step exceeded 5 times.', -- KIN_LINESEARCH_BCFAIL: 'Linesearch beta-test failed, probably because of too poor progress.', -- KIN_LINSOLV_NO_RECOVERY: 'psolve encountered an error, but preconditioner is current', -- KIN_LINIT_FAIL: 'Init of linear solver failed.', -- KIN_LSETUP_FAIL: 'pset (preconditioner setup fct) encountered an error', -- KIN_LSOLVE_FAIL: 'Error in either psolve or in linear solver routine.', -- KIN_SYSFUNC_FAIL: 'Call to RHS failed.', -- KIN_FIRST_SYSFUNC_ERR: 'Call to RHS failed on first call', -- KIN_REPTD_SYSFUNC_ERR: 'Call to RHS failed multiple times.', -- KIN_STEP_LT_STPTOL: 'Scaled step length too small. Either an approximate solution or a local minimum is reached. Check value of residual.', -- 234: 'KINSOL not compiled with SuperLU, sparse functionality not available.'} -- value = 0 -- def __init__(self, value): -- self.value = value -- -- def __str__(self): -- try: -- return repr(self.msg[self.value]) -- except KeyError: -- return repr('KINSOL failed with flag %s.'%(self.value)) -- --cdef class KINSOL_wrap: -- -- cdef: -- void* solver -- ProblemData pData # A struct containing problem data -- #ProblemData *ppData # Pointer to above mentioned struct -- N_Vector x_cur,x_init, x_scale, f_scale, con_nv -- booleantype noInitSetup, sparse -- realtype norm_of_res,reg_param -- ndarray con -- int print_level -- ndarray fscale -- -- def __cinit__(self): -- self.pData = ProblemData() #Create a new problem struct -- #self.ppData = &self.pData -- -- cpdef KINSOL_set_problem_info(self, RHS, dim, JAC= None): -- """ -- Sets the problem info to the desired struct -- """ -- -- # Sets residual function -- self.pData.RHS = <void*>RHS -- self.pData.dim = dim -- if JAC != None: -- self.pData.JAC = <void*>JAC -- else: -- self.pData.JAC = NULL -- -- def KINSOL_init(self,RHS,x0,dim, JAC = None, con = None,sparse = False, print_level = 0,norm_of_res = 0,reg_param=0,fscale=None): -- """ -- Initializes solver -- """ -- cdef int flag -- # Set problem info -- self.noInitSetup = False -- self.print_level = print_level -- self.norm_of_res = norm_of_res -- self.reg_param = reg_param -- self.KINSOL_set_problem_info(RHS,dim,JAC) -- self.con = con -- self.sparse = sparse -- self.fscale = fscale -- -- -- # Create initial guess from the supplied numpy array -- # print "x0 got in KINSOL wrapper: ", x0 -- self.x_cur = arr2nv(x0) -- self.x_init = arr2nv(x0) -- -- -- if self.solver == NULL: # solver runs for the first time -- -- # Create solver -- self.solver = KINCreate() -- if self.solver == NULL: -- raise KINError(KIN_MEM_FAIL) -- if self.print_level >0: -- print "KINSOL created" -- -- # Stop solver from performing precond setup -- flag = KINSetNoInitSetup(self.solver, self.noInitSetup) -- if flag < 0: -- raise KINError(flag) -- -- flag = KINSetPrintLevel(self.solver, self.print_level) -- if flag < 0: -- raise KINError(flag) -- -- if self.norm_of_res >0: -- flag = KINSetMaxNewtonStep(self.solver, self.norm_of_res) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetMaxNewtonStep(self.solver, 1e20) -- if flag < 0: -- raise KINError(flag) -- -- # If the user has specified constraints, connect them -- if self.con is not None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(self.con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- -- # Allocate internal memory -- flag = KINInit(self.solver,kin_res,self.x_cur) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "KINSOL initialized" -- -- # Link to linear solver -- if self.sparse: -- raise KINError(234) -- else: -- flag = KINPinv(self.solver,self.pData.dim) -- #flag = KINDense(self.solver,self.pData.dim) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Linear solver connected" -- -- # Set regularization parameter, if necessary -- if self.reg_param != 0.0: -- if self.sparse: -- raise KINError(234) -- else: -- flag = KINPinvSetRegParam(self.solver, self.reg_param) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Reg param ", self.reg_param, " set" -- else: -- # If the user has specified constraints, connect them -- if con is not None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetConstraints(self.solver, NULL) -- if flag < 0: -- raise KINError(flag) -- -- # If the user supplied a Jacobien, link it to the solver -- if self.pData.JAC != NULL: -- if sparse: -- raise KINError(234) -- else: -- flag = KINPinvSetJacFn(self.solver,kin_jac) -- #flag = KINDlsSetDenseJacFn(self.solver,kin_jac) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Jacobian supplied by user connected" -- else: -- if sparse: -- raise KINError(234) -- else: -- flag = KINPinvSetJacFn(self.solver,NULL) -- #flag = KINDlsSetDenseJacFn(self.solver,NULL) -- -- if flag < 0: -- raise KINError(flag) -- -- # Link the solver to the supplied problem data -- flag = KINSetUserData(self.solver,<void*>self.pData) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "User data set" -- -- # Create scaling vectors filled with ones -- # since the problem is assumed to be scaled -- self.x_scale = arr2nv(np.ones(self.pData.dim)) -- if self.fscale is not None: -- self.f_scale = arr2nv(self.fscale) -- else: -- self.f_scale = arr2nv(np.ones(self.pData.dim)) -- -- -- def KINSOL_solve(self,wo_LineSearch): -- """ -- Function that should be called after a call to KINSOL_init -- solves the function supplied as RHS -- """ -- -- cdef: -- int flag, DLSflag, lsflag, count,ind -- long int nb_iters -- ndarray x0 -- if self.print_level >0: -- print "Calling solver..." -- -- if not wo_LineSearch: -- # Use LineSearch -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_LINESEARCH,self.x_scale,self.f_scale) -- count = 0 -- if flag > 1: -- # Probably stuck at a minimum -- self.x_cur = self.x_init -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- -- if flag > 1: -- # Stuck at minimum again -- raise KINError(flag) -- else: -- # Don't use linesearch -- flag = -7 -- -- while flag <0: -- # Get nbr of non-linear iterations -- lsflag = KINGetNumNonlinSolvIters(<void*>self.solver, &nb_iters) -- if lsflag != 0: -- raise KINError(lsflag) -- -- -- if (flag == -5 and nb_iters == 1): -- # Cannot take a step bigger than zero on the first iteration -- # Try with heuristic initial guess -- raise KINError(42) -- -- if (flag == -7 or flag == -6 or (flag == -5 and nb_iters > 1)): -- if count == 0: -- # First time -- print "Calling solver again as simple newton." -- self.x_cur = self.x_init -- DLSflag = KINSetMaxSetupCalls(self.solver, 1); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 1: -- # second time -- print "Turning of residual monitoring" -- DLSflag = KINSetNoResMon(self.solver, True); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 3: -- # fourth time, release the bound on the step -- -- print "Remove limit on stepsize" -- DLSflag = KINSetMaxNewtonStep(self.solver, 1.0e34) -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- count += 1 -- if count >= 10: -- if self.print_level >0: -- print "Tried with simple newton ten times" -- raise KINError(-7) -- -- else: -- -- print "Error in solve, flag: ", flag -- lsflag = KINPinvGetLastFlag(self.solver, &lsflag) -- #lsflag = KINDlsGetLastFlag(self.solver, &lsflag) -- if lsflag != 0: -- if lsflag <0: -- print "Last flag from Linear solver: ", lsflag -- else: -- print "Jacobian singular at column ", lsflag -- raise KINError(flag) -- -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- if self.print_level >0: -- print "Problem solved, returning result" -- -- return nv2arr(self.x_cur) -Index: src/lib/kinsol_jmod.c -=================================================================== ---- assimulo/lib/kinsol_jmod.c (revision 838) -+++ assimulo/lib/kinsol_jmod.c (nonexistent) -@@ -1,481 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --/* -- * ================================================================= -- * IMPORTED HEADER FILES -- * ================================================================= -- */ -- --#include <stdio.h> --#include <stdlib.h> -- --#include "kinsol/kinsol_impl.h" --#include "kinsol_jmod_impl.h" --#include <sundials/sundials_math.h> -- --/* -- * ================================================================= -- * FUNCTION SPECIFIC CONSTANTS -- * ================================================================= -- */ -- --/* Constant for DQ Jacobian approximation */ --#define MIN_INC_MULT RCONST(1000.0) -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define uround (kin_mem->kin_uround) --#define func (kin_mem->kin_func) --#define user_data (kin_mem->kin_user_data) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#define sfdotJp (kin_mem->kin_sfdotJp) --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define bjac (kinpinv_mem->d_bjac) --#define J (kinpinv_mem->d_J) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define RTR (kinpinv_mem->d_RTR) --#define reg_param (kinpinv_mem->d_reg_param) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPinvSetJacFn -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (jac != NULL) { -- jacDQ = FALSE; -- djac = jac; -- } else { -- jacDQ = TRUE; -- } -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvSetRegParam -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetRegParam(void *kinmem, realtype reg_p) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- reg_param = reg_p; -- -- return(KINPINV_SUCCESS); --} -- -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetWorkSpace -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetWorkSpace(void *kinmem, long int *lenrwLS, long int *leniwLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (mtype == SUNDIALS_DENSE) { -- *lenrwLS = n*n; -- *leniwLS = n; -- } else if (mtype == SUNDIALS_BAND) { -- *lenrwLS = n*(smu + mu + 2*ml + 2); -- *leniwLS = n; -- } -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumJacEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumJacEvals(void *kinmem, long int *njevals) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *njevals = nje; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumFuncEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumFuncEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumGuncEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *nfevalsLS = nfeDQ; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetLastFlag -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetLastFlag(void *kinmem, int *flag) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *flag = last_flag; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetReturnFlagName -- * ----------------------------------------------------------------- -- */ -- --char *KINPinvGetReturnFlagName(int flag) --{ -- char *name; -- -- name = (char *)malloc(30*sizeof(char)); -- -- switch(flag) { -- case KINPINV_SUCCESS: -- sprintf(name, "KINPINV_SUCCESS"); -- break; -- case KINPINV_MEM_NULL: -- sprintf(name, "KINPINV_MEM_NULL"); -- break; -- case KINPINV_LMEM_NULL: -- sprintf(name, "KINPINV_LMEM_NULL"); -- break; -- case KINPINV_ILL_INPUT: -- sprintf(name, "KINPINV_ILL_INPUT"); -- break; -- case KINPINV_MEM_FAIL: -- sprintf(name, "KINPINV_MEM_FAIL"); -- break; -- default: -- sprintf(name, "NONE"); -- } -- -- return(name); --} -- -- --/* -- * ================================================================= -- * DQ JACOBIAN APPROXIMATIONS -- * ================================================================= -- */ -- -- -- -- --/* -- * ----------------------------------------------------------------- -- * kinPinvDQJac -- * ----------------------------------------------------------------- -- * This routine generates a dense difference quotient approximation to -- * the Jacobian of F(u). It assumes that a dense matrix of type -- * DlsMat is stored column-wise, and that elements within each column -- * are contiguous. The address of the jth column of J is obtained via -- * the macro DENSE_COL and this pointer is associated with an N_Vector -- * using the N_VGetArrayPointer/N_VSetArrayPointer functions. -- * Finally, the actual computation of the jth column of the Jacobian is -- * done with a call to N_VLinearSum. -- * -- * The increment used in the finite-difference approximation -- * J_ij = ( F_i(u+sigma_j * e_j) - F_i(u) ) / sigma_j -- * is -- * sigma_j = max{|u_j|, |1/uscale_j|} * sqrt(uround) -- * -- * Note: uscale_j = 1/typ(u_j) -- * -- * NOTE: Any type of failure of the system function her leads to an -- * unrecoverable failure of the Jacobian function and thus -- * of the linear solver setup function, stopping KINSOL. -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2) --{ -- realtype inc, inc_inv, ujsaved, ujscale, sign; -- realtype *tmp2_data, *u_data, *uscale_data; -- N_Vector ftemp, jthCol; -- long int j; -- int retval = 0; -- -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* data points to kin_mem */ -- kin_mem = (KINMem) data; -- kinpinv_mem = (KINPinvMem) lmem; -- -- /* Save pointer to the array in tmp2 */ -- tmp2_data = N_VGetArrayPointer(tmp2); -- -- /* Rename work vectors for readibility */ -- ftemp = tmp1; -- jthCol = tmp2; -- -- /* Obtain pointers to the data for u and uscale */ -- u_data = N_VGetArrayPointer(u); -- uscale_data = N_VGetArrayPointer(uscale); -- -- /* This is the only for loop for 0..N-1 in KINSOL */ -- -- for (j = 0; j < N; j++) { -- -- /* Generate the jth col of Jac(u) */ -- -- N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol); -- -- ujsaved = u_data[j]; -- ujscale = ONE/uscale_data[j]; -- sign = (ujsaved >= ZERO) ? ONE : -ONE; --#if SUNDIALS_26 -- inc = sqrt_relfunc*SUNMAX(SUNRabs(ujsaved), ujscale)*sign; --#else -- inc = sqrt_relfunc*MAX(ABS(ujsaved), ujscale)*sign; --#endif -- u_data[j] += inc; -- -- retval = func(u, ftemp, user_data); -- nfeDQ++; -- if(retval > 0) { -- /* try to recover by stepping in the opposite direction */ -- inc = -inc; -- u_data[j] = ujsaved + inc; -- -- retval = func(u, ftemp, user_data); -- nfeDQ++; -- } -- if (retval != 0) break; -- -- u_data[j] = ujsaved; -- -- inc_inv = ONE/inc; -- N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol); -- -- } -- -- /* Restore original array pointer in tmp2 */ -- N_VSetArrayPointer(tmp2_data, tmp2); -- -- return(retval); --} -Index: src/lib/kinsol_jmod_impl_wSLU.h -=================================================================== ---- assimulo/lib/kinsol_jmod_impl_wSLU.h (revision 838) -+++ assimulo/lib/kinsol_jmod_impl_wSLU.h (nonexistent) -@@ -1,168 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifndef _KINJMOD_IMPL_H --#define _KINJMOD_IMPL_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include "kinsol_jmod_wSLU.h" --#include "slu_ddefs.h" -- -- --/* -- * ----------------------------------------------------------------- -- * Types: KINPinvMemRec, KINPinvMem -- * ----------------------------------------------------------------- -- * The type KINPinvMem is pointer to a KINPinvMemRec. -- * This structure contains KINPinv/KINSLUG solver-specific data. -- * ----------------------------------------------------------------- -- */ -- --typedef struct KINPinvMemRec { -- -- int d_type; /* SUNDIALS_DENSE or SUNDIALS_BAND */ -- -- int d_n; /* problem dimension */ -- -- int d_ml; /* lower bandwidth of Jacobian */ -- int d_mu; /* upper bandwidth of Jacobian */ -- int d_smu; /* upper bandwith of M = MIN(N-1,d_mu+d_ml) */ -- -- booleantype d_jacDQ; /* TRUE if using internal DQ Jacobian approx. */ -- KINPinvJacFn d_djac; /* dense Jacobian routine to be called */ -- KINSLUGJacFn d_spjac; /* sparse Jacobian routine to be called */ -- -- void *d_J_data; /* J_data is passed to djac or bjac */ -- -- DlsMat d_J; /* problem Jacobian */ -- -- SuperMatrix *sp_J; /* sparse problem jacobian */ -- SuperMatrix *sp_L; /* L in the sparse LU factorization */ -- SuperMatrix *sp_U; /* U in the sparse LU factorization */ -- SuperMatrix *sp_B; /* sparse right hand side */ -- int *sp_perm_r; /* row permutations from partial pivoting */ -- int *sp_perm_c; /* column permutation vector */ -- -- SuperMatrix *sp_JTJ; /* Matrix needed for regularisation */ -- -- superlu_options_t *sp_options; /* options struct for SuperLU */ -- SuperLUStat_t *sp_stat; /* statistcis struct for SuperLU */ -- -- int *d_pivots; /* pivot array for PM = LU */ -- realtype *d_beta; -- realtype d_reg_param; /* Regularization parameter */ -- -- long int d_nje; /* no. of calls to jac */ -- -- long int d_nfeDQ; /* no. of calls to F due to DQ Jacobian approx. */ -- -- int d_last_flag; /* last error return flag */ -- -- DlsMat d_JTJ; -- booleantype d_regularized; /* Boolean set to true if problem is regularized*/ -- booleantype d_redojac; -- --} *KINPinvMem; -- -- --/* -- * ----------------------------------------------------------------- -- * Prototypes of internal functions -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2); -- --/* -- * ----------------------------------------------------------------- -- * Error Messages -- * ----------------------------------------------------------------- -- */ -- --#define MSGD_KINMEM_NULL "KINSOL memory is NULL." --#define MSGD_BAD_NVECTOR "A required vector operation is not implemented." --#define MSGD_MEM_FAIL "A memory request failed." --#define MSGD_LMEM_NULL "Linear solver memory is NULL." --#define MSGD_BAD_SIZES "Illegal bandwidth parameter(s). Must have 0 <= ml, mu <= N-1." --#define MSGD_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner." -- --#ifdef __cplusplus --} --#endif -- --#endif -Index: src/lib/reg_routines.h -=================================================================== ---- assimulo/lib/reg_routines.h (revision 838) -+++ assimulo/lib/reg_routines.h (nonexistent) -@@ -1,43 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * This program is distributed in the hope that it will be useful, -- * but WITHOUT ANY WARRANTY; without even the implied warranty of -- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -- * GNU Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- * Johan Ylikiiskilä - johan.ylikiiskila@gmail.com -- * -- */ -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#ifndef _REG_ROUTINES_H --#define _REG_ROUTINES_H -- -- --#include "slu_ddefs.h" -- -- /* -- * Routines used in regularisation -- */ -- --SuperMatrix* regSparseMatrix( SuperMatrix *jac, double h); -- --SuperMatrix* getRegRHS( SuperMatrix *jac, SuperMatrix *B); -- --double getRegParam(SuperMatrix *jac, SuperMatrix *B); --#endif -- --#ifdef __cplusplus --} --#endif -Index: src/lib/kinslug.c -=================================================================== ---- assimulo/lib/kinslug.c (revision 838) -+++ assimulo/lib/kinslug.c (nonexistent) -@@ -1,678 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * This program is distributed in the hope that it will be useful, -- * but WITHOUT ANY WARRANTY; without even the implied warranty of -- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -- * GNU Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- */ -- --/* -- * This is a linear solver using Tikhonov regularization and the sparse -- * solver SuperLU -- */ --#include <stdio.h> --#include <stdlib.h> --#include <math.h> -- --#include "slu_ddefs.h" --#include "supermatrix.h" --#include "reg_routines.h" --#include "kinsol_jmod_impl_wSLU.h" --#include "kinslug.h" --#include "kinsol/kinsol_impl.h" -- --#include <nvector/nvector_serial.h> --#include <sundials/sundials_direct.h> --#include <sundials/sundials_math.h> -- --/* Constants */ -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --#define DEBUG -- --/* -- * ================================================================= -- * PROTOTYPES FOR PRIVATE FUNCTIONS -- * ================================================================= -- */ -- -- --/* help functions */ --void DenSLU2Vector(N_Vector res, SuperMatrix *x); --void Vector2DenSLU(N_Vector x, SuperMatrix *res); -- -- --/* KINPinv linit, lsetup, lsolve, and lfree routines */ -- --static int kinSLUGInit(KINMem kin_mem); --static int kinSLUGSetup(KINMem kin_mem); --static int kinSLUGSolve(KINMem kin_mem, N_Vector x, N_Vector b, -- realtype *res_norm); --static void kinSLUGFree(KINMem kin_mem); -- -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define func (kin_mem->kin_func) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#define sfdotJp (kin_mem->kin_sfdotJp) --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define spjac (kinpinv_mem->d_spjac) --#define J (kinpinv_mem->d_J) --#define spJ (kinpinv_mem->sp_J) --#define L (kinpinv_mem->sp_L) --#define U (kinpinv_mem->sp_U) --#define B (kinpinv_mem->sp_B) --#define perm_r (kinpinv_mem->sp_perm_r) --#define perm_c (kinpinv_mem->sp_perm_c) --#define spJTJ (kinpinv_mem->sp_JTJ) --#define options (kinpinv_mem->sp_options) --#define stat (kinpinv_mem->sp_stat) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define JTJ (kinpinv_mem->d_JTJ) --#define regularized (kinpinv_mem->d_regularized) --#define redojac (kinpinv_mem->d_redojac) --#define beta (kinpinv_mem->d_beta) --#define reg_param (kinpinv_mem->d_reg_param) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINSLUG -- * ----------------------------------------------------------------- -- * This routine initializes the memory record and sets various function -- * fields specific to the sparse regularization solver module. -- * KINSLUG sets the kin_linit, kin_lsetup, kin_lsolve, kin_lfree fields -- * in *kinmem to be kinSLUGInit, kinSLUGSetup, kinSLUGSolve, and -- * kinSLUGFree, respectively. -- * It allocates memory for a structure of type KINPinvMemRec and sets -- * the kin_lmem field in *kinmem to the address of this structure. -- * It sets setupNonNull in *kinmem to TRUE, and the djac field to the -- * default kinPinvDenseDQJac. -- * Finally, it allocates memory for J and pivots. -- * -- * NOTE: The sparse regularized linear solver assumes a serial implementation -- * of the NVECTOR package. Therefore, KINSLUG will first -- * test for compatible a compatible N_Vector internal -- * representation by checking that N_VGetArrayPointer and -- * N_VSetArrayPointer exist. -- * ----------------------------------------------------------------- -- */ -- --int KINSLUG(void *kinmem, int N) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- /* Check if kinmem is different from NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinv", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- /* Test if the NVECTOR package is present */ -- if (vec_tmpl->ops->nvgetarraypointer == NULL || -- vec_tmpl->ops->nvsetarraypointer == NULL) { -- KINProcessError(kin_mem, KINPINV_ILL_INPUT, "KINPINV", "KINPinv", MSGD_BAD_NVECTOR); -- return(KINPINV_ILL_INPUT); -- } -- -- if (lfree !=NULL) lfree(kin_mem); -- -- /* Set four main function fields in kin_mem */ -- linit = kinSLUGInit; -- lsetup = kinSLUGSetup; -- lsolve = kinSLUGSolve; -- lfree = kinSLUGFree; -- -- /* Get memory for KINDlsMemRec */ -- kinpinv_mem = NULL; -- kinpinv_mem = (KINPinvMem) malloc(sizeof(struct KINPinvMemRec)); -- if (kinpinv_mem == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- return(KINPINV_MEM_FAIL); -- } -- -- /* Set matrix type */ -- mtype = SUNDIALS_DENSE; -- -- /* Set default Jacobian routine and Jacobian data */ -- jacDQ = TRUE; -- spjac = NULL; -- J_data = NULL; -- last_flag = KINPINV_SUCCESS; -- -- setupNonNull = TRUE; -- -- /* Set problem dimension */ -- n = N; -- -- /* Allocate memory for J,JTJ and pivots */ -- -- spJ = NULL; -- spJ = calloc(1,sizeof(SuperMatrix)) ; -- /*spJ = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );*/ -- if (spJ == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- L = NULL; -- L = calloc(1,sizeof(SuperMatrix)) ; -- if (L == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- U = NULL; -- U = calloc(1,sizeof(SuperMatrix)) ; -- if (U == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- B = NULL; -- B = calloc(1,sizeof(SuperMatrix)) ; -- if (B == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- perm_r = NULL; -- perm_r = intMalloc(N); -- if (perm_r == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- perm_c = NULL; -- perm_c = calloc(N,sizeof(int)); -- if (perm_c == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- options = NULL; -- options = (superlu_options_t *) calloc(1,sizeof(superlu_options_t)); -- if (options == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(perm_c); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- stat = NULL; -- stat = (SuperLUStat_t *) calloc(1,sizeof(SuperLUStat_t)); -- if (stat == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(perm_c); -- free(options); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- spJTJ = NULL; -- spJTJ= (SuperMatrix *) calloc(1,sizeof(SuperMatrix)); -- if (spJTJ == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINSLUG", "KINSLUG", MSGD_MEM_FAIL); -- free(spJ); -- free(L); -- free(U); -- free(B); -- free(perm_r); -- free(perm_c); -- free(options); -- free(stat); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- if (printfl>0) printf("Sparse data structures created \n"); -- -- /* Set the default input options. */ -- set_default_options(options); -- options->ColPerm = NATURAL; -- if (printfl>0) printf("Sparse options set \n"); -- -- /* Initialize the statistics variables. */ -- StatInit(stat); -- if (printfl>0) printf("Sparse stat initialized \n"); -- -- /* This is a direct linear solver */ -- inexact_ls = FALSE; -- -- /* Attach linear solver memory to integrator memory */ -- lmem = kinpinv_mem; -- -- /* Set reg_param 'not set' */ -- reg_param = 0; -- -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ================================================================= -- * PRIVATE FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGInit -- * ----------------------------------------------------------------- -- * This routine does remaining initializations specific to the Sparse -- * linear solver. -- * ----------------------------------------------------------------- -- */ -- --static int kinSLUGInit(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- nje = 0; -- nfeDQ = 0; -- -- /* -- * Set where to find jacobian data -- * if jacDQ=True it will be calculated by finite differenses -- * otherwise a user-supplied jacobian will be used -- */ -- -- -- if (jacDQ) { -- djac = kinPinvDQJac; -- J_data = kin_mem; -- } else { -- J_data = kin_mem->kin_user_data; -- } -- -- /* Set regularization parameter */ -- if (reg_param == 0){ -- reg_param = 1; -- } -- last_flag = KINPINV_SUCCESS; -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGSetup -- * ----------------------------------------------------------------- -- * This routine does the setup operations for the linear solver and -- * prepares J transpose J + h^2 I if necessary for regularization. -- * ----------------------------------------------------------------- -- */ -- --static int kinSLUGSetup(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- SuperMatrix spJC ; /* Matrix postmultiplied by P_c */ -- int retval,relax,panel_size,*etree,permc_spec; -- double rp = 0; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- regularized = FALSE; -- -- /* Calculate value of jacobian */ -- nje++; -- -- retval = spjac(n, uu, fval, spJ, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- last_flag = -1; -- return(-1); -- } -- -- /* Set flag that the matrix needs to be factorized */ -- -- options->Fact = DOFACT; -- -- -- /* Get column permutation vector perm_c */ -- permc_spec = options->ColPerm; -- get_perm_c(permc_spec,spJ,perm_c); -- -- /* Permute Jacobian and calculate and calculate etree */ -- etree = calloc(n,sizeof(int)); -- -- sp_preorder(options, spJ, perm_c, etree, &spJC); -- -- /* Get sizes for supernodes and panels from SUperLU */ -- panel_size = sp_ienv(1); -- relax = sp_ienv(2); -- -- -- /* Try to do a LU factorization of spJ */ -- dgstrf(options, &spJC, relax, panel_size, etree, -- NULL, 0,perm_c, perm_r, L, U, stat, &retval); -- if (retval != 0){ -- /* Factorization failed, determine reason */ -- if (retval < 0) { -- printf("Argument %i has illegal value\n",(-retval)); -- } else if (retval > n) { -- printf("%i nb of bytes allocated when memory allocation failed \n",(retval+n)); -- } else { -- /* Jacobian is singular, perform regularization */ -- if (printfl > 0) printf("Jacobian is singular at %i, regularizing \n ",retval); -- -- -- /* Recalculate jacobian */ -- -- Destroy_CompCol_Matrix(spJ); -- nje++; -- retval = spjac(n, uu, fval, spJ, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- last_flag = -1; -- free(etree); -- return(-1); -- } -- -- /* Calculate regularization parameter */ -- /* Put rhs in B2 */ -- B = calloc(1,sizeof(SuperMatrix)); -- Vector2DenSLU(fval,B); -- -- rp = getRegParam(spJ,B); -- -- printf("rp: %f \n",rp*rp); -- if (rp > 1) { -- reg_param = 1; -- } else { -- reg_param = rp; -- } -- if (printfl > 0) printf("Regparam: %f, \n ",reg_param); -- /* Calculate regularized matrix */ -- spJTJ = regSparseMatrix(spJ,sqrt(reg_param)); -- -- /* Get column permutation vector perm_c */ -- permc_spec = options->ColPerm; -- get_perm_c(permc_spec,spJTJ,perm_c); -- -- /* Do a preorder on the column etree */ -- sp_preorder(options, spJTJ, perm_c, etree, &spJC); -- -- /* Get sizes for supernodes and panels from SUperLU */ -- panel_size = sp_ienv(1); -- relax = sp_ienv(2); -- -- -- /* Try to do a LU factorization of spJ */ -- dgstrf(options, &spJC, relax, panel_size, etree, -- NULL, 0,perm_c, perm_r, L, U, stat, &retval); -- -- regularized = TRUE; -- if (retval == 0) { -- if (printfl > 0) printf("Jacobian regularized \n"); -- /* set flag that the matrix is factorized */ -- options->Fact = FACTORED; -- } else { -- if (printfl >0) { -- if (retval > n) { -- printf("Error in Jacobian regularization, %d nb of bytes allocated when allocation failed. n: %d",retval,n); -- } else { -- printf("Error in Jacobian regularization, retval: %d, n: %d,reg_param: %f",retval,n,reg_param); -- } -- } -- } -- } -- -- -- } else { -- /* set flag that the matrix is factorized */ -- options->Fact = FACTORED; -- fflush(stdout); -- } -- -- redojac = FALSE; -- /* Free allocated memory and return flags */ -- SUPERLU_FREE(etree); -- Destroy_CompCol_Permuted(&spJC); -- -- last_flag = retval; -- if (retval > 0) return(-1); -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * Vector2SDenSLU -- * ----------------------------------------------------------------- -- * Routine used to get the information from a NVector into the structure -- * needed by SuperLU. -- * ----------------------------------------------------------------- -- */ --void Vector2DenSLU(N_Vector x, SuperMatrix *res) { -- realtype *data; -- int length = NV_LENGTH_S(x); -- Stype_t S = SLU_DN ; -- Dtype_t D = SLU_D ; -- Mtype_t M = SLU_GE; -- -- /* Get data from the N_Vector */ -- data = N_VGetArrayPointer(x); -- -- /* Create dense Supermatrix */ -- dCreate_Dense_Matrix(res,length,1,data,length,S,D,M); -- --} -- --/* -- * ----------------------------------------------------------------- -- * DenSLU2Vector -- * ----------------------------------------------------------------- -- * Routine used to get the information from a SuperMatrix into a -- * N_Vector. -- * ----------------------------------------------------------------- -- */ --void DenSLU2Vector(N_Vector res, SuperMatrix *x) { -- realtype *data; -- DNformat *storage; -- -- /* Get data pointer from SuperMatrix */ -- storage = x->Store; -- data = storage->nzval; -- -- /* Set data pointer in N_Vector */ -- N_VSetArrayPointer(data,res); -- --} -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGSolve -- * ----------------------------------------------------------------- -- * This routine handles the solve operation for the sparse linear solver -- * by calling the sparse backsolve routine. The returned value is 0. -- * ----------------------------------------------------------------- -- */ -- --static int kinSLUGSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm) --{ -- KINPinvMem kinpinv_mem; -- SuperMatrix *B2; -- int info; -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (redojac) { -- return 1; -- } -- /* Copy the right-hand side into x*/ -- N_VScale(ONE, b, x); -- -- -- if (regularized) { -- -- /* Put rhs in B2 */ -- B2 = calloc(1,sizeof(SuperMatrix)); -- Vector2DenSLU(x,B2); -- -- /* Get regularized rhs */ -- B = getRegRHS(spJ,B2); -- -- regularized = FALSE; -- redojac = TRUE; -- -- /* Take out the garbage */ -- free(B2); -- } else { -- /* put rhs in B */ -- Vector2DenSLU(x,B); -- -- } -- -- /* Back-solve, get solution in B and propagate it to x*/ -- dgstrs (NOTRANS, L, U, perm_c, perm_r, B, stat, &info); -- -- if (info < 0) { -- printf("Solving failed with flag %i \n",info); -- last_flag = info; -- -- } else { -- DenSLU2Vector(x,B); -- } -- -- -- -- -- /* Compute the terms Jpnorm and sfdotJp for use in the global strategy -- routines and in KINForcingTerm. Both of these terms are subsequently -- corrected if the step is reduced by constraints or the line search. -- -- sJpnorm is the norm of the scaled product (scaled by fscale) of -- the current Jacobian matrix J and the step vector p. -- -- sfdotJp is the dot product of the scaled f vector and the scaled -- vector J*p, where the scaling uses fscale. */ -- -- sJpnorm = N_VWL2Norm(b,fscale); -- N_VProd(b, fscale, b); -- N_VProd(b, fscale, b); -- sfdotJp = N_VDotProd(fval, b); -- -- -- last_flag = KINPINV_SUCCESS; -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinSLUGFree -- * ----------------------------------------------------------------- -- * This routine frees memory specific to the dense linear solver. -- * ----------------------------------------------------------------- -- */ -- --static void kinSLUGFree(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- -- kinpinv_mem = (KINPinvMem) lmem; -- /* De-allocate storage */ -- -- SUPERLU_FREE (perm_r); -- SUPERLU_FREE (perm_c); -- Destroy_CompCol_Matrix(spJ); -- Destroy_Dense_Matrix(B); -- Destroy_SuperNode_Matrix(L); -- Destroy_CompCol_Matrix(U); -- Destroy_CompCol_Matrix(spJTJ); -- StatFree(stat); -- free(stat); -- free(options); -- free(spJ); -- free(B); -- free(L); -- free(U); -- free(spJTJ); -- -- free(kinpinv_mem); kinpinv_mem = NULL; --} -- -Index: src/lib/kinsol_jmod.h -=================================================================== ---- assimulo/lib/kinsol_jmod.h (revision 838) -+++ assimulo/lib/kinsol_jmod.h (nonexistent) -@@ -1,241 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- -- --#ifndef _KINJMOD_H --#define _KINJMOD_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include <sundials/sundials_direct.h> --#include <sundials/sundials_nvector.h> -- --/* -- * ================================================================= -- * K I N P I N V C O N S T A N T S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPINV return values -- * ----------------------------------------------------------------- -- */ -- --#define KINPINV_SUCCESS 0 --#define KINPINV_MEM_NULL -1 --#define KINPINV_LMEM_NULL -2 --#define KINPINV_ILL_INPUT -3 --#define KINPINV_MEM_FAIL -4 -- --/* Additional last_flag values */ -- --#define KINPINV_JACFUNC_UNRECVR -5 --#define KINPINV_JACFUNC_RECVR -6 -- --/* -- * ================================================================= -- * F U N C T I O N T Y P E S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Type: KINPinvJacFn -- * ----------------------------------------------------------------- -- * -- * A dense Jacobian approximation function Jac must be of type -- * KINDlsDenseJacFn. Its parameters are: -- * -- * N - problem size. -- * -- * u - current iterate (unscaled) [input] -- * -- * fu - vector (type N_Vector) containing result of nonlinear -- * system function evaluated at current iterate: -- * fu = F(u) [input] -- * -- * J - dense matrix (of type DlsMat) that will be loaded -- * by a KINDlsDenseJacFn with an approximation to the -- * Jacobian matrix J = (dF_i/dy_j). -- * -- * user_data - pointer to user data - the same as the user_data -- * parameter passed to KINSetFdata. -- * -- * tmp1, tmp2 - available scratch vectors (volatile storage) -- * -- * A KINPinvJacFn should return 0 if successful, a positive -- * value if a recoverable error occurred, and a negative value if -- * an unrecoverable error occurred. -- * -- * ----------------------------------------------------------------- -- * -- * NOTE: The following are two efficient ways to load a dense Jac: -- * (1) (with macros - no explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = DENSE_COL(Jac,j); -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * (2) (without macros - explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = (Jac->data)[j]; -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * A third way, using the DENSE_ELEM(A,i,j) macro, is much less -- * efficient in general. It is only appropriate for use in small -- * problems in which efficiency of access is NOT a major concern. -- * -- * ----------------------------------------------------------------- -- */ -- -- --typedef int (*KINPinvJacFn)(int N, -- N_Vector u, N_Vector fu, -- DlsMat J, void *user_data, -- N_Vector tmp1, N_Vector tmp2); -- -- --/* -- * ================================================================= -- * E X P O R T E D F U N C T I O N S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Optional inputs to the KINPinv linear solver -- * ----------------------------------------------------------------- -- * -- * KINDlsSetDenseJacFn specifies the dense Jacobian approximation -- * routine to be used for a direct dense linear solver. -- * -- * By default, a difference quotient approximation, supplied with -- * the solver is used. -- * -- * The return value is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- --SUNDIALS_EXPORT int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac); -- --/* -- Set regularization parameter -- */ --SUNDIALS_EXPORT int KINPinvSetRegParam(void *kinmem, realtype reg_p); -- --/* -- * ----------------------------------------------------------------- -- * Optional outputs from a KINDLS linear solver -- * ----------------------------------------------------------------- -- * -- * KINPinvGetWorkSpace returns the real and integer workspace used -- * by the KINDLS linear solver. -- * KINPinvGetNumJacEvals returns the number of calls made to the -- * Jacobian evaluation routine. -- * KINPinvGetNumFuncEvals returns the number of calls to the user's F -- * routine due to finite difference Jacobian -- * evaluation. -- * KINPinvGetLastFlag returns the last error flag set by any of -- * the KINDLS interface functions. -- * KINPinvGetReturnFlagName returns the name of the constant -- * associated with a KINDLS return flag -- * -- * The return value of KINPinvGet* is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- -- --SUNDIALS_EXPORT int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB); --SUNDIALS_EXPORT int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB); --SUNDIALS_EXPORT int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB); --SUNDIALS_EXPORT int KINPinvGetLastFlag(void *kinmem, int *flag); --SUNDIALS_EXPORT char *KINPinvGetReturnFlagName(int flag); -- --#ifdef __cplusplus --} --#endif -- --#endif -Index: src/lib/kinpinv.c -=================================================================== ---- assimulo/lib/kinpinv.c (revision 838) -+++ assimulo/lib/kinpinv.c (nonexistent) -@@ -1,578 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --/* -- *This is a linear solver using Tikhonov regularization -- */ --#include <stdio.h> --#include <stdlib.h> --#include <math.h> -- --#include "kinsol_jmod_impl.h" --#include "kinpinv.h" --#include "kinsol/kinsol_impl.h" -- --#include <sundials/sundials_direct.h> --#include <nvector/nvector_serial.h> --#include <sundials/sundials_math.h> -- --/* Constants */ -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --/* -- * ================================================================= -- * PROTOTYPES FOR PRIVATE FUNCTIONS -- * ================================================================= -- */ -- -- --/* help functions */ --void regMatrix(realtype **JTJ_c, realtype **jac, realtype h,int size); -- --/* KINPinv linit, lsetup, lsolve, and lfree routines */ -- --static int kinPinvInit(KINMem kin_mem); --static int kinPinvSetup(KINMem kin_mem); --static int kinPinvSolve(KINMem kin_mem, N_Vector x, N_Vector b, -- realtype *res_norm); --static void kinPinvFree(KINMem kin_mem); -- -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define func (kin_mem->kin_func) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#if SUNDIALS_26 -- #define sFdotJp (kin_mem->kin_sFdotJp) --#else -- #define sfdotJp (kin_mem->kin_sfdotJp) --#endif --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define J (kinpinv_mem->d_J) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define JTJ (kinpinv_mem->d_JTJ) --#define regularized (kinpinv_mem->d_regularized) -- --#define redojac (kinpinv_mem->d_redojac) --#define beta (kinpinv_mem->d_beta) --#define reg_param (kinpinv_mem->d_reg_param) -- --#define ihfun (kin_mem->kin_ihfun) --#define ih_data (kin_mem->kin_ih_data) --#define ehfun (kin_mem->kin_ehfun) --#define eh_data (kin_mem->kin_eh_data) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPinv -- * ----------------------------------------------------------------- -- * This routine initializes the memory record and sets various function -- * fields specific to the dense pseude-inverse linear solver module. -- * KINPinv sets the kin_linit, kin_lsetup, kin_lsolve, kin_lfree fields -- * in *kinmem to be kinPinvInit, kinPinvSetup, kinPinvSolve, and -- * kinPinvFree, respectively. -- * It allocates memory for a structure of type KINDlsMemRec and sets -- * the kin_lmem field in *kinmem to the address of this structure. -- * It sets setupNonNull in *kinmem to TRUE, and the djac field to the -- * default kinDlsDenseDQJac. -- * Finally, it allocates memory for J and pivots. -- * -- * NOTE: The dense pseudo-inverse linear solver assumes a serial implementation -- * of the NVECTOR package. Therefore, KINPinv will first -- * test for compatible a compatible N_Vector internal -- * representation by checking that N_VGetArrayPointer and -- * N_VSetArrayPointer exist. -- * ----------------------------------------------------------------- -- */ -- --int KINPinv(void *kinmem, int N) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- /* Check if kinmem is different from NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinv", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- /* Test if the NVECTOR package is present */ -- if (vec_tmpl->ops->nvgetarraypointer == NULL || -- vec_tmpl->ops->nvsetarraypointer == NULL) { -- KINProcessError(kin_mem, KINPINV_ILL_INPUT, "KINPINV", "KINPinv", MSGD_BAD_NVECTOR); -- return(KINPINV_ILL_INPUT); -- } -- -- if (lfree !=NULL) lfree(kin_mem); -- -- /* Set four main function fields in kin_mem */ -- linit = kinPinvInit; -- lsetup = kinPinvSetup; -- lsolve = kinPinvSolve; -- lfree = kinPinvFree; -- -- /* Get memory for KINDlsMemRec */ -- kinpinv_mem = NULL; -- kinpinv_mem = (KINPinvMem) malloc(sizeof(struct KINPinvMemRec)); -- if (kinpinv_mem == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- return(KINPINV_MEM_FAIL); -- } -- -- /* Set matrix type */ -- mtype = SUNDIALS_DENSE; -- -- /* Set default Jacobian routine and Jacobian data */ -- jacDQ = TRUE; -- djac = NULL; -- J_data = NULL; -- last_flag = KINPINV_SUCCESS; -- -- setupNonNull = TRUE; -- -- /* Set problem dimension */ -- n = N; -- -- /* Allocate memory for J,RTR and pivots */ -- -- J = NULL; -- J = NewDenseMat(N, N); -- if (J == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- JTJ = NULL; -- JTJ = NewDenseMat(n, n); -- if (JTJ == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- DestroyMat(J); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- pivots = NULL; -- pivots = NewIntArray(N); -- if (pivots == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- DestroyMat(J); -- DestroyMat(JTJ); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- beta = NULL; -- beta = NewRealArray(N); -- if (beta == NULL) { -- KINProcessError(kin_mem, KINPINV_MEM_FAIL, "KINPINV", "KINPinv", MSGD_MEM_FAIL); -- DestroyMat(J); -- DestroyMat(JTJ); -- DestroyArray(pivots); -- free(kinpinv_mem); kinpinv_mem = NULL; -- return(KINPINV_MEM_FAIL); -- } -- -- /* This is a direct linear solver */ -- inexact_ls = FALSE; -- -- /* Attach linear solver memory to integrator memory */ -- lmem = kinpinv_mem; -- /* Set reg_param 'not set' */ -- reg_param = 0; -- redojac=1 ; -- regularized = FALSE; -- nje = 0; -- nfeDQ = 0; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ================================================================= -- * PRIVATE FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * kinPinvInit -- * ----------------------------------------------------------------- -- * This routine does remaining initializations specific to the dense -- * linear solver. -- * ----------------------------------------------------------------- -- */ -- --static int kinPinvInit(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- nje = 0; -- nfeDQ = 0; -- -- /* -- * Set where to find jacobian data -- * if jacDQ=True it will be calculated by finite differenses -- * otherwise a user-supplied jacobian will be used -- */ -- -- if (jacDQ) { -- djac = kinPinvDQJac; -- J_data = kin_mem; -- } else { -- J_data = kin_mem->kin_user_data; -- } -- -- /* Set regularization parameter */ -- if (reg_param == 0){ -- reg_param = 1; -- } -- -- last_flag = KINPINV_SUCCESS; -- return(0); --} -- -- -- --/* -- * Function that calculates the regularized matrix for a given h -- * The result is returned in res while h is the regularization parameter -- * and J is the current prblem jacobian. -- */ -- --void regMatrix(realtype **JTJ_c, realtype **jac, realtype h, int size) --{ -- int i,j,k; -- -- -- i = 0; -- j = 0; -- k = 0; -- for (i=0;i<size;i++) { -- for (j=0;j<size;j++){ -- -- /*Calculate value at RTR(i,j) */ -- JTJ_c[j][i] = 0; -- for (k=0;k<size;k++) JTJ_c[j][i] += jac[j][k]*jac[i][k]; -- -- /* add the regularization parameter on the diagonal */ -- if (i==j)JTJ_c[j][i] += h*h; -- } -- } --} -- --/* -- * ----------------------------------------------------------------- -- * kinPinvSetup -- * ----------------------------------------------------------------- -- * This routine does the setup operations for the linear solver and -- * prepares J transpose J + h^2 I if necessary for regularization. -- * ----------------------------------------------------------------- -- */ -- --static int kinPinvSetup(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- long int ier; -- int retval; -- int i,j; -- -- realtype **JTJ_c ; -- realtype **jac ; -- realtype rp ; -- double data; -- realtype *b; -- -- char msg[200]; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- /* Calculate value of jacobian */ -- nje++; -- SetToZero(J); -- retval = djac(n, uu, fval, J, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- ehfun(KIN_LSETUP_FAIL,"KINPINV", "kinPinvSetup", "Jacobian evaluation function failed", eh_data); -- last_flag = -1; -- return(-1); -- } -- -- /* Try to do a LU factorization of J */ -- ier = DenseGETRF(J, (long int*)pivots); -- -- /* If the LU factorization failed, perform regularization */ -- if (ier > 0) { -- /* Calculate value of jacobian */ -- SetToZero(J); -- -- /* SetToZero(pivots);*/ -- for (i=0; i<n; i++) pivots[i] = 0; -- nje++; -- retval = djac(n, uu, fval, J, J_data, vtemp1, vtemp2); -- if (retval != 0) { -- ehfun(KIN_LSETUP_FAIL, "KINPINV", "kinPinvSetup", "Jacobian evaluation function failed", eh_data); -- last_flag = -1; -- return(-1); -- } -- -- /* Calculate J tranpose J */ -- SetToZero(JTJ); -- -- /* Calculate the regularization parameter */ -- -- jac = J->cols; -- b = N_VGetArrayPointer(fval); -- rp = 0; -- for (i = 0;i<n;i++) { -- data = 0; -- for(j = 0;j<n;j++) data += jac[i][j]*b[j]; -- rp += data*data ; -- } -- /*printf("rp: %f \n",rp);*/ -- if (sqrt(rp)<1) { -- reg_param = sqrt(rp); -- } else { -- reg_param = 1; -- } -- if (printfl > 0) { -- sprintf(msg, "Singular jacobian detected, regparam: %f",reg_param); -- ihfun("KINPINV", "kinPinvSetup", msg, ih_data); -- } -- /* calculate a regularized matrix */ -- JTJ_c = JTJ->cols; -- -- regMatrix(JTJ_c,jac,reg_param,n); -- -- /* LU-factorize the regularized matrix*/ -- ier = DenseGETRF(JTJ,(long int*)pivots); -- regularized = TRUE; -- } -- else { -- regularized = FALSE; -- } -- -- redojac = FALSE; -- /* Return 0 if the LU was complete; otherwise return -1 */ -- last_flag = ier; -- if (ier > 0) return(-1); -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinPinvSolve -- * ----------------------------------------------------------------- -- * This routine handles the solve operation for the dense linear solver -- * by calling the dense backsolve routine. The returned value is 0. -- * ----------------------------------------------------------------- -- */ -- --static int kinPinvSolve(KINMem kin_mem, N_Vector x, N_Vector b, realtype *res_norm) --{ -- KINPinvMem kinpinv_mem; -- realtype **jac; -- realtype *xd; -- realtype *bx; -- int i,j; -- -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (redojac) { -- return(1); -- } -- -- if (regularized) { -- if (printfl > 0) -- ihfun("KINPINV", "kinPinvSetup", "Solving regularized problem", ih_data); -- /* Calculate new right hand side b = J transpose * b */ -- bx = N_VGetArrayPointer(b); -- xd = N_VGetArrayPointer(x); -- jac = J->cols; -- i=0; -- j=0; -- for (i=0;i<n;i++){ -- xd[i] = 0; -- for (j=0;j<n;j++) xd[i] += jac[i][j]*bx[j]; -- } -- /* Back-solve and get solution in x */ -- -- /*N_VScale(ONE,x,b);*/ -- -- -- DenseGETRS(JTJ, (long int*)pivots, xd); -- -- /* Calculate a fresh jacobian - not really needed */ -- redojac = TRUE; -- } else { -- -- /* Copy the right-hand side into x */ -- -- N_VScale(ONE, b, x); -- xd = N_VGetArrayPointer(x); -- -- /* Back-solve and get solution in x */ -- DenseGETRS(J, (long int*)pivots, xd); -- -- } -- -- /* Compute the terms Jpnorm and sfdotJp for use in the global strategy -- routines and in KINForcingTerm. Both of these terms are subsequently -- corrected if the step is reduced by constraints or the line search. -- -- sJpnorm is the norm of the scaled product (scaled by fscale) of -- the current Jacobian matrix J and the step vector p. -- -- sfdotJp is the dot product of the scaled f vector and the scaled -- vector J*p, where the scaling uses fscale. */ -- -- -- sJpnorm = N_VWL2Norm(b,fscale); -- N_VProd(b, fscale, b); -- N_VProd(b, fscale, b); --#if SUNDIALS_26 -- sFdotJp = N_VDotProd(fval, b); --#else -- sfdotJp = N_VDotProd(fval, b); --#endif -- -- last_flag = KINPINV_SUCCESS; -- -- return(0); --} -- --/* -- * ----------------------------------------------------------------- -- * kinPinvFree -- * ----------------------------------------------------------------- -- * This routine frees memory specific to the dense linear solver. -- * ----------------------------------------------------------------- -- */ -- --static void kinPinvFree(KINMem kin_mem) --{ -- KINPinvMem kinpinv_mem; -- kinpinv_mem = (KINPinvMem) lmem; -- -- DestroyMat(J); -- DestroyMat(JTJ); -- DestroyArray(pivots); -- DestroyArray(beta); -- free(kinpinv_mem); kinpinv_mem = NULL; --} -- -Index: src/lib/sundials_kinsol_core_wSLU.pyx -=================================================================== ---- assimulo/lib/sundials_kinsol_core_wSLU.pyx (revision 838) -+++ assimulo/lib/sundials_kinsol_core_wSLU.pyx (nonexistent) -@@ -1,413 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# This program is distributed in the hope that it will be useful, --# but WITHOUT ANY WARRANTY; without even the implied warranty of --# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --# GNU Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --""" --from __future__ import division -- --import scipy.sparse as ss -- --include "sundials_kinsol_core_wSLU.pxi" # Includes fcts from other header files --include "sundials_kinsol_core.pxi" # Includes the constants (textual include) -- -- --cdef int kin_res(N_Vector xv, N_Vector fval, void *problem_data): -- """ -- Residual fct called by KINSOL -- """ -- -- cdef: -- int i -- ProblemData pData = <ProblemData>problem_data -- ndarray[realtype, ndim = 1, mode = 'c'] rhs # Container holding return from user fct -- -- ndarray x = nv2arr(xv) -- realtype* resptr = (<N_VectorContent_Serial>fval.content).data -- -- try: -- rhs = (<object>pData.RHS)(x) -- for i in range(pData.dim): -- -- resptr[i] = rhs[i] -- -- return KIN_SUCCESS -- except: -- return KIN_REC_ERR -- --cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian, -- void *problem_data, N_Vector tmp1, N_Vector tmp2): -- """ -- This method is used to connect the Assimulo.Problem.jac to the Sundials -- Jacobian function. -- """ -- cdef: -- ProblemData pData = <ProblemData>problem_data -- realtype* col_i=DENSE_COL(Jacobian,0) -- ndarray x = nv2arr(xv) -- -- try: -- -- jac=(<object>pData.JAC)(x) -- #This needs further investigations: -- #memcpy(Jacobian.data,<realtype*>jac.data, pData.memSizeJac) -- -- for i in range(Neq): -- col_i = DENSE_COL(Jacobian, i) -- for j in range(Neq): -- col_i[j] = jac[j,i] -- -- return KINDLS_SUCCESS -- except: -- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description) -- -- --cdef int kin_sp_jac(int Neq, N_Vector xv, N_Vector fval, SuperMatrix *Jacobian, -- void *problem_data, N_Vector tmp1, N_Vector tmp2): -- """ -- This method is used to connect the Assimulo.Problem.sparsejac to the KINSOL -- Jacobian function. -- """ -- cdef: -- ProblemData pData = <ProblemData>problem_data -- ndarray x = nv2arr(xv) -- -- int m,n,nnz,i -- realtype* nzval = NULL -- int* rowind = NULL -- int* colptr = NULL -- Stype_t S = SLU_NC -- Dtype_t D = SLU_D -- Mtype_t M = SLU_GE -- -- try: -- # Get jacobian and extract data -- jac=(<object>pData.JAC)(x) -- nnz = jac.nnz -- m = jac.shape[0] -- n = jac.shape[1] -- nzval = <realtype*>calloc(len(jac.data),sizeof(realtype)) -- rowind = <int*>calloc(len(jac.indices),sizeof(int)) -- colptr = <int*>calloc(len(jac.indptr),sizeof(int)) -- arr2realtype(jac.data,nzval,len(jac.data)) -- arr2int(jac.indices,rowind,len(jac.indices)) -- arr2int(jac.indptr,colptr,len(jac.indptr)) -- -- # Create the necessary matrix -- dCreate_CompCol_Matrix(Jacobian, m, n, nnz, nzval, rowind, colptr, S, D, M) -- -- -- return KINDLS_SUCCESS -- except: -- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description) -- -- --class KINError(Exception): -- """ -- Kinsol exception -- """ -- -- msg = {KIN_MEM_FAIL : 'Memory allocation failed.', -- KIN_MEM_NULL : 'KINSOL not properly created.', -- KIN_ILL_INPUT: 'There was an input error.', -- KIN_NO_MALLOC: 'Memory not allocated, call KINSOL_init(...)', -- KIN_LINESEARCH_NONCONV: 'Linesearch could not find a step big enough.', -- KIN_MAXITER_REACHED: 'Max number of iterations reached', -- KIN_MXNEWT_5X_EXCEEDED: 'Max size of newton step exceeded 5 times.', -- KIN_LINESEARCH_BCFAIL: 'Linesearch beta-test failed, probably because of too poor progress.', -- KIN_LINSOLV_NO_RECOVERY: 'psolve encountered an error, but preconditioner is current', -- KIN_LINIT_FAIL: 'Init of linear solver failed.', -- KIN_LSETUP_FAIL: 'pset (preconditioner setup fct) encountered an error', -- KIN_LSOLVE_FAIL: 'Error in either psolve or in linear solver routine.', -- KIN_SYSFUNC_FAIL: 'Call to RHS failed.', -- KIN_FIRST_SYSFUNC_ERR: 'Call to RHS failed on first call', -- KIN_REPTD_SYSFUNC_ERR: 'Call to RHS failed multiple times.', -- KIN_STEP_LT_STPTOL: 'Scaled step length too small. Either an approximate solution or a local minimum is reached. Check value of residual.'} -- value = 0 -- def __init__(self, value): -- self.value = value -- -- def __str__(self): -- try: -- return repr(self.msg[self.value]) -- except KeyError: -- return repr('KINSOL failed with flag %s.'%(self.value)) -- --cdef class KINSOL_wrap: -- -- cdef: -- void* solver -- ProblemData pData # A struct containing problem data -- #ProblemData *ppData # Pointer to above mentioned struct -- N_Vector x_cur,x_init, x_scale, f_scale, con_nv -- booleantype noInitSetup, sparse -- realtype norm_of_res,reg_param -- ndarray con -- int print_level -- ndarray fscale -- -- def __cinit__(self): -- self.pData = ProblemData() #Create a new problem struct -- #self.ppData = &self.pData -- -- cpdef KINSOL_set_problem_info(self, RHS, dim, JAC= None): -- """ -- Sets the problem info to the desired struct -- """ -- -- # Sets residual function -- self.pData.RHS = <void*>RHS -- self.pData.dim = dim -- if JAC != None: -- self.pData.JAC = <void*>JAC -- else: -- self.pData.JAC = NULL -- -- -- def KINSOL_init(self,RHS,x0,dim, JAC = None, con = None,sparse = False, print_level = 0,norm_of_res = 0,reg_param=0,fscale=None): -- """ -- Initializes solver -- """ -- cdef int flag -- # Set problem info -- self.noInitSetup = False -- self.print_level = print_level -- self.norm_of_res = norm_of_res -- self.reg_param = reg_param -- self.KINSOL_set_problem_info(RHS,dim,JAC) -- self.con = con -- self.sparse = sparse -- self.fscale = fscale -- -- -- # Create initial guess from the supplied numpy array -- # print "x0 got in KINSOL wrapper: ", x0 -- self.x_cur = arr2nv(x0) -- self.x_init = arr2nv(x0) -- -- -- if self.solver == NULL: # solver runs for the first time -- -- # Create solver -- self.solver = KINCreate() -- if self.solver == NULL: -- raise KINError(KIN_MEM_FAIL) -- if self.print_level >0: -- print "KINSOL created" -- -- # Stop solver from performing precond setup -- flag = KINSetNoInitSetup(self.solver, self.noInitSetup) -- if flag < 0: -- raise KINError(flag) -- -- flag = KINSetPrintLevel(self.solver, self.print_level) -- if flag < 0: -- raise KINError(flag) -- -- if self.norm_of_res >0: -- flag = KINSetMaxNewtonStep(self.solver, self.norm_of_res) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetMaxNewtonStep(self.solver, 1e20) -- if flag < 0: -- raise KINError(flag) -- -- # If the user has specified constraints, connect them -- if self.con != None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(self.con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- -- # Allocate internal memory -- flag = KINInit(self.solver,kin_res,self.x_cur) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "KINSOL initialized" -- -- # Link to linear solver -- if self.sparse: -- flag = KINSLUG(self.solver,self.pData.dim) -- else: -- flag = KINPinv(self.solver,self.pData.dim) -- #flag = KINDense(self.solver,self.pData.dim) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Linear solver connected" -- -- # Set regularization parameter, if necessary -- if self.reg_param != 0.0: -- if self.sparse: -- flag = KINSLUGSetRegParam(self.solver, self.reg_param) -- else: -- flag = KINPinvSetRegParam(self.solver, self.reg_param) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Reg param ", self.reg_param, " set" -- else: -- # If the user has specified constraints, connect them -- if con != None: -- if self.print_level >0: -- print "Applying constraints" -- self.con_nv = arr2nv(con) -- flag = KINSetConstraints(self.solver, self.con_nv) -- if flag < 0: -- raise KINError(flag) -- else: -- flag = KINSetConstraints(self.solver, NULL) -- if flag < 0: -- raise KINError(flag) -- -- # If the user supplied a Jacobien, link it to the solver -- if self.pData.JAC != NULL: -- if sparse: -- flag = KINSLUGSetJacFn(self.solver,kin_sp_jac) -- else: -- flag = KINPinvSetJacFn(self.solver,kin_jac) -- #flag = KINDlsSetDenseJacFn(self.solver,kin_jac) -- -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "Jacobian supplied by user connected" -- else: -- if sparse: -- flag = KINSLUGSetJacFn(self.solver,NULL) -- else: -- flag = KINPinvSetJacFn(self.solver,NULL) -- #flag = KINDlsSetDenseJacFn(self.solver,NULL) -- -- if flag < 0: -- raise KINError(flag) -- -- # Link the solver to the supplied problem data -- flag = KINSetUserData(self.solver,<void*>self.pData) -- if flag < 0: -- raise KINError(flag) -- if self.print_level >0: -- print "User data set" -- -- # Create scaling vectors filled with ones -- # since the problem is assumed to be scaled -- self.x_scale = arr2nv(np.ones(self.pData.dim)) -- if self.fscale != None: -- self.f_scale = arr2nv(self.fscale) -- else: -- self.f_scale = arr2nv(np.ones(self.pData.dim)) -- -- -- def KINSOL_solve(self,wo_LineSearch): -- """ -- Function that should be called after a call to KINSOL_init -- solves the function supplied as RHS -- """ -- -- cdef: -- int flag, DLSflag, lsflag, count,ind -- long int nb_iters -- ndarray x0 -- if self.print_level >0: -- print "Calling solver..." -- -- if not wo_LineSearch: -- # Use LineSearch -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_LINESEARCH,self.x_scale,self.f_scale) -- count = 0 -- if flag > 1: -- # Probably stuck at a minimum -- self.x_cur = self.x_init -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- -- if flag > 1: -- # Stuck at minimum again -- raise KINError(flag) -- else: -- # Don't use linesearch -- flag = -7 -- -- while flag <0: -- # Get nbr of non-linear iterations -- lsflag = KINGetNumNonlinSolvIters(<void*>self.solver, &nb_iters) -- if lsflag != 0: -- raise KINError(lsflag) -- -- -- if (flag == -5 and nb_iters == 1): -- # Cannot take a step bigger than zero on the first iteration -- # Try with heuristic initial guess -- raise KINError(42) -- -- if (flag == -7 or flag == -6 or (flag == -5 and nb_iters > 1)): -- if count == 0: -- # First time -- print "Calling solver again as simple newton." -- self.x_cur = self.x_init -- DLSflag = KINSetMaxSetupCalls(self.solver, 1); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 1: -- # second time -- print "Turning of residual monitoring" -- DLSflag = KINSetNoResMon(self.solver, True); -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- elif count == 3: -- # fourth time, release the bound on the step -- -- print "Remove limit on stepsize" -- DLSflag = KINSetMaxNewtonStep(self.solver, 1.0e34) -- if DLSflag < 0: -- raise KINError(DLSflag) -- -- count += 1 -- if count >= 10: -- if self.print_level >0: -- print "Tried with simple newton ten times" -- raise KINError(-7) -- -- else: -- -- print "Error in solve, flag: ", flag -- lsflag = KINPinvGetLastFlag(self.solver, &lsflag) -- #lsflag = KINDlsGetLastFlag(self.solver, &lsflag) -- if lsflag != 0: -- if lsflag <0: -- print "Last flag from Linear solver: ", lsflag -- else: -- print "Jacobian singular at column ", lsflag -- raise KINError(flag) -- -- flag = KINSol(<void*>self.solver,self.x_cur,KIN_NONE,self.x_scale,self.f_scale) -- if self.print_level >0: -- print "Problem solved, returning result" -- -- return nv2arr(self.x_cur) -Index: src/lib/kinslug.h -=================================================================== ---- assimulo/lib/kinslug.h (revision 838) -+++ assimulo/lib/kinslug.h (nonexistent) -@@ -1,58 +0,0 @@ --/* -- * Copyright (C) 2010 Modelon AB / Copyright (c) 2002, The Regents of the University of California. -- * -- * This program is free software: you can redistribute it and/or modify -- * it under the terms of the GNU Lesser General Public License as published by -- * the Free Software Foundation, version 3 of the License. -- -- * This program is distributed in the hope that it will be useful, -- * but WITHOUT ANY WARRANTY; without even the implied warranty of -- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -- * GNU Lesser General Public License for more details. -- -- * You should have received a copy of the GNU Lesser General Public License -- * along with this program. If not, see <http://www.gnu.org/licenses/>. -- * -- * This file is a modification of the file kinsol_dense.h , revision 1.5, -- * from the SUNDIALS suite. The file is modified by: -- * -- * Johan Ylikiiskilä - johan.ylikiiskila@gmail.com -- * -- */ -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#ifndef _KINSLU_H --#define _KINSLU_H -- --#include "kinsol_jmod_wSLU.h" --#include <sundials/sundials_dense.h> -- --/* -- * ----------------------------------------------------------------- -- * Function : KINSLUG -- * ----------------------------------------------------------------- -- * A call to the KINSLUG function links the main solver with the -- * SLUG (SuperLU reGularization implementation) linear solver. -- * Its arguments are as follows: -- * -- * kinmem - pointer to an internal memory block allocated during a -- * prior call to KINCreate -- * -- * N - problem size -- * -- * The return value of KINSLUG is one of: -- * 0 if successful -- * int different from zero otherwise -- * ----------------------------------------------------------------- -- */ -- -- SUNDIALS_EXPORT int KINSLUG(void *kinmem, int N); -- --#endif -- --#ifdef __cplusplus --} --#endif -Index: src/lib/kinsol_jmod_wSLU.c -=================================================================== ---- assimulo/lib/kinsol_jmod_wSLU.c (revision 838) -+++ assimulo/lib/kinsol_jmod_wSLU.c (nonexistent) -@@ -1,530 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --/* -- * ================================================================= -- * IMPORTED HEADER FILES -- * ================================================================= -- */ -- -- --#include <stdio.h> --#include <stdlib.h> -- --#include "kinsol/kinsol_impl.h" --#include "kinsol_jmod_impl_wSLU.h" --#include <sundials/sundials_math.h> -- --/* -- * ================================================================= -- * FUNCTION SPECIFIC CONSTANTS -- * ================================================================= -- */ -- --/* Constant for DQ Jacobian approximation */ --#define MIN_INC_MULT RCONST(1000.0) -- --#define ZERO RCONST(0.0) --#define ONE RCONST(1.0) --#define TWO RCONST(2.0) -- --/* -- * ================================================================= -- * READIBILITY REPLACEMENTS -- * ================================================================= -- */ -- --#define lrw1 (kin_mem->kin_lrw1) --#define liw1 (kin_mem->kin_liw1) --#define uround (kin_mem->kin_uround) --#define func (kin_mem->kin_func) --#define user_data (kin_mem->kin_user_data) --#define printfl (kin_mem->kin_printfl) --#define linit (kin_mem->kin_linit) --#define lsetup (kin_mem->kin_lsetup) --#define lsolve (kin_mem->kin_lsolve) --#define lfree (kin_mem->kin_lfree) --#define lmem (kin_mem->kin_lmem) --#define inexact_ls (kin_mem->kin_inexact_ls) --#define uu (kin_mem->kin_uu) --#define fval (kin_mem->kin_fval) --#define uscale (kin_mem->kin_uscale) --#define fscale (kin_mem->kin_fscale) --#define sqrt_relfunc (kin_mem->kin_sqrt_relfunc) --#define sJpnorm (kin_mem->kin_sJpnorm) --#define sfdotJp (kin_mem->kin_sfdotJp) --#define errfp (kin_mem->kin_errfp) --#define infofp (kin_mem->kin_infofp) --#define setupNonNull (kin_mem->kin_setupNonNull) --#define vtemp1 (kin_mem->kin_vtemp1) --#define vec_tmpl (kin_mem->kin_vtemp1) --#define vtemp2 (kin_mem->kin_vtemp2) -- --#define mtype (kinpinv_mem->d_type) --#define n (kinpinv_mem->d_n) --#define ml (kinpinv_mem->d_ml) --#define mu (kinpinv_mem->d_mu) --#define smu (kinpinv_mem->d_smu) --#define jacDQ (kinpinv_mem->d_jacDQ) --#define djac (kinpinv_mem->d_djac) --#define spjac (kinpinv_mem->d_spjac) --#define bjac (kinpinv_mem->d_bjac) --#define J (kinpinv_mem->d_J) --#define pivots (kinpinv_mem->d_pivots) --#define nje (kinpinv_mem->d_nje) --#define nfeDQ (kinpinv_mem->d_nfeDQ) --#define J_data (kinpinv_mem->d_J_data) --#define last_flag (kinpinv_mem->d_last_flag) --#define RTR (kinpinv_mem->d_RTR) --#define reg_param (kinpinv_mem->d_reg_param) -- --/* -- * ================================================================= -- * EXPORTED FUNCTIONS -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINDlsSetJacFn -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (jac != NULL) { -- jacDQ = FALSE; -- djac = jac; -- } else { -- jacDQ = TRUE; -- } -- -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvSetRegParam -- * ----------------------------------------------------------------- -- */ -- --int KINPinvSetRegParam(void *kinmem, realtype reg_p) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- reg_param = reg_p; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINSLUGSetRegParam -- * ----------------------------------------------------------------- -- */ -- --int KINSLUGSetRegParam(void *kinmem, realtype reg_p) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- reg_param = reg_p; -- -- return(KINPINV_SUCCESS); --} -- --int KINSLUGSetJacFn(void *kinmem, KINSLUGJacFn jac) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINSLUG", "KINSLUGSetJacFn", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINSLUG", "KINSLUGSetJacFn", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (jac != NULL) { -- jacDQ = FALSE; -- spjac = jac; -- } else { -- jacDQ = TRUE; -- } -- -- return(KINPINV_SUCCESS); --} --/* -- * ----------------------------------------------------------------- -- * KINDlsGetWorkSpace -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetWorkSpace(void *kinmem, long int *lenrwLS, long int *leniwLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetWorkSpace", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- if (mtype == SUNDIALS_DENSE) { -- *lenrwLS = n*n; -- *leniwLS = n; -- } else if (mtype == SUNDIALS_BAND) { -- *lenrwLS = n*(smu + mu + 2*ml + 2); -- *leniwLS = n; -- } -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumJacEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumJacEvals(void *kinmem, long int *njevals) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumJacEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *njevals = nje; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetNumFuncEvals -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsLS) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetNumFuncEvals", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetNumGuncEvals", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *nfevalsLS = nfeDQ; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetLastFlag -- * ----------------------------------------------------------------- -- */ -- --int KINPinvGetLastFlag(void *kinmem, int *flag) --{ -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* Return immediately if kinmem is NULL */ -- if (kinmem == NULL) { -- KINProcessError(NULL, KINPINV_MEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_KINMEM_NULL); -- return(KINPINV_MEM_NULL); -- } -- kin_mem = (KINMem) kinmem; -- -- if (lmem == NULL) { -- KINProcessError(kin_mem, KINPINV_LMEM_NULL, "KINPINV", "KINPinvGetLastFlag", MSGD_LMEM_NULL); -- return(KINPINV_LMEM_NULL); -- } -- kinpinv_mem = (KINPinvMem) lmem; -- -- *flag = last_flag; -- -- return(KINPINV_SUCCESS); --} -- --/* -- * ----------------------------------------------------------------- -- * KINPinvGetReturnFlagName -- * ----------------------------------------------------------------- -- */ -- --char *KINPinvGetReturnFlagName(int flag) --{ -- char *name; -- -- name = (char *)malloc(30*sizeof(char)); -- -- switch(flag) { -- case KINPINV_SUCCESS: -- sprintf(name, "KINPINV_SUCCESS"); -- break; -- case KINPINV_MEM_NULL: -- sprintf(name, "KINPINV_MEM_NULL"); -- break; -- case KINPINV_LMEM_NULL: -- sprintf(name, "KINPINV_LMEM_NULL"); -- break; -- case KINPINV_ILL_INPUT: -- sprintf(name, "KINPINV_ILL_INPUT"); -- break; -- case KINPINV_MEM_FAIL: -- sprintf(name, "KINPINV_MEM_FAIL"); -- break; -- default: -- sprintf(name, "NONE"); -- } -- -- return(name); --} -- --/* -- * ================================================================= -- * DQ JACOBIAN APPROXIMATIONS -- * ================================================================= -- */ -- -- -- -- --/* -- * ----------------------------------------------------------------- -- * kinPinvDQJac -- * ----------------------------------------------------------------- -- * This routine generates a dense difference quotient approximation to -- * the Jacobian of F(u). It assumes that a dense matrix of type -- * DlsMat is stored column-wise, and that elements within each column -- * are contiguous. The address of the jth column of J is obtained via -- * the macro DENSE_COL and this pointer is associated with an N_Vector -- * using the N_VGetArrayPointer/N_VSetArrayPointer functions. -- * Finally, the actual computation of the jth column of the Jacobian is -- * done with a call to N_VLinearSum. -- * -- * The increment used in the finite-difference approximation -- * J_ij = ( F_i(u+sigma_j * e_j) - F_i(u) ) / sigma_j -- * is -- * sigma_j = max{|u_j|, |1/uscale_j|} * sqrt(uround) -- * -- * Note: uscale_j = 1/typ(u_j) -- * -- * NOTE: Any type of failure of the system function her leads to an -- * unrecoverable failure of the Jacobian function and thus -- * of the linear solver setup function, stopping KINSOL. -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2) --{ -- realtype inc, inc_inv, ujsaved, ujscale, sign; -- realtype *tmp2_data, *u_data, *uscale_data; -- N_Vector ftemp, jthCol; -- long int j; -- int retval = 0; -- -- KINMem kin_mem; -- KINPinvMem kinpinv_mem; -- -- /* data points to kin_mem */ -- kin_mem = (KINMem) data; -- kinpinv_mem = (KINPinvMem) lmem; -- -- /* Save pointer to the array in tmp2 */ -- tmp2_data = N_VGetArrayPointer(tmp2); -- -- /* Rename work vectors for readibility */ -- ftemp = tmp1; -- jthCol = tmp2; -- -- /* Obtain pointers to the data for u and uscale */ -- u_data = N_VGetArrayPointer(u); -- uscale_data = N_VGetArrayPointer(uscale); -- -- /* This is the only for loop for 0..N-1 in KINSOL */ -- -- for (j = 0; j < N; j++) { -- -- /* Generate the jth col of Jac(u) */ -- -- N_VSetArrayPointer(DENSE_COL(Jac,j), jthCol); -- -- ujsaved = u_data[j]; -- ujscale = ONE/uscale_data[j]; -- sign = (ujsaved >= ZERO) ? ONE : -ONE; --#if SUNDIALS_26 -- inc = sqrt_relfunc*SUNMAX(SUNRabs(ujsaved), ujscale)*sign; --#else -- inc = sqrt_relfunc*MAX(ABS(ujsaved), ujscale)*sign; --#endif -- u_data[j] += inc; -- -- retval = func(u, ftemp, user_data); -- nfeDQ++; -- if (retval != 0) break; -- -- u_data[j] = ujsaved; -- -- inc_inv = ONE/inc; -- N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol); -- -- } -- -- /* Restore original array pointer in tmp2 */ -- N_VSetArrayPointer(tmp2_data, tmp2); -- -- return(retval); --} -Index: src/lib/kinpinv.h -=================================================================== ---- assimulo/lib/kinpinv.h (revision 838) -+++ assimulo/lib/kinpinv.h (nonexistent) -@@ -1,107 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#ifndef _PINV_H --#define _PINV_H -- --#include "kinsol_jmod.h" --#include <sundials/sundials_dense.h> -- --/* -- * ----------------------------------------------------------------- -- * Function : KINpinv -- * ----------------------------------------------------------------- -- * A call to the KINpinv function links the main solver with the -- * pinv linear solver. Its arguments are as follows: -- * -- * kinmem - pointer to an internal memory block allocated during a -- * prior call to KINCreate -- * -- * N - problem size -- * -- * The return value of KINpinv is one of: -- * 0 if successful -- * int different from zero otherwise -- * ----------------------------------------------------------------- -- */ -- -- SUNDIALS_EXPORT int KINPinv(void *kinmem, int N); -- --#endif -- --#ifdef __cplusplus --} --#endif -Index: src/lib/kinsol_jmod_wSLU.h -=================================================================== ---- assimulo/lib/kinsol_jmod_wSLU.h (revision 838) -+++ assimulo/lib/kinsol_jmod_wSLU.h (nonexistent) -@@ -1,255 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifndef _KINJMOD_H --#define _KINJMOD_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include <sundials/sundials_direct.h> --#include <sundials/sundials_nvector.h> -- --#include "slu_ddefs.h" -- --/* -- * ================================================================= -- * K I N P I N V C O N S T A N T S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * KINPINV return values -- * ----------------------------------------------------------------- -- */ -- --#define KINPINV_SUCCESS 0 --#define KINPINV_MEM_NULL -1 --#define KINPINV_LMEM_NULL -2 --#define KINPINV_ILL_INPUT -3 --#define KINPINV_MEM_FAIL -4 -- --/* Additional last_flag values */ -- --#define KINPINV_JACFUNC_UNRECVR -5 --#define KINPINV_JACFUNC_RECVR -6 -- --/* -- * ================================================================= -- * F U N C T I O N T Y P E S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Type: KINPinvJacFn -- * ----------------------------------------------------------------- -- * -- * A dense Jacobian approximation function Jac must be of type -- * KINDlsDenseJacFn. Its parameters are: -- * -- * N - problem size. -- * -- * u - current iterate (unscaled) [input] -- * -- * fu - vector (type N_Vector) containing result of nonlinear -- * system function evaluated at current iterate: -- * fu = F(u) [input] -- * -- * J - dense matrix (of type DlsMat) that will be loaded -- * by a KINDlsDenseJacFn with an approximation to the -- * Jacobian matrix J = (dF_i/dy_j). -- * -- * user_data - pointer to user data - the same as the user_data -- * parameter passed to KINSetFdata. -- * -- * tmp1, tmp2 - available scratch vectors (volatile storage) -- * -- * A KINPinvJacFn should return 0 if successful, a positive -- * value if a recoverable error occurred, and a negative value if -- * an unrecoverable error occurred. -- * -- * ----------------------------------------------------------------- -- * -- * NOTE: The following are two efficient ways to load a dense Jac: -- * (1) (with macros - no explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = DENSE_COL(Jac,j); -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * (2) (without macros - explicit data structure references) -- * for (j=0; j < Neq; j++) { -- * col_j = (Jac->data)[j]; -- * for (i=0; i < Neq; i++) { -- * generate J_ij = the (i,j)th Jacobian element -- * col_j[i] = J_ij; -- * } -- * } -- * A third way, using the DENSE_ELEM(A,i,j) macro, is much less -- * efficient in general. It is only appropriate for use in small -- * problems in which efficiency of access is NOT a major concern. -- * -- * ----------------------------------------------------------------- -- * ----------------------------------------------------------------- -- * Type: KINSLUGJacFn -- * ----------------------------------------------------------------- -- * -- * Basically the same as above except that instead of passing a -- * DlsMat to store the Jacobian SuperMatrix is passed. -- * -- * ----------------------------------------------------------------- -- */ -- -- --typedef int (*KINPinvJacFn)(int N, -- N_Vector u, N_Vector fu, -- DlsMat J, void *user_data, -- N_Vector tmp1, N_Vector tmp2); -- --typedef int (*KINSLUGJacFn)(int N, -- N_Vector u, N_Vector fu, -- SuperMatrix *J, void *user_data, -- N_Vector tmp1, N_Vector tmp2); -- -- --/* -- * ================================================================= -- * E X P O R T E D F U N C T I O N S -- * ================================================================= -- */ -- --/* -- * ----------------------------------------------------------------- -- * Optional inputs to the KINPinv linear solver -- * ----------------------------------------------------------------- -- * -- * KINDlsSetDenseJacFn specifies the dense Jacobian approximation -- * routine to be used for a direct dense linear solver. -- * -- * By default, a difference quotient approximation, supplied with -- * the solver is used. -- * -- * The return value is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- --SUNDIALS_EXPORT int KINPinvSetJacFn(void *kinmem, KINPinvJacFn jac); -- --SUNDIALS_EXPORT int KINSLUGSetJacFn(void *kinmem, KINSLUGJacFn jac); -- --SUNDIALS_EXPORT int KINSLUGSetRegParam(void *kinmem, realtype reg_p); -- --SUNDIALS_EXPORT int KINPinvSetRegParam(void *kinmem, realtype reg_p); -- --/* -- * ----------------------------------------------------------------- -- * Optional outputs from a KINDLS linear solver -- * ----------------------------------------------------------------- -- * -- * KINPinvGetWorkSpace returns the real and integer workspace used -- * by the KINDLS linear solver. -- * KINPinvGetNumJacEvals returns the number of calls made to the -- * Jacobian evaluation routine. -- * KINPinvGetNumFuncEvals returns the number of calls to the user's F -- * routine due to finite difference Jacobian -- * evaluation. -- * KINPinvGetLastFlag returns the last error flag set by any of -- * the KINDLS interface functions. -- * KINPinvGetReturnFlagName returns the name of the constant -- * associated with a KINDLS return flag -- * -- * The return value of KINPinvGet* is one of: -- * KINPINV_SUCCESS if successful -- * KINPINV_MEM_NULL if the KINSOL memory was NULL -- * KINPINV_LMEM_NULL if the linear solver memory was NULL -- * ----------------------------------------------------------------- -- */ -- --SUNDIALS_EXPORT int KINPinvGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB); --SUNDIALS_EXPORT int KINPinvGetNumJacEvals(void *kinmem, long int *njevalsB); --SUNDIALS_EXPORT int KINPinvGetNumFuncEvals(void *kinmem, long int *nfevalsB); --SUNDIALS_EXPORT int KINPinvGetLastFlag(void *kinmem, int *flag); --SUNDIALS_EXPORT char *KINPinvGetReturnFlagName(int flag); -- --#ifdef __cplusplus --} --#endif -- --#endif -Index: src/lib/sundials_kinsol_core.pxi -=================================================================== ---- assimulo/lib/sundials_kinsol_core.pxi (revision 838) -+++ assimulo/lib/sundials_kinsol_core.pxi (nonexistent) -@@ -1,144 +0,0 @@ --#!/usr/bin/env python --# -*- coding: utf-8 -*- -- --# Copyright (C) 2010 Modelon AB --# --# This program is free software: you can redistribute it and/or modify --# it under the terms of the GNU Lesser General Public License as published by --# the Free Software Foundation, version 3 of the License. --# --# This program is distributed in the hope that it will be useful, --# but WITHOUT ANY WARRANTY; without even the implied warranty of --# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the --# GNU Lesser General Public License for more details. --# --# You should have received a copy of the GNU Lesser General Public License --# along with this program. If not, see <http://www.gnu.org/licenses/>. -- --""" --Cython Wrapper for interfacing Python with KINSOL (Sundials Version 2.4.0) --Johan Ylikiiskilä Modelon AB -- --see also Jon Olav Vik: --http://codespeak.net/pipermail/cython-dev/2009-June/005947.html -- --This file contains the constants from Sundials. --""" --# ----------------------------------------------------------------- --# KINSOL return flags --# ----------------------------------------------------------------- -- --DEF KIN_SUCCESS =0 --DEF KIN_INITIAL_GUESS_OK =1 --DEF KIN_STEP_LT_STPTOL =2 -- --DEF KIN_REC_ERR =5 -- --DEF KIN_WARNING =99 -- --DEF KIN_MEM_NULL =-1 --DEF KIN_ILL_INPUT =-2 --DEF KIN_NO_MALLOC =-3 --DEF KIN_MEM_FAIL =-4 --DEF KIN_LINESEARCH_NONCONV =-5 --DEF KIN_MAXITER_REACHED =-6 --DEF KIN_MXNEWT_5X_EXCEEDED =-7 --DEF KIN_LINESEARCH_BCFAIL =-8 --DEF KIN_LINSOLV_NO_RECOVERY =-9 --DEF KIN_LINIT_FAIL =-10 --DEF KIN_LSETUP_FAIL =-11 --DEF KIN_LSOLVE_FAIL =-12 -- --DEF KIN_SYSFUNC_FAIL =-13 --DEF KIN_FIRST_SYSFUNC_ERR =-14 --DEF KIN_REPTD_SYSFUNC_ERR =-15 -- -- --# ----------------------------------------------------------------- --# Enumeration for inputs to KINSetEtaForm (eta choice) --# ----------------------------------------------------------------- --# KIN_ETACONSTANT : use constant value for eta (default value is --# 0.1 but a different value can be specified via --# a call to KINSetEtaConstValue) --# --# KIN_ETACHOICE1 : use choice #1 as given in Eisenstat and Walker's --# paper of SIAM J.Sci.Comput.,17 (1996), pp 16-32, --# wherein eta is defined to be: --# --# eta(k+1) = ABS(||F(u_k+1)||_L2-||F(u_k)+J(u_k)*p_k||_L2) --# --------------------------------------------- --# ||F(u_k)||_L2 --# --# 1+sqrt(5) --# eta_safe = eta(k)^ealpha where ealpha = --------- --# 2 --# --# KIN_ETACHOICE2 : use choice #2 as given in Eisenstat and Walker's --# paper wherein eta is defined to be: --# --# [ ||F(u_k+1)||_L2 ]^ealpha --# eta(k+1) = egamma * [ --------------- ] --# [ ||F(u_k)||_L2 ] --# --# where egamma = [0,1] and ealpha = (1,2] --# --# eta_safe = egamma*(eta(k)^ealpha) --# --# Note: The default values of the scalar --# coefficients egamma and ealpha (both required) --# are egamma = 0.9 and ealpha = 2.0, but the --# routine KINSetEtaParams can be used to specify --# different values. --# --# When using either KIN_ETACHOICE1 or KIN_ETACHOICE2, if --# eta_safe > 0.1 then the following safeguard is applied: --# --# eta(k+1) = MAX {eta(k+1), eta_safe} --# --# The following safeguards are always applied when using either --# KIN_ETACHOICE1 or KIN_ETACHOICE2 so that eta_min <= eta <= eta_max: --# --# eta(k+1) = MAX {eta(k+1), eta_min} --# eta(k+1) = MIN {eta(k+1), eta_max} --# --# where eta_min = 1.0e-4 and eta_max = 0.9 (see KINForcingTerm). --# ----------------------------------------------------------------- -- --DEF KIN_ETACHOICE1 =1 --DEF KIN_ETACHOICE2 =2 --DEF KIN_ETACONSTANT =3 -- --# ----------------------------------------------------------------- --# Enumeration for global strategy --# ----------------------------------------------------------------- --# Choices are KIN_NONE and KIN_LINESEARCH. --# ----------------------------------------------------------------- -- --DEF KIN_NONE =0 --DEF KIN_LINESEARCH =1 -- --# ----------------------------------------------------------------- --# KINDirect constants --# ----------------------------------------------------------------- -- --DEF KINDLS_SUCCESS =0 -- --DEF KINDLS_MEM_NULL =-1 --DEF KINDLS_LMEM_NULL =-2 --DEF KINDLS_ILL_INPUT =-3 --DEF KINDLS_MEM_FAIL =-4 --DEF KINDLS_JACFUNC_UNRECVR =-5 --DEF KINDLS_JACFUNC_RECVR =-6 -- -- --# ----------------------------------------------------------------- --# KINPinv constants --# ----------------------------------------------------------------- -- --DEF KINPINV_SUCCESS = 0 --DEF KINPINV_MEM_NULL = -1 --DEF KINPINV_LMEM_NULL = -2 --DEF KINPINV_ILL_INPUT = -3 --DEF KINPINV_MEM_FAIL = -4 --DEF KINPINV_JACFUNC_UNRECVR = -5 --DEF KINPINV_JACFUNC_RECVR = -6 -Index: src/lib/kinsol_jmod_impl.h -=================================================================== ---- assimulo/lib/kinsol_jmod_impl.h (revision 838) -+++ assimulo/lib/kinsol_jmod_impl.h (nonexistent) -@@ -1,153 +0,0 @@ --/* -- * This file has been modified by Modelon AB, changes are Copyright -- * (c) 2014 Modelon AB. -- * -- * Original copyright notice: -- * -- * Copyright (c) 2002, The Regents of the University of California. -- * Produced at the Lawrence Livermore National Laboratory. -- * Written by S.D. Cohen, A.C. Hindmarsh, R. Serban, -- * D. Shumaker, and A.G. Taylor. -- * UCRL-CODE-155951 (CVODE) -- * UCRL-CODE-155950 (CVODES) -- * UCRL-CODE-155952 (IDA) -- * UCRL-CODE-237203 (IDAS) -- * UCRL-CODE-155953 (KINSOL) -- * All rights reserved. -- -- * This file is part of SUNDIALS. -- -- * Redistribution and use in source and binary forms, with or without -- * modification, are permitted provided that the following conditions -- * are met: -- * -- * 1. Redistributions of source code must retain the above copyright -- * notice, this list of conditions and the disclaimer below. -- * -- * 2. Redistributions in binary form must reproduce the above copyright -- * notice, this list of conditions and the disclaimer (as noted below) -- * in the documentation and/or other materials provided with the -- * distribution. -- * -- * 3. Neither the name of the UC/LLNL nor the names of its contributors -- * may be used to endorse or promote products derived from this software -- * without specific prior written permission. -- * -- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -- * REGENTS OF THE UNIVERSITY OF CALIFORNIA, THE U.S. DEPARTMENT OF ENERGY -- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -- * -- * Additional BSD Notice -- * --------------------- -- * 1. This notice is required to be provided under our contract with -- * the U.S. Department of Energy (DOE). This work was produced at the -- * University of California, Lawrence Livermore National Laboratory -- * under Contract No. W-7405-ENG-48 with the DOE. -- * -- * 2. Neither the United States Government nor the University of -- * California nor any of their employees, makes any warranty, express -- * or implied, or assumes any liability or responsibility for the -- * accuracy, completeness, or usefulness of any information, apparatus, -- * product, or process disclosed, or represents that its use would not -- * infringe privately-owned rights. -- * -- * 3. Also, reference herein to any specific commercial products, -- * process, or services by trade name, trademark, manufacturer or -- * otherwise does not necessarily constitute or imply its endorsement, -- * recommendation, or favoring by the United States Government or the -- * University of California. The views and opinions of authors expressed -- * herein do not necessarily state or reflect those of the United States -- * Government or the University of California, and shall not be used for -- * advertising or product endorsement purposes. -- */ -- --#ifndef _KINJMOD_IMPL_H --#define _KINJMOD_IMPL_H -- --#ifdef __cplusplus /* wrapper to enable C++ usage */ --extern "C" { --#endif -- --#include "kinsol_jmod.h" -- -- --/* -- * ----------------------------------------------------------------- -- * Types: KINPinvMemRec, KINPinvMem -- * ----------------------------------------------------------------- -- * The type KINPinvMem is pointer to a KINPinvMemRec. -- * This structure contains KINPinv solver-specific data. -- * ----------------------------------------------------------------- -- */ -- --typedef struct KINPinvMemRec { -- -- int d_type; /* SUNDIALS_DENSE or SUNDIALS_BAND */ -- -- int d_n; /* problem dimension */ -- -- int d_ml; /* lower bandwidth of Jacobian */ -- int d_mu; /* upper bandwidth of Jacobian */ -- int d_smu; /* upper bandwith of M = MIN(N-1,d_mu+d_ml) */ -- -- booleantype d_jacDQ; /* TRUE if using internal DQ Jacobian approx. */ -- KINPinvJacFn d_djac; /* dense Jacobian routine to be called */ -- -- void *d_J_data; /* J_data is passed to djac or bjac */ -- -- DlsMat d_J; /* problem Jacobian */ -- -- int *d_pivots; /* pivot array for PM = LU */ -- realtype *d_beta; -- realtype d_reg_param; /* Regularization parameter */ -- long int d_nje; /* no. of calls to jac */ -- -- long int d_nfeDQ; /* no. of calls to F due to DQ Jacobian approx. */ -- -- int d_last_flag; /* last error return flag */ -- -- DlsMat d_JTJ; -- booleantype d_regularized; /* Flag if the regularized problem should be solved*/ -- booleantype d_redojac; -- --} *KINPinvMem; -- -- --/* -- * ----------------------------------------------------------------- -- * Prototypes of internal functions -- * ----------------------------------------------------------------- -- */ -- --int kinPinvDQJac(int N, -- N_Vector u, N_Vector fu, -- DlsMat Jac, void *data, -- N_Vector tmp1, N_Vector tmp2); -- --/* -- * ----------------------------------------------------------------- -- * Error Messages -- * ----------------------------------------------------------------- -- */ -- --#define MSGD_KINMEM_NULL "KINSOL memory is NULL." --#define MSGD_BAD_NVECTOR "A required vector operation is not implemented." --#define MSGD_MEM_FAIL "A memory request failed." --#define MSGD_LMEM_NULL "Linear solver memory is NULL." --#define MSGD_BAD_SIZES "Illegal bandwidth parameter(s). Must have 0 <= ml, mu <= N-1." --#define MSGD_JACFUNC_FAILED "The Jacobian routine failed in an unrecoverable manner." -- --#ifdef __cplusplus --} --#endif -- --#endif diff --git a/r840.patch b/r840.patch deleted file mode 100644 index eede3cd2b80e..000000000000 --- a/r840.patch +++ /dev/null @@ -1,206 +0,0 @@ -Index: setup.py -=================================================================== ---- setup.py (revision 839) -+++ setup.py (revision 840) -@@ -24,10 +24,17 @@ - import argparse - import numpy.distutils as nd - from distutils.version import StrictVersion -+from os import listdir -+from os.path import isfile, join - - def str2bool(v): -- return v.lower() in ("yes", "true", "t", "1") -+ return v.lower() in ("yes", "true", "t", "1") - -+def remove_prefix(name, prefix): -+ if name.startswith(prefix): -+ return name[len(prefix):] -+ return name -+ - parser = argparse.ArgumentParser(description='Assimulo setup script.') - parser.register('type','bool',str2bool) - package_arguments=['plugins','sundials','blas','superlu','lapack'] -@@ -36,7 +43,8 @@ - parser.add_argument("--{}-home".format(pg), - help="Location of the {} directory".format(pg.upper()),type=str,default='') - parser.add_argument("--blas-name", help="name of the blas package",default='blas') --parser.add_argument("--extra-c-flags", help='Extra C-flags (a list enclosed in " ")',default='') -+parser.add_argument("--extra-c-flags", help='Extra C-flags (a list enclosed in " ")',default='') -+parser.add_argument("--with_openmp", type='bool', help="set to true if present",default=False) - parser.add_argument("--is_static", type='bool', help="set to true if present",default=False) - parser.add_argument("--sundials-with-superlu", type='bool', help="set to true if Sundials has been compiled with SuperLU",default=False) - parser.add_argument("--debug", type='bool', help="set to true if present",default=False) -@@ -120,7 +128,6 @@ - self.static = args[0].is_static - self.static_link_gcc = ["-static-libgcc"] if self.static else [] - self.static_link_gfortran = ["-static-libgfortran"] if self.static else [] -- self.debug_flag = args[0].debug - self.force_32bit = args[0].force_32bit - self.flag_32bit = ["-m32"] if self.force_32bit else [] - self.no_mvscr = args[0].no_msvcr -@@ -128,7 +135,7 @@ - self.extra_fortran_link_flags = args[0].extra_fortran_link_flags.split() - self.extra_fortran_link_files = args[0].extra_fortran_link_files.split() - self.thirdparty_methods = thirdparty_methods -- -+ self.with_openmp = args[0].with_openmp - - - if self.no_mvscr: -@@ -270,12 +277,10 @@ - """ - Check if SuperLU package installed - """ -- self.with_SLU = self.with_BLAS -- kinsol_msg='KINSOL will not be compiled with support for SuperLU' -+ self.with_SLU = True -+ sundials_msg='SUNDIALS will not be compiled with support for SuperLU' - -- if not self.with_BLAS: -- L.warning(kinsol_msg+' as BLAS is missing.') -- elif self.SLUdir != "": -+ if self.SLUdir != "": - self.SLUincdir = os.path.join(self.SLUdir,'SRC') - self.SLUlibdir = os.path.join(self.SLUdir,'lib') - if not os.path.exists(os.path.join(self.SLUincdir,'supermatrix.h')): -@@ -283,16 +288,18 @@ - L.warning("Could not find SuperLU, disabling support. View more information using --log=DEBUG") - L.debug("Could not find SuperLU at the given path {}.".format(self.SLUdir)) - L.debug("usage: --superlu-home path") -- L.debug(kinsol_msg+'.') -+ L.debug(sundials_msg+'.') - else: - L.debug("SuperLU found in {} and {}: ".format(self.SLUincdir, self.SLUlibdir)) -+ -+ self.superLUFiles = [remove_prefix(f.split(".")[0],"lib") for f in listdir(self.SLUlibdir) if isfile(join(self.SLUlibdir, f)) and f.endswith(".a")] -+ self.superLUFiles.sort(reverse=True) - else: - L.warning("No path to SuperLU supplied, disabling support. View more information using --log=DEBUG") -- L.debug("No path to SuperLU supplied, KINSOL will not be compiled with support for SUperLU.") -+ L.debug("No path to SuperLU supplied, SUNDIALS will not be compiled with support for SuperLU.") - L.debug("usage: --superlu-home=path") - L.debug("Note: the path required is to the folder where the folders 'SRC' and 'lib' are found.") - self.with_SLU = False -- L.debug(kinsol_msg+'.') - - def check_SUNDIALS(self): - """ -@@ -301,13 +308,22 @@ - if os.path.exists(os.path.join(os.path.join(self.incdirs,'cvodes'), 'cvodes.h')): - self.with_SUNDIALS=True - L.debug('SUNDIALS found.') -+ sundials_version = None - -- if os.path.exists(os.path.join(os.path.join(self.incdirs,'arkode'), 'arkode.h')): #This was added in 2.6 -- sundials_version = (2,6,0) -- L.debug('SUNDIALS 2.6 found.') -- else: -- sundials_version = (2,5,0) -- L.debug('SUNDIALS 2.5 found.') -+ try: -+ if os.path.exists(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')): -+ with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f: -+ for line in f: -+ if "SUNDIALS_PACKAGE_VERSION" in line: -+ sundials_version = tuple([int(f) for f in line.split()[-1][1:-1].split(".")]) -+ L.debug('SUNDIALS %d.%d found.'%(sundials_version[0], sundials_version[1])) -+ except Exception: -+ if os.path.exists(os.path.join(os.path.join(self.incdirs,'arkode'), 'arkode.h')): #This was added in 2.6 -+ sundials_version = (2,6,0) -+ L.debug('SUNDIALS 2.6 found.') -+ else: -+ sundials_version = (2,5,0) -+ L.debug('SUNDIALS 2.5 found.') - - self.SUNDIALS_version = sundials_version - -@@ -367,6 +383,10 @@ - ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs] - ext_list[-1].library_dirs = [self.libdirs] - ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"] -+ if self.with_SLU: #If SUNDIALS is compiled with support for SuperLU -+ ext_list[-1].include_dirs.append(self.SLUincdir) -+ ext_list[-1].library_dirs.append(self.SLUlibdir) -+ ext_list[-1].libraries.extend(self.superLUFiles) - - #Kinsol - ext_list += cythonize(["assimulo"+os.path.sep+"solvers"+os.path.sep+"kinsol.pyx"], -@@ -386,6 +406,9 @@ - el.extra_compile_args = ["-O2", "-fno-strict-aliasing"] - if self.platform == "mac": - el.extra_compile_args += ["-Wno-error=return-type"] -+ if self.with_openmp: -+ el.extra_link_args.append("-fopenmp") -+ el.extra_compile_args.append("-fopenmp") - el.extra_compile_args += self.flag_32bit + self.extra_c_flags - - for el in ext_list: -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;
-@@ -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]
-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
diff --git a/r841.patch b/r841.patch deleted file mode 100644 index d7aa11b3ea6f..000000000000 --- a/r841.patch +++ /dev/null @@ -1,13 +0,0 @@ -Index: setup.py -=================================================================== ---- setup.py (revision 840) -+++ setup.py (revision 841) -@@ -383,7 +383,7 @@ - ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs] - ext_list[-1].library_dirs = [self.libdirs] - ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"] -- if self.with_SLU: #If SUNDIALS is compiled with support for SuperLU -+ if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU - ext_list[-1].include_dirs.append(self.SLUincdir) - ext_list[-1].library_dirs.append(self.SLUlibdir) - ext_list[-1].libraries.extend(self.superLUFiles) diff --git a/r845.patch b/r845.patch deleted file mode 100644 index dec23fc92e53..000000000000 --- a/r845.patch +++ /dev/null @@ -1,1689 +0,0 @@ -Index: setup.py -=================================================================== ---- setup.py (revision 844) -+++ setup.py (revision 845) -@@ -309,15 +309,30 @@ - self.with_SUNDIALS=True - L.debug('SUNDIALS found.') - sundials_version = None -+ sundials_vector_type_size = None - - try: - if os.path.exists(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')): - with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f: - for line in f: -- if "SUNDIALS_PACKAGE_VERSION" in line: -+ if "SUNDIALS_PACKAGE_VERSION" in line or "SUNDIALS_VERSION" in line: - sundials_version = tuple([int(f) for f in line.split()[-1][1:-1].split(".")]) - L.debug('SUNDIALS %d.%d found.'%(sundials_version[0], sundials_version[1])) -- except Exception: -+ break -+ with open(os.path.join(os.path.join(self.incdirs,'sundials'), 'sundials_config.h')) as f: -+ for line in f: -+ if "SUNDIALS_INT32_T" in line and line.startswith("#define"): -+ sundials_vector_type_size = "32" -+ L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size)) -+ break -+ if "SUNDIALS_INT64_T" in line and line.startswith("#define"): -+ sundials_vector_type_size = "64" -+ L.debug('SUNDIALS vector type size %s bit found.'%(sundials_vector_type_size)) -+ if self.with_SLU: -+ L.warning("It is recommended to set the SUNDIALS_INDEX_TYPE to an 32bit integer when using SUNDIALS together with SuperLU.") -+ L.warning("SuperLU may not function properly.") -+ break -+ except Exception as e: - if os.path.exists(os.path.join(os.path.join(self.incdirs,'arkode'), 'arkode.h')): #This was added in 2.6 - sundials_version = (2,6,0) - L.debug('SUNDIALS 2.6 found.') -@@ -326,6 +341,7 @@ - L.debug('SUNDIALS 2.5 found.') - - self.SUNDIALS_version = sundials_version -+ self.SUNDIALS_vector_size = sundials_vector_type_size - - else: - L.warning(("Could not find Sundials, check the provided path (--sundials-home={}) "+ -@@ -375,7 +391,8 @@ - # SUNDIALS - if self.with_SUNDIALS: - compile_time_env = {'SUNDIALS_VERSION': self.SUNDIALS_version, -- 'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu} -+ 'SUNDIALS_WITH_SUPERLU': self.sundials_with_superlu, -+ 'SUNDIALS_VECTOR_SIZE': self.SUNDIALS_vector_size} - #CVode and IDA - ext_list += cythonize(["assimulo" + os.path.sep + "solvers" + os.path.sep + "sundials.pyx"], - include_path=[".","assimulo","assimulo" + os.sep + "lib"], -@@ -382,11 +399,19 @@ - compile_time_env=compile_time_env, force=True) - ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs] - ext_list[-1].library_dirs = [self.libdirs] -- ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"] -+ -+ if self.SUNDIALS_version >= (3,0,0): -+ ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas", "sundials_sunlinsoldense", "sundials_sunlinsolspgmr", "sundials_sunmatrixdense", "sundials_sunmatrixsparse"] -+ else: -+ ext_list[-1].libraries = ["sundials_cvodes", "sundials_nvecserial", "sundials_idas"] - if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU -+ if self.SUNDIALS_version >= (3,0,0): -+ ext_list[-1].libraries.extend(["sundials_sunlinsolsuperlumt"]) -+ - ext_list[-1].include_dirs.append(self.SLUincdir) - ext_list[-1].library_dirs.append(self.SLUlibdir) - ext_list[-1].libraries.extend(self.superLUFiles) -+ - - #Kinsol - ext_list += cythonize(["assimulo"+os.path.sep+"solvers"+os.path.sep+"kinsol.pyx"], -@@ -395,7 +420,11 @@ - ext_list[-1].include_dirs = [np.get_include(), "assimulo","assimulo"+os.sep+"lib", self.incdirs] - ext_list[-1].library_dirs = [self.libdirs] - ext_list[-1].libraries = ["sundials_kinsol", "sundials_nvecserial"] -- -+ -+ if self.sundials_with_superlu and self.with_SLU: #If SUNDIALS is compiled with support for SuperLU -+ ext_list[-1].include_dirs.append(self.SLUincdir) -+ ext_list[-1].library_dirs.append(self.SLUlibdir) -+ ext_list[-1].libraries.extend(self.superLUFiles) - - for el in ext_list: - #Debug -Index: assimulo/solvers/__init__.py -=================================================================== ---- assimulo/solvers/__init__.py (revision 844) -+++ assimulo/solvers/__init__.py (revision 845) -@@ -25,7 +25,7 @@ - from .euler import ExplicitEuler
- from .euler import ImplicitEuler
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .radau5 import Radau5ODE
- from .radau5 import Radau5DAE
-@@ -32,43 +32,43 @@ - from .radau5 import _Radau5ODE
- from .radau5 import _Radau5DAE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .sundials import IDA
- from .sundials import CVode
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .kinsol import KINSOL
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .runge_kutta import RungeKutta34
- from .runge_kutta import RungeKutta4
- from .runge_kutta import Dopri5
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .rosenbrock import RodasODE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .odassl import ODASSL
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .odepack import LSODAR
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .radar5 import Radar5ODE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .dasp3 import DASP3ODE
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
- try:
- from .glimda import GLIMDA
- except ImportError as ie:
-- sys.stderr.write("Could not find {}\n".format(ie.args[0].split("'")[1]))
-+ sys.stderr.write("Could not find " + str(ie) + "\n")
-Index: src/solvers/kinsol.pyx -=================================================================== ---- assimulo/solvers/kinsol.pyx (revision 844) -+++ assimulo/solvers/kinsol.pyx (revision 845) -@@ -30,8 +30,8 @@ - cimport sundials_includes as SUNDIALS
-
- #Various C includes transfered to namespace
--from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL
--from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat
-+from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype
-+from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse
- from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial
- from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial
-
-@@ -57,6 +57,8 @@ -
- cdef object pt_fcn, pt_jac, pt_jacv, pt_prec_setup, pt_prec_solve
- cdef object _added_linear_solver
-+ cdef SUNDIALS.SUNMatrix sun_matrix
-+ cdef SUNDIALS.SUNLinearSolver sun_linearsolver
-
- def __init__(self, problem):
- Algebraic.__init__(self, problem) #Calls the base class
-@@ -86,6 +88,7 @@ - self.options["no_min_epsilon"] = False #Specifies wheter the scaled linear residual is bounded from below
- self.options["max_beta_fails"] = 10
- self.options["max_krylov"] = 0
-+ self.options["precond"] = PREC_NONE
-
- #Statistics
- self.statistics["nfevals"] = 0 #Function evaluations
-@@ -109,6 +112,13 @@ - if self.kinsol_mem != NULL:
- #Free Memory
- SUNDIALS.KINFree(&self.kinsol_mem)
-+
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ if self.sun_matrix != NULL:
-+ SUNDIALS.SUNMatDestroy(self.sun_matrix)
-+
-+ if self.sun_linearsolver != NULL:
-+ SUNDIALS.SUNLinSolFree(self.sun_linearsolver)
-
- def update_variable_scaling(self, value="Automatic"):
- """
-@@ -170,7 +180,7 @@ -
- cdef initialize_kinsol(self):
- cdef int flag #Used for return
--
-+
- self.y_temp = arr2nv(self.y)
- self.y_scale = arr2nv([1.0]*self.problem_info["dim"])
- self.f_scale = arr2nv([1.0]*self.problem_info["dim"])
-@@ -208,16 +218,34 @@ -
- cpdef add_linear_solver(self):
- if self.options["linear_solver"] == "DENSE":
-- flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ #Create a dense Sundials matrix
-+ self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim)
-+ #Create a dense Sundials linear solver
-+ self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.y_temp, self.sun_matrix)
-+ #Attach it to Kinsol
-+ flag = SUNDIALS.KINDlsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver, self.sun_matrix)
-+ ELSE:
-+ flag = SUNDIALS.KINDense(self.kinsol_mem, self.problem_info["dim"])
- if flag < 0:
- raise KINSOLError(flag)
-
- if self.problem_info["jac_fcn"]:
-- flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ flag = SUNDIALS.KINDlsSetJacFn(self.kinsol_mem, kin_jac);
-+ ELSE:
-+ flag = SUNDIALS.KINDlsSetDenseJacFn(self.kinsol_mem, kin_jac);
- if flag < 0:
- raise KINSOLError(flag)
- elif self.options["linear_solver"] == "SPGMR":
-- flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
-+ IF SUNDIALS_VERSION >= (3,0,0):
-+ #Create the linear solver
-+ self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.y_temp, self.options["precond"], self.options["max_krylov"])
-+ #Attach it to Kinsol
-+ flag = SUNDIALS.KINSpilsSetLinearSolver(self.kinsol_mem, self.sun_linearsolver)
-+ ELSE:
-+ #Specify the use of KINSpgmr linear solver.
-+ flag = SUNDIALS.KINSpgmr(self.kinsol_mem, self.options["max_krylov"])
- if flag < 0:
- raise KINSOLError(flag)
-
-Index: src/solvers/sundials.pyx -=================================================================== ---- assimulo/solvers/sundials.pyx (revision 844) -+++ assimulo/solvers/sundials.pyx (revision 845) -@@ -34,8 +34,8 @@ - cimport sundials_includes as SUNDIALS - - #Various C includes transfered to namespace --from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL --from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat -+from sundials_includes cimport N_Vector, realtype, N_VectorContent_Serial, DENSE_COL, sunindextype -+from sundials_includes cimport memcpy, N_VNew_Serial, DlsMat, SlsMat, SUNMatrix, SUNMatrixContent_Dense, SUNMatrixContent_Sparse - from sundials_includes cimport malloc, free, realtype, N_VCloneVectorArray_Serial - from sundials_includes cimport N_VConst_Serial, N_VDestroy_Serial - -@@ -70,6 +70,8 @@ - cdef public N.ndarray yS0 - #cdef N.ndarray _event_info - cdef public N.ndarray g_old -+ cdef SUNDIALS.SUNMatrix sun_matrix -+ cdef SUNDIALS.SUNLinearSolver sun_linearsolver - - def __init__(self, problem): - Implicit_ODE.__init__(self, problem) #Calls the base class -@@ -100,7 +102,9 @@ - self.options["dqrhomax"] = 0.0 - self.options["pbar"] = [1]*self.problem_info["dimSens"] - self.options["external_event_detection"] = False #Sundials rootfinding is used for event location as default -+ self.options["precond"] = PREC_NONE - -+ - #Solver support - self.supports["report_continuously"] = True - self.supports["interpolated_output"] = True -@@ -173,6 +177,13 @@ - if self.ida_mem != NULL: - #Free Memory - SUNDIALS.IDAFree(&self.ida_mem) -+ -+ IF SUNDIALS_VERSION >= (3,0,0): -+ if self.sun_matrix != NULL: -+ SUNDIALS.SUNMatDestroy(self.sun_matrix) -+ -+ if self.sun_linearsolver != NULL: -+ SUNDIALS.SUNLinSolFree(self.sun_linearsolver) - - cpdef state_event_info(self): - """ -@@ -255,21 +266,30 @@ - if flag < 0: - raise IDAError(flag, self.t) - -- #Specify the use of the internal dense linear algebra functions. -- flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim) -- if flag < 0: -- raise IDAError(flag, self.t) -- -- #Choose a linear solver if and only if NEWTON is choosen -+ #Choose a linear solver if and only if NEWTON is choosen - if self.options["linear_solver"] == 'DENSE': -- #Specify the use of the internal dense linear algebra functions. -- flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create a dense Sundials matrix -+ self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim) -+ #Create a dense Sundials linear solver -+ self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix) -+ #Attach it to IDA -+ flag = SUNDIALS.IDADlsSetLinearSolver(self.ida_mem, self.sun_linearsolver, self.sun_matrix); -+ ELSE: -+ #Specify the use of the internal dense linear algebra functions. -+ flag = SUNDIALS.IDADense(self.ida_mem, self.pData.dim) - if flag < 0: - raise IDAError(flag, self.t) - - elif self.options["linear_solver"] == 'SPGMR': -- #Specify the use of SPGMR linear solver. -- flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create the linear solver -+ self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], 0) -+ #Attach it to IDAS -+ flag = SUNDIALS.IDASpilsSetLinearSolver(self.ida_mem, self.sun_linearsolver) -+ ELSE: -+ #Specify the use of SPGMR linear solver. -+ flag = SUNDIALS.IDASpgmr(self.ida_mem, 0) #0 == Default krylov iterations - if flag < 0: - raise IDAError(flag, self.t) - -@@ -310,12 +330,17 @@ - if self.options["linear_solver"] == 'DENSE': - #Specify the jacobian to the solver - if self.pData.JAC != NULL and self.options["usejac"]: -- -- flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, ida_jac) -+ ELSE: -+ flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, ida_jac) - if flag < 0: - raise IDAError(flag,self.t) - else: -- flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDADlsSetJacFn(self.ida_mem, NULL) -+ ELSE: -+ flag = SUNDIALS.IDADlsSetDenseJacFn(self.ida_mem, NULL) - if flag < 0: - raise IDAError(flag,self.t) - -@@ -322,11 +347,17 @@ - elif self.options["linear_solver"] == 'SPGMR': - #Specify the jacobian times vector function - if self.pData.JACV != NULL and self.options["usejac"]: -- flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv); -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, SUNDIALS.ida_spils_jtsetup_dummy, ida_jacv); -+ ELSE: -+ flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, ida_jacv); - if flag < 0: - raise IDAError(flag, self.t) - else: -- flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL); -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.IDASpilsSetJacTimes(self.ida_mem, NULL, NULL); -+ ELSE: -+ flag = SUNDIALS.IDASpilsSetJacTimesVecFn(self.ida_mem, NULL); - if flag < 0: - raise IDAError(flag, self.t) - else: -@@ -1415,6 +1446,8 @@ - cdef public N.ndarray yS0 - #cdef N.ndarray _event_info - cdef public N.ndarray g_old -+ cdef SUNDIALS.SUNMatrix sun_matrix -+ cdef SUNDIALS.SUNLinearSolver sun_linearsolver - - def __init__(self, problem): - Explicit_ODE.__init__(self, problem) #Calls the base class -@@ -1479,6 +1512,13 @@ - if self.cvode_mem != NULL: - #Free Memory - SUNDIALS.CVodeFree(&self.cvode_mem) -+ -+ IF SUNDIALS_VERSION >= (3,0,0): -+ if self.sun_matrix != NULL: -+ SUNDIALS.SUNMatDestroy(self.sun_matrix) -+ -+ if self.sun_linearsolver != NULL: -+ SUNDIALS.SUNLinSolFree(self.sun_linearsolver) - - cpdef get_local_errors(self): - """ -@@ -2019,29 +2059,51 @@ - Updates the simulation options. - """ - cdef flag -- -+ - #Choose a linear solver if and only if NEWTON is choosen - if self.options["linear_solver"] == 'DENSE' and self.options["iter"] == "Newton": -- #Specify the use of the internal dense linear algebra functions. -- flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim) -- if flag < 0: -- raise CVodeError(flag) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create a dense Sundials matrix -+ self.sun_matrix = SUNDIALS.SUNDenseMatrix(self.pData.dim, self.pData.dim) -+ #Create a dense Sundials linear solver -+ self.sun_linearsolver = SUNDIALS.SUNDenseLinearSolver(self.yTemp, self.sun_matrix) -+ #Attach it to CVode -+ flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix); -+ if flag < 0: -+ raise CVodeError(flag) -+ ELSE: -+ #Specify the use of the internal dense linear algebra functions. -+ flag = SUNDIALS.CVDense(self.cvode_mem, self.pData.dim) -+ if flag < 0: -+ raise CVodeError(flag) - - #Specify the jacobian to the solver - if self.pData.JAC != NULL and self.options["usejac"]: -- flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac); -+ ELSE: -+ flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, cv_jac) - if flag < 0: - raise CVodeError(flag) - else: -- flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, NULL); -+ ELSE: -+ flag = SUNDIALS.CVDlsSetDenseJacFn(self.cvode_mem, NULL) - if flag < 0: - raise CVodeError(flag) - - elif self.options["linear_solver"] == 'SPGMR' and self.options["iter"] == "Newton": -- #Specify the use of CVSPGMR linear solver. -- flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"]) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ #Create the linear solver -+ self.sun_linearsolver = SUNDIALS.SUNSPGMR(self.yTemp, self.options["precond"], self.options["maxkrylov"]) -+ #Attach it to CVode -+ flag = SUNDIALS.CVSpilsSetLinearSolver(self.cvode_mem, self.sun_linearsolver) -+ ELSE: -+ #Specify the use of CVSPGMR linear solver. -+ flag = SUNDIALS.CVSpgmr(self.cvode_mem, self.options["precond"], self.options["maxkrylov"]) - if flag < 0: -- raise CVodeError(flag) -+ raise CVodeError(flag) - - if self.pData.PREC_SOLVE != NULL: - if self.pData.PREC_SETUP != NULL: -@@ -2055,11 +2117,17 @@ - - #Specify the jacobian times vector function - if self.pData.JACV != NULL and self.options["usejac"]: -- flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv) -- if flag < 0: -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, SUNDIALS.cv_spils_jtsetup_dummy, cv_jacv); -+ ELSE: -+ flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, cv_jacv) -+ if flag < 0: - raise CVodeError(flag) - else: -- flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVSpilsSetJacTimes(self.cvode_mem, NULL, NULL); -+ ELSE: -+ flag = SUNDIALS.CVSpilsSetJacTimesVecFn(self.cvode_mem, NULL) - if flag < 0: - raise CVodeError(flag) - elif self.options["linear_solver"] == 'SPARSE' and self.options["iter"] == "Newton": -@@ -2066,17 +2134,28 @@ - - if SUNDIALS.version() < (2,6,0): - raise AssimuloException("Not supported with this SUNDIALS version.") -- -+ if SUNDIALS.with_superlu() == 0: -+ raise AssimuloException("No support for SuperLU was detected, please verify that SuperLU and SUNDIALS has been installed correctly.") -+ - #Specify the use of CVSPGMR linear solver. - if self.problem_info["jac_fcn_nnz"] == -1: - raise AssimuloException("Need to specify the number of non zero elements in the Jacobian via the option 'jac_nnz'") -- flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"]) -+ -+ IF SUNDIALS_VERSION >= (3,0,0): -+ self.sun_matrix = SUNDIALS.SUNSparseMatrix(self.pData.dim, self.pData.dim, self.problem_info["jac_fcn_nnz"], CSC_MAT) -+ self.sun_linearsolver = SUNDIALS.SUNSuperLUMT(self.yTemp, self.sun_matrix, self.options["num_threads"]) -+ flag = SUNDIALS.CVDlsSetLinearSolver(self.cvode_mem, self.sun_linearsolver, self.sun_matrix) -+ ELSE: -+ flag = SUNDIALS.CVSuperLUMT(self.cvode_mem, self.options["num_threads"], self.pData.dim, self.problem_info["jac_fcn_nnz"]) - if flag < 0: - raise CVodeError(flag) - - #Specify the jacobian to the solver - if self.pData.JAC != NULL and self.options["usejac"]: -- flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsSetJacFn(self.cvode_mem, cv_jac_sparse) -+ ELSE: -+ flag = SUNDIALS.CVSlsSetSparseJacFn(self.cvode_mem, cv_jac_sparse) - if flag < 0: - raise CVodeError(flag) - else: -@@ -2899,7 +2978,10 @@ - flag = SUNDIALS.CVSpilsGetNumRhsEvals(self.cvode_mem, &nfevalsLS) #Number of rhs due to jac*vector - self.statistics["njacvecs"] += njvevals - elif self.options["linear_solver"] == "SPARSE": -- flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals) -+ IF SUNDIALS_VERSION >= (3,0,0): -+ flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals) -+ ELSE: -+ flag = SUNDIALS.CVSlsGetNumJacEvals(self.cvode_mem, &njevals) - self.statistics["njacs"] += njevals - else: - flag = SUNDIALS.CVDlsGetNumJacEvals(self.cvode_mem, &njevals) #Number of jac evals -Index: assimulo/lib/sundials_includes.pxd -=================================================================== ---- assimulo/lib/sundials_includes.pxd (revision 844) -+++ assimulo/lib/sundials_includes.pxd (revision 845) -@@ -37,19 +37,6 @@ - ctypedef double realtype
- ctypedef bint booleantype # should be bool instead of bint, but there is a bug in Cython
-
--#==============================================
--# C headers
--#==============================================
--cdef extern from "string.h":
-- void *memcpy(void *s1, void *s2, int n)
--cdef extern from "stdlib.h":
-- void *malloc(int size)
-- void free(void *ptr)
--
--#==============================================
--#External definitions from Sundials headers
--#==============================================
--
- cdef extern from "sundials/sundials_nvector.h":
- ctypedef _generic_N_Vector* N_Vector
-
-@@ -77,6 +64,107 @@ - void N_VDestroy_Serial(N_Vector v)
- void N_VPrint_Serial(N_Vector v)
-
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "sundials/sundials_types.h":
-+ IF SUNDIALS_VECTOR_SIZE == "64":
-+ ctypedef long int sunindextype
-+ ELSE:
-+ ctypedef int sunindextype
-+ cdef extern from "sundials/sundials_matrix.h":
-+ ctypedef _generic_SUNMatrix *SUNMatrix;
-+ void SUNMatDestroy(SUNMatrix A);
-+
-+ cdef struct _generic_SUNMatrix_Ops:
-+ SUNMatrix_ID (*getid)(SUNMatrix);
-+ SUNMatrix (*clone)(SUNMatrix);
-+ void (*destroy)(SUNMatrix);
-+ int (*zero)(SUNMatrix);
-+ int (*copy)(SUNMatrix, SUNMatrix);
-+ int (*scaleadd)(realtype, SUNMatrix, SUNMatrix);
-+ int (*scaleaddi)(realtype, SUNMatrix);
-+ int (*matvec)(SUNMatrix, N_Vector, N_Vector);
-+ int (*space)(SUNMatrix, long int*, long int*);
-+
-+ cdef struct _generic_SUNMatrix:
-+ void *content;
-+ _generic_SUNMatrix_Ops *ops;
-+
-+ cdef enum SUNMatrix_ID:
-+ SUNMATRIX_DENSE,
-+ SUNMATRIX_BAND,
-+ SUNMATRIX_SPARSE,
-+ SUNMATRIX_CUSTOM
-+
-+ cdef extern from "sundials/sundials_linearsolver.h":
-+ ctypedef _generic_SUNLinearSolver *SUNLinearSolver;
-+ int SUNLinSolFree(SUNLinearSolver S);
-+
-+ cdef struct _generic_SUNLinearSolver_Ops:
-+ SUNLinearSolver_Type (*gettype)(SUNLinearSolver);
-+ int (*setatimes)(SUNLinearSolver, void*, ATimesFn);
-+ int (*setpreconditioner)(SUNLinearSolver, void*,
-+ PSetupFn, PSolveFn);
-+ int (*setscalingvectors)(SUNLinearSolver,
-+ N_Vector, N_Vector);
-+ int (*initialize)(SUNLinearSolver);
-+ int (*setup)(SUNLinearSolver, SUNMatrix);
-+ int (*solve)(SUNLinearSolver, SUNMatrix, N_Vector,
-+ N_Vector, realtype);
-+ int (*numiters)(SUNLinearSolver);
-+ realtype (*resnorm)(SUNLinearSolver);
-+ long int (*lastflag)(SUNLinearSolver);
-+ int (*space)(SUNLinearSolver, long int*, long int*);
-+ N_Vector (*resid)(SUNLinearSolver);
-+ int (*free)(SUNLinearSolver);
-+
-+ cdef struct _generic_SUNLinearSolver:
-+ void *content;
-+ _generic_SUNLinearSolver_Ops *ops;
-+
-+ cdef enum SUNLinearSolver_Type:
-+ SUNLINEARSOLVER_DIRECT,
-+ SUNLINEARSOLVER_ITERATIVE,
-+ SUNLINEARSOLVER_CUSTOM
-+
-+ cdef extern from "sunmatrix/sunmatrix_dense.h":
-+ ctypedef _SUNMatrixContent_Dense *SUNMatrixContent_Dense;
-+ cdef struct _SUNMatrixContent_Dense:
-+ sunindextype M;
-+ sunindextype N;
-+ realtype *data;
-+ sunindextype ldata;
-+ realtype **cols;
-+ SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N);
-+ cdef extern from "sunmatrix/sunmatrix_sparse.h":
-+ ctypedef _SUNMatrixContent_Sparse *SUNMatrixContent_Sparse;
-+ cdef struct _SUNMatrixContent_Sparse:
-+ sunindextype M;
-+ sunindextype N;
-+ sunindextype NNZ;
-+ sunindextype NP;
-+ realtype *data;
-+ int sparsetype;
-+ sunindextype *indexvals;
-+ sunindextype *indexptrs;
-+ sunindextype **rowvals;
-+ sunindextype **colptrs;
-+ sunindextype **colvals;
-+ sunindextype **rowptrs;
-+ SUNMatrix SUNSparseMatrix(sunindextype M, sunindextype N, sunindextype NNZ, int sparsetype);
-+ cdef extern from "sunlinsol/sunlinsol_dense.h":
-+ SUNLinearSolver SUNDenseLinearSolver(N_Vector y, SUNMatrix A);
-+ cdef extern from "sunlinsol/sunlinsol_spgmr.h":
-+ SUNLinearSolver SUNSPGMR(N_Vector y, int pretype, int maxl);
-+
-+ELSE:
-+ #Dummy defines
-+ ctypedef void *SUNLinearSolver
-+ ctypedef void *SUNMatrix
-+ ctypedef void *SUNMatrixContent_Dense
-+ ctypedef void *SUNMatrixContent_Sparse
-+ ctypedef int sunindextype
-+
-+
- #Struct for handling the Jacobian data
- cdef extern from "sundials/sundials_direct.h":
- cdef struct _DlsMat:
-@@ -92,7 +180,7 @@ - realtype **cols
- ctypedef _DlsMat* DlsMat
- cdef realtype* DENSE_COL(DlsMat A, int j)
--
-+
- IF SUNDIALS_VERSION >= (2,6,3):
- cdef extern from "sundials/sundials_sparse.h":
- cdef struct _SlsMat:
-@@ -128,7 +216,26 @@ - int *rowvals
- int *colptrs
- ctypedef _SlsMat* SlsMat
-+
-+#==============================================
-+# C headers
-+#==============================================
-+cdef extern from "string.h":
-+ void *memcpy(void *s1, void *s2, int n)
-+cdef extern from "stdlib.h":
-+ void *malloc(int size)
-+ void free(void *ptr)
-+
-+#==============================================
-+#External definitions from Sundials headers
-+#==============================================
-
-+IF SUNDIALS_WITH_SUPERLU:
-+ cdef inline int with_superlu(): return 1
-+ELSE:
-+ cdef inline int with_superlu(): return 0
-+
-+
- cdef extern from "cvodes/cvodes.h":
- void* CVodeCreate(int lmm, int iter)
- ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void *f_data)
-@@ -225,45 +332,67 @@ - int CVodeGetSensNonlinSolvStats(void *cvode_mem, long int *nSniters, long int *nSncfails)
- int CVodeGetStgrSensNumNonlinSolvIters(void *cvode_mem, long int *nSTGR1niters)
- int CVodeGetStgrSensNumNonlinSolvConvFails(void *cvode_mem, long int *nSTGR1ncfails)
-+
-+cdef extern from "cvodes/cvodes_spils.h":
-+ ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
-+ N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "cvodes/cvodes_direct.h":
-+ ctypedef int (*CVDlsDenseJacFn)(realtype t, N_Vector y, N_Vector fy,
-+ SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int CVDlsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS, SUNMatrix A);
-+ int CVDlsSetJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
-+ cdef extern from "cvodes/cvodes_spils.h":
-+ int CVSpilsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS);
-+ ctypedef int (*CVSpilsJacTimesSetupFn)(realtype t, N_Vector y, N_Vector fy, void *user_data);
-+ int CVSpilsSetJacTimes(void *cvode_mem, CVSpilsJacTimesSetupFn jtsetup, CVSpilsJacTimesVecFn jtimes);
-
-
--cdef extern from "cvodes/cvodes_dense.h":
-- int CVDense(void *cvode_mem, long int n)
-- ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy,
-- DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-- int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
-+ IF SUNDIALS_WITH_SUPERLU:
-+ cdef extern from "sunlinsol/sunlinsol_superlumt.h":
-+ SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads)
-+ ELSE:
-+ cdef inline SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads): return NULL
-+
-+ cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0
-+ cdef inline tuple version(): return (3,0,0)
-+ELSE:
-+ cdef extern from "cvodes/cvodes_dense.h":
-+ int CVDense(void *cvode_mem, long int n)
-+ ctypedef int (*CVDlsDenseJacFn)(int n, realtype t, N_Vector y, N_Vector fy,
-+ DlsMat Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int CVDlsSetDenseJacFn(void *cvode_mem, CVDlsDenseJacFn djac)
-
--cdef extern from "cvodes/cvodes_spgmr.h":
-- int CVSpgmr(void *cvode_mem, int pretype, int max1)
-+ cdef extern from "cvodes/cvodes_spgmr.h":
-+ int CVSpgmr(void *cvode_mem, int pretype, int max1)
-
--IF SUNDIALS_VERSION >= (2,6,0):
-- cdef extern from "cvodes/cvodes_sparse.h":
-+ cdef extern from "cvodes/cvodes_spils.h":
-+ int CVSpilsSetJacTimesVecFn(void *cvode_mem, CVSpilsJacTimesVecFn jtv)
-+
-+ IF SUNDIALS_VERSION >= (2,6,0):
-+ cdef extern from "cvodes/cvodes_sparse.h":
-+ ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
-+ SlsMat Jac, void *user_data, N_Vector tmp1,
-+ N_Vector tmp2, N_Vector tmp3)
-+ int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
-+ int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
-+ cdef inline tuple version(): return (2,6,0)
-+ IF SUNDIALS_WITH_SUPERLU:
-+ cdef extern from "cvodes/cvodes_sparse.h":
-+ int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
-+ ELSE:
-+ cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
-+ ELSE:
-+ cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
- ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
- SlsMat Jac, void *user_data, N_Vector tmp1,
- N_Vector tmp2, N_Vector tmp3)
-- int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac)
-- int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals)
-- #cdef inline char* version(): return "2.6.0"
-- cdef inline tuple version(): return (2,6,0)
-- IF SUNDIALS_WITH_SUPERLU:
-- cdef extern from "cvodes/cvodes_sparse.h":
-- int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz)
-- ELSE:
-- cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
--ELSE:
-- cdef inline int CVSuperLUMT(void *cvode_mem, int numthreads, int n, int nnz): return -1
-- ctypedef int (*CVSlsSparseJacFn)(realtype t, N_Vector y, N_Vector fy,
-- SlsMat Jac, void *user_data, N_Vector tmp1,
-- N_Vector tmp2, N_Vector tmp3)
-- cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
-- cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
-- #cdef inline char* version(): return "2.5.0"
-- cdef inline tuple version(): return (2,5,0)
-+ cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
-+ cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
-+ cdef inline tuple version(): return (2,5,0)
-
- cdef extern from "cvodes/cvodes_spils.h":
-- ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
-- N_Vector y, N_Vector fy,
-- void *user_data, N_Vector tmp)
- ctypedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy,
- booleantype jok, booleantype *jcurPtr,
- realtype gamma, void *user_data,
-@@ -273,13 +402,12 @@ - N_Vector r, N_Vector z,
- realtype gamma, realtype delta,
- int lr, void *user_data, N_Vector tmp)
-- int CVSpilsSetJacTimesVecFn(void *cvode_mem, CVSpilsJacTimesVecFn jtv)
-+
- int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
- int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
- int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
- int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals)
- int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves)
--
-
- cdef extern from "idas/idas.h":
- ctypedef int (*IDAResFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
-@@ -380,19 +508,40 @@ -
- #End Sensitivities
- #=================
--
--cdef extern from "idas/idas_dense.h":
-- int IDADense(void *ida_mem, long int n)
-- ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy,
-- N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data,
-- N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-- int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
-
--cdef extern from "idas/idas_dense.h":
-- int IDASpgmr(void *ida_mem, int max1)
-- ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
-- int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
-+cdef extern from "idas/idas_spils.h":
-+ ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr,
-+ N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
-
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "idas/idas_direct.h":
-+ ctypedef int (*IDADlsDenseJacFn)(realtype tt, realtype cj, N_Vector yy,
-+ N_Vector yp, N_Vector rr, SUNMatrix Jac, void *user_data,
-+ N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int IDADlsSetJacFn(void *ida_mem, IDADlsDenseJacFn djac)
-+ int IDADlsSetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix A);
-+
-+ cdef extern from "idas/idas_spils.h":
-+ int IDASpilsSetLinearSolver(void *ida_mem, SUNLinearSolver LS);
-+ ctypedef int (*IDASpilsJacTimesSetupFn)(realtype tt, N_Vector yy,
-+ N_Vector yp, N_Vector rr, realtype c_j, void *user_data);
-+ int IDASpilsSetJacTimes(void *ida_mem,
-+ IDASpilsJacTimesSetupFn jtsetup, IDASpilsJacTimesVecFn jtimes);
-+
-+ cdef inline int ida_spils_jtsetup_dummy(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, realtype c_j, void *user_data): return 0
-+ELSE:
-+ cdef extern from "idas/idas_dense.h":
-+ int IDADense(void *ida_mem, long int n)
-+ ctypedef int (*IDADlsDenseJacFn)(int Neq, realtype tt, realtype cj, N_Vector yy,
-+ N_Vector yp, N_Vector rr, DlsMat Jac, void *user_data,
-+ N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-+ int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn djac)
-+
-+ int IDASpgmr(void *ida_mem, int max1)
-+
-+ cdef extern from "idas/idas_spils.h":
-+ int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
-+
- cdef extern from "idas/idas_spils.h":
- int IDASpilsGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
- int IDASpilsGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
-@@ -455,14 +604,31 @@ - # fuction used to deallocate memory used by KINSOL
- void KINFree(void **kinmem)
-
--# functions used for supplying jacobian, and receiving info from linear solver
--cdef extern from "kinsol/kinsol_direct.h":
-- # user functions
-- ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef extern from "kinsol/kinsol_direct.h":
-+ ctypedef int (*KINDlsDenseJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
-+ int KINDlsSetLinearSolver(void *kinmem, SUNLinearSolver LS, SUNMatrix A);
-+ int KINDlsSetJacFn(void *kinmem, KINDlsDenseJacFn djac)
-
-- # function used to link user functions to KINSOL
-- int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
-+ cdef extern from "kinsol/kinsol_spils.h":
-+ int KINSpilsSetLinearSolver(void *kinsol_mem, SUNLinearSolver LS);
-+ELSE:
-+ # functions used for supplying jacobian, and receiving info from linear solver
-+ cdef extern from "kinsol/kinsol_direct.h":
-+ # user functions
-+ ctypedef int (*KINDlsDenseJacFn)(int dim, N_Vector u, N_Vector fu, DlsMat J, void *user_data, N_Vector tmp1, N_Vector tmp2)
-+
-+ # function used to link user functions to KINSOL
-+ int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
-+
-+ cdef extern from "kinsol/kinsol_dense.h":
-+ int KINDense(void *kinmem, int dim)
-+
-+ cdef extern from "kinsol/kinsol_spgmr.h":
-+ int KINSpgmr(void *kinmem, int maxl)
-
-+cdef extern from "kinsol/kinsol_direct.h":
- # optional output fcts for linear direct solver
- int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB)
- int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB)
-@@ -470,12 +636,6 @@ - int KINDlsGetLastFlag(void *kinmem, int *flag)
- char *KINDlsGetReturnFlagName(int flag)
-
--cdef extern from "kinsol/kinsol_dense.h":
-- int KINDense(void *kinmem, int dim)
--
--cdef extern from "kinsol/kinsol_spgmr.h":
-- int KINSpgmr(void *kinmem, int maxl)
--
- cdef extern from "kinsol/kinsol_spils.h":
- ctypedef int (*KINSpilsJacTimesVecFn)(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
- void *problem_data)
-Index: src/lib/sundials_constants.pxi -=================================================================== ---- assimulo/lib/sundials_constants.pxi (revision 844) -+++ assimulo/lib/sundials_constants.pxi (revision 845) -@@ -95,6 +95,9 @@ - DEF CV_BAD_IS = -45
-
-
-+DEF CSC_MAT = 0
-+DEF CSR_MAT = 1
-+
- #==========
- # IDA
- #==========
-Index: src/lib/sundials_callbacks.pxi -=================================================================== ---- assimulo/lib/sundials_callbacks.pxi (revision 844) -+++ assimulo/lib/sundials_callbacks.pxi (revision 845) -@@ -94,98 +94,55 @@ - traceback.print_exc()
- return CV_UNREC_RHSFUNC_ERR
-
--@cython.boundscheck(False)
--@cython.wraparound(False)
--cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
-- void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-- """
-- This method is used to connect the Assimulo.Problem.jac to the Sundials
-- Sparse Jacobian function.
-- """
-- cdef ProblemData pData = <ProblemData>problem_data
-- #cdef N.ndarray y = nv2arr(yv)
-- cdef N.ndarray y = pData.work_y
-- cdef int i
-- cdef int nnz = Jacobian.NNZ
-- cdef int ret_nnz
-- cdef int dim = Jacobian.N
-- cdef realtype* data = Jacobian.data
--
-- IF SUNDIALS_VERSION >= (2,6,3):
-- cdef int* rowvals = Jacobian.rowvals[0]
-- cdef int* colptrs = Jacobian.colptrs[0]
-- ELSE:
-- cdef int* rowvals = Jacobian.rowvals
-- cdef int* colptrs = Jacobian.colptrs
--
-- nv2arr_inplace(yv, y)
-- """
-- realtype *data;
-- int *rowvals;
-- int *colptrs;
-- """
-- try:
-- if pData.dimSens > 0: #Sensitivity activated
-- p = realtype2arr(pData.p,pData.dimSens)
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-- else:
-- jac=(<object>pData.JAC)(t,y,p=p)
-- else:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-- else:
-- jac=(<object>pData.JAC)(t,y)
--
-- if not isinstance(jac, sparse.csc.csc_matrix):
-- jac = sparse.csc.csc_matrix(jac)
-- raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-- ret_nnz = jac.nnz
-- if ret_nnz > nnz:
-- raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")
-
-- for i in range(min(ret_nnz,nnz)):
-- data[i] = jac.data[i]
-- rowvals[i] = jac.indices[i]
-- for i in range(dim+1):
-- colptrs[i] = jac.indptr[i]
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ @cython.boundscheck(False)
-+ @cython.wraparound(False)
-+ cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Sparse Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef SUNMatrixContent_Sparse Jacobian = <SUNMatrixContent_Sparse>Jac.content
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i
-+ cdef sunindextype nnz = Jacobian.NNZ
-+ cdef int ret_nnz
-+ cdef sunindextype dim = Jacobian.N
-+ cdef realtype* data = Jacobian.data
-+ cdef sunindextype* rowvals = Jacobian.rowvals[0]
-+ cdef sunindextype* colptrs = Jacobian.colptrs[0]
-
-- return CVDLS_SUCCESS
-- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-- return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-- except:
-- traceback.print_exc()
-- return CVDLS_JACFUNC_UNRECVR
-+ nv2arr_inplace(yv, y)
-
--cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian,
-- void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-- """
-- This method is used to connect the Assimulo.Problem.jac to the Sundials
-- Jacobian function.
-- """
-- cdef ProblemData pData = <ProblemData>problem_data
-- #cdef ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-- cdef realtype* col_i=DENSE_COL(Jacobian,0)
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #cdef N.ndarray y = nv2arr(yv)
-- cdef N.ndarray y = pData.work_y
-- cdef int i,j
--
-- nv2arr_inplace(yv, y)
--
-- if pData.dimSens>0: #Sensitivity activated
-- p = realtype2arr(pData.p,pData.dimSens)
- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
-+ if pData.dimSens > 0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p=p)
- else:
-- jac=(<object>pData.JAC)(t,y,p)
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-+ if not isinstance(jac, sparse.csc.csc_matrix):
-+ jac = sparse.csc.csc_matrix(jac)
-+ raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-+ ret_nnz = jac.nnz
-+ if ret_nnz > nnz:
-+ raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")
-
-+ for i in range(min(ret_nnz,nnz)):
-+ data[i] = jac.data[i]
-+ rowvals[i] = jac.indices[i]
-+ for i in range(dim+1):
-+ colptrs[i] = jac.indptr[i]
-+
- return CVDLS_SUCCESS
- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
- return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-@@ -192,18 +149,62 @@ - except:
- traceback.print_exc()
- return CVDLS_JACFUNC_UNRECVR
-- else:
-+ELSE:
-+ @cython.boundscheck(False)
-+ @cython.wraparound(False)
-+ cdef int cv_jac_sparse(realtype t, N_Vector yv, N_Vector fy, SlsMat Jacobian,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Sparse Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i
-+ cdef int nnz = Jacobian.NNZ
-+ cdef int ret_nnz
-+ cdef int dim = Jacobian.N
-+ cdef realtype* data = Jacobian.data
-+
-+ IF SUNDIALS_VERSION >= (2,6,3):
-+ cdef int* rowvals = Jacobian.rowvals[0]
-+ cdef int* colptrs = Jacobian.colptrs[0]
-+ ELSE:
-+ cdef int* rowvals = Jacobian.rowvals
-+ cdef int* colptrs = Jacobian.colptrs
-+
-+ nv2arr_inplace(yv, y)
-+ """
-+ realtype *data;
-+ int *rowvals;
-+ int *colptrs;
-+ """
- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ if pData.dimSens > 0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,p=p,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p=p)
- else:
-- jac=(<object>pData.JAC)(t,y)
--
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-+
-+ if not isinstance(jac, sparse.csc.csc_matrix):
-+ jac = sparse.csc.csc_matrix(jac)
-+ raise AssimuloException("The Jacobian must be stored on Scipy's CSC format.")
-+ ret_nnz = jac.nnz
-+ if ret_nnz > nnz:
-+ raise AssimuloException("The Jacobian has more entries than supplied to the problem class via 'jac_nnz'")
-
-+ for i in range(min(ret_nnz,nnz)):
-+ data[i] = jac.data[i]
-+ rowvals[i] = jac.indices[i]
-+ for i in range(dim+1):
-+ colptrs[i] = jac.indptr[i]
-+
- return CVDLS_SUCCESS
- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
- return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-@@ -210,7 +211,112 @@ - except:
- traceback.print_exc()
- return CVDLS_JACFUNC_UNRECVR
-+
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef int cv_jac(realtype t, N_Vector yv, N_Vector fy, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef realtype* col_i=Jacobian.cols[0]
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i,j, Neq = pData.dim
-
-+ nv2arr_inplace(yv, y)
-+
-+ if pData.dimSens>0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+ELSE:
-+ cdef int cv_jac(int Neq, realtype t, N_Vector yv, N_Vector fy, DlsMat Jacobian,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef realtype* col_i=DENSE_COL(Jacobian,0)
-+ cdef N.ndarray y = pData.work_y
-+ cdef int i,j
-+
-+ nv2arr_inplace(yv, y)
-+
-+ if pData.dimSens>0: #Sensitivity activated
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw,p=p)
-+ else:
-+ jac=(<object>pData.JAC)(t,y,p)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(t,y,sw=<list>pData.sw)
-+ else:
-+ jac=(<object>pData.JAC)(t,y)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return CVDLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return CVDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ except:
-+ traceback.print_exc()
-+ return CVDLS_JACFUNC_UNRECVR
-+
- cdef int cv_jacv(N_Vector vv, N_Vector Jv, realtype t, N_Vector yv, N_Vector fyv,
- void *problem_data, N_Vector tmp):
- """
-@@ -337,9 +443,6 @@ - Root-finding function.
- """
- cdef ProblemData pData = <ProblemData>problem_data
-- #cdef ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #cdef N.ndarray y = nv2arr(yv)
- cdef N.ndarray y = pData.work_y
- cdef int i
-
-@@ -366,12 +469,8 @@ - """
- cdef ProblemData pData = <ProblemData>problem_data
- cdef N.ndarray[realtype, ndim=1, mode='c'] res #Used for return from the user function
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #(<ndarray>pData.yd).data = <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
- cdef N.ndarray y = pData.work_y
- cdef N.ndarray yd = pData.work_yd
-- # cdef N.ndarray y = nv2arr(yv)
-- # cdef N.ndarray yd = nv2arr(yvdot)
- cdef realtype* resptr=(<N_VectorContent_Serial>residual.content).data
- cdef int i
-
-@@ -414,63 +513,113 @@ - except:
- traceback.print_exc()
- return IDA_RES_FAIL
--
--cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
-+
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef int ida_jac(realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-+ cdef realtype* col_i=Jacobian.cols[0]
-+ cdef N.ndarray y = pData.work_y
-+ cdef N.ndarray yd = pData.work_yd
-+ cdef int i,j, Neq = pData.dim
-+
-+ nv2arr_inplace(yv, y)
-+ nv2arr_inplace(yvdot, yd)
-+
-+ if pData.dimSens!=0: #SENSITIVITY
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd,p=p)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd)
-+
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-+ELSE:
-+ cdef int ida_jac(int Neq, realtype t, realtype c, N_Vector yv, N_Vector yvdot, N_Vector residual, DlsMat Jacobian,
- void* problem_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3):
-- """
-- This method is used to connect the Assimulo.Problem.jac to the Sundials
-- Jacobian function.
-- """
-- cdef ProblemData pData = <ProblemData>problem_data
-- cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-- cdef realtype* col_i=DENSE_COL(Jacobian,0)
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #(<ndarray>pData.yd).data = <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
-- cdef N.ndarray y = pData.work_y
-- cdef N.ndarray yd = pData.work_yd
-- #cdef N.ndarray y = nv2arr(yv)
-- #cdef N.ndarray yd = nv2arr(yvdot)
-- cdef int i,j
--
-- nv2arr_inplace(yv, y)
-- nv2arr_inplace(yvdot, yd)
--
-- if pData.dimSens!=0: #SENSITIVITY
-- p = realtype2arr(pData.p,pData.dimSens)
-- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p) # call to the python residual function
-- else:
-- jac=(<object>pData.JAC)(c,t,y,yd,p=p)
-+ """
-+ This method is used to connect the Assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef ProblemData pData = <ProblemData>problem_data
-+ cdef N.ndarray[realtype, ndim=2, mode='c'] jac #Used for return from the user function
-+ cdef realtype* col_i=DENSE_COL(Jacobian,0)
-+ cdef N.ndarray y = pData.work_y
-+ cdef N.ndarray yd = pData.work_yd
-+ cdef int i,j
-+
-+ nv2arr_inplace(yv, y)
-+ nv2arr_inplace(yvdot, yd)
-+
-+ if pData.dimSens!=0: #SENSITIVITY
-+ p = realtype2arr(pData.p,pData.dimSens)
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,sw=<list>pData.sw,p=p) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd,p=p)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-+ else:
-+ try:
-+ if pData.sw != NULL:
-+ jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw) # call to the python residual function
-+ else:
-+ jac=(<object>pData.JAC)(c,t,y,yd)
-+
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+ return IDADLS_SUCCESS
-+ except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-+ return IDADLS_JACFUNC_RECVR #Recoverable Error
-+ except:
-+ traceback.print_exc()
-+ return IDADLS_JACFUNC_UNRECVR
-
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-- return IDADLS_SUCCESS
-- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-- return IDADLS_JACFUNC_RECVR #Recoverable Error
-- except:
-- traceback.print_exc()
-- return IDADLS_JACFUNC_UNRECVR
-- else:
-- try:
-- if pData.sw != NULL:
-- jac=(<object>pData.JAC)(c,t,y,yd,<list>pData.sw) # call to the python residual function
-- else:
-- jac=(<object>pData.JAC)(c,t,y,yd)
--
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-- return IDADLS_SUCCESS
-- except(N.linalg.LinAlgError,ZeroDivisionError,AssimuloRecoverableError):
-- return IDADLS_JACFUNC_RECVR #Recoverable Error
-- except:
-- traceback.print_exc()
-- return IDADLS_JACFUNC_UNRECVR
--
-
- cdef int ida_root(realtype t, N_Vector yv, N_Vector yvdot, realtype *gout, void* problem_data):
- """
-@@ -479,12 +628,8 @@ - """
- cdef ProblemData pData = <ProblemData>problem_data
- cdef N.ndarray[realtype, ndim=1, mode='c'] root #Used for return from the user function
-- #(<ndarray>pData.y).data = <realtype*>((<N_VectorContent_Serial>yv.content).data)
-- #(<ndarray>pData.yd).data = <realtype*>((<N_VectorContent_Serial>yvdot.content).data)
- cdef N.ndarray y = pData.work_y
- cdef N.ndarray yd = pData.work_yd
-- #cdef N.ndarray y = nv2arr(yv)
-- #cdef N.ndarray yd = nv2arr(yvdot)
- cdef int i
-
- nv2arr_inplace(yv, y)
-@@ -553,30 +698,54 @@ - traceback.print_exc()
- return SPGMR_PSOLVE_FAIL_UNREC
-
-+IF SUNDIALS_VERSION >= (3,0,0):
-+ cdef int kin_jac(N_Vector xv, N_Vector fval, SUNMatrix Jac,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2):
-+ """
-+ This method is used to connect the assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef SUNMatrixContent_Dense Jacobian = <SUNMatrixContent_Dense>Jac.content
-+ cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-+ cdef realtype* col_i=Jacobian.cols[0]
-+ cdef N.ndarray x = nv2arr(xv)
-+ cdef int i,j, Neq = pData.dim
-+
-+ try:
-+ jac=(<object>pData.JAC)(x)
-
--cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian,
-- void *problem_data, N_Vector tmp1, N_Vector tmp2):
-- """
-- This method is used to connect the assimulo.Problem.jac to the Sundials
-- Jacobian function.
-- """
-- cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-- cdef realtype* col_i=DENSE_COL(Jacobian,0)
-- cdef N.ndarray x = nv2arr(xv)
-- cdef int i,j
--
-- try:
-- jac=(<object>pData.JAC)(x)
-+ for i in range(Neq):
-+ col_i = Jacobian.cols[i]
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-
-- for i in range(Neq):
-- col_i = DENSE_COL(Jacobian, i)
-- for j in range(Neq):
-- col_i[j] = jac[j,i]
-+ return KINDLS_SUCCESS
-+ except:
-+ return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+ELSE:
-+ cdef int kin_jac(int Neq, N_Vector xv, N_Vector fval, DlsMat Jacobian,
-+ void *problem_data, N_Vector tmp1, N_Vector tmp2):
-+ """
-+ This method is used to connect the assimulo.Problem.jac to the Sundials
-+ Jacobian function.
-+ """
-+ cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-+ cdef realtype* col_i=DENSE_COL(Jacobian,0)
-+ cdef N.ndarray x = nv2arr(xv)
-+ cdef int i,j
-+
-+ try:
-+ jac=(<object>pData.JAC)(x)
-
-- return KINDLS_SUCCESS
-- except:
-- return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
--
-+ for i in range(Neq):
-+ col_i = DENSE_COL(Jacobian, i)
-+ for j in range(Neq):
-+ col_i[j] = jac[j,i]
-+
-+ return KINDLS_SUCCESS
-+ except:
-+ return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-+
- cdef int kin_jacv(N_Vector vv, N_Vector Jv, N_Vector vx, bint new_u,
- void *problem_data):
- cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
-Index: examples/cvode_with_jac.py -=================================================================== ---- assimulo/examples/cvode_with_jac.py (revision 844) -+++ assimulo/examples/cvode_with_jac.py (revision 845) -@@ -46,7 +46,6 @@ - def f(t,y): - yd_0 = y[1] - yd_1 = -9.82 -- #print y, yd_0, yd_1 - return N.array([yd_0,yd_1]) - - #Defines the Jacobian -Index: examples/__init__.py -=================================================================== ---- assimulo/examples/__init__.py (revision 844) -+++ assimulo/examples/__init__.py (revision 845) -@@ -24,6 +24,6 @@ - "mech_system_pendulum", "euler_vanderpol", "cvode_with_parameters_modified", - "cvode_basic_backward","ida_basic_backward","dasp3_basic", "cvode_with_preconditioning", - "kinsol_basic","kinsol_with_jac", "radau5dae_time_events", "kinsol_ors", "lsodar_bouncing_ball", -- "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result"] -+ "cvode_with_parameters_fcn", "ida_with_user_defined_handle_result", "cvode_with_jac_sparse"] - - -Index: examples/cvode_with_jac_sparse.py -=================================================================== ---- assimulo/examples/cvode_with_jac_sparse.py (nonexistent) -+++ assimulo/examples/cvode_with_jac_sparse.py (revision 845) -@@ -0,0 +1,102 @@ -+#!/usr/bin/env python -+# -*- coding: utf-8 -*- -+ -+# Copyright (C) 2010 Modelon AB -+# -+# This program is free software: you can redistribute it and/or modify -+# it under the terms of the GNU Lesser General Public License as published by -+# the Free Software Foundation, version 3 of the License. -+# -+# This program is distributed in the hope that it will be useful, -+# but WITHOUT ANY WARRANTY; without even the implied warranty of -+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -+# GNU Lesser General Public License for more details. -+# -+# You should have received a copy of the GNU Lesser General Public License -+# along with this program. If not, see <http://www.gnu.org/licenses/>. -+ -+import numpy as N -+import scipy.sparse as SP -+import pylab as P -+import nose -+from assimulo.solvers import CVode -+from assimulo.problem import Explicit_Problem -+ -+ -+def run_example(with_plots=True): -+ r""" -+ Example for demonstrating the use of a user supplied Jacobian (sparse). -+ Note that this will only work if Assimulo has been configured with -+ Sundials + SuperLU. Based on the SUNDIALS example cvRoberts_sps.c -+ -+ ODE: -+ -+ .. math:: -+ -+ \dot y_1 &= -0.04y_1 + 1e4 y_2 y_3 \\ -+ \dot y_2 &= - \dot y_1 - \dot y_3 \\ -+ \dot y_3 &= 3e7 y_2^2 -+ -+ -+ on return: -+ -+ - :dfn:`exp_mod` problem instance -+ -+ - :dfn:`exp_sim` solver instance -+ -+ """ -+ -+ #Defines the rhs -+ def f(t,y): -+ yd_0 = -0.04*y[0] + 1e4*y[1]*y[2] -+ yd_2 = 3e7*y[1]*y[1] -+ yd_1 = -yd_0 - yd_2 -+ return N.array([yd_0,yd_1,yd_2]) -+ -+ #Defines the Jacobian -+ def jac(t,y): -+ -+ colptrs = [0,3,6,9] -+ rowvals = [0, 1, 2, 0, 1, 2, 0, 1, 2] -+ data = [-0.04, 0.04, 0.0, 1e4*y[2], -1e4*y[2]-6e7*y[1], 6e7*y[1], 1e4*y[1], -1e4*y[1], 0.0] -+ -+ J = SP.csc_matrix((data, rowvals, colptrs)) -+ return J -+ -+ #Defines an Assimulo explicit problem -+ y0 = [1.0,0.0,0.0] #Initial conditions -+ -+ exp_mod = Explicit_Problem(f,y0, name = 'Example using analytic (sparse) Jacobian') -+ -+ exp_mod.jac = jac #Sets the Jacobian -+ exp_mod.jac_nnz = 9 -+ -+ -+ exp_sim = CVode(exp_mod) #Create a CVode solver -+ -+ #Set the parameters -+ exp_sim.iter = 'Newton' #Default 'FixedPoint' -+ exp_sim.discr = 'BDF' #Default 'Adams' -+ exp_sim.atol = [1e-8,1e-14,1e-6] #Default 1e-6 -+ exp_sim.rtol = 1e-4 #Default 1e-6 -+ exp_sim.linear_solver = "sparse" -+ -+ #Simulate -+ t, y = exp_sim.simulate(0.4) #Simulate 0.4 seconds -+ -+ #Basic tests -+ nose.tools.assert_almost_equal(y[-1][0],0.9851,3) -+ -+ #Plot -+ if with_plots: -+ P.plot(t,y[:,1],linestyle="dashed",marker="o") #Plot the solution -+ P.xlabel('Time') -+ P.ylabel('State') -+ P.title(exp_mod.name) -+ P.show() -+ -+ return exp_mod, exp_sim -+ -+ -+if __name__=='__main__': -+ mod,sim = run_example() -Index: tests/test_examples.py -=================================================================== ---- assimulo/tests/test_examples.py (revision 844) -+++ assimulo/tests/test_examples.py (revision 845) -@@ -22,7 +22,15 @@ -
- class Test_Examples:
-
-+
- @testattr(stddist = True)
-+ def test_cvode_with_jac_sparse(self):
-+ try:
-+ cvode_with_jac_sparse.run_example(with_plots=False)
-+ except AssimuloException:
-+ pass #Handle the case when SuperLU is not installed
-+
-+ @testattr(stddist = True)
- def test_ida_with_user_defined_handle_result(self):
- ida_with_user_defined_handle_result.run_example(with_plots=False)
-
diff --git a/sundials4.patch b/sundials4.patch new file mode 100644 index 000000000000..05a302cadae8 --- /dev/null +++ b/sundials4.patch @@ -0,0 +1,26 @@ +Index: src/lib/sundials_includes.pxd +=================================================================== +--- src/lib/sundials_includes.pxd (revision 875) ++++ src/lib/sundials_includes.pxd (working copy) +@@ -237,7 +237,7 @@ +
+
+ cdef extern from "cvodes/cvodes.h":
+- void* CVodeCreate(int lmm, int iter)
++ void* CVodeCreate(int lmm)
+ ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void *f_data)
+ int CVodeInit(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0)
+ int CVodeReInit(void *cvode_mem, realtype t0, N_Vector y0)
+Index: src/solvers/sundials.pyx +=================================================================== +--- src/solvers/sundials.pyx (revision 875) ++++ src/solvers/sundials.pyx (working copy) +@@ -1699,7 +1699,7 @@ + if self.cvode_mem == NULL: #The solver is not initialized + + #Create the solver +- self.cvode_mem= SUNDIALS.CVodeCreate(CV_BDF if self.options["discr"] == "BDF" else CV_ADAMS, CV_NEWTON if self.options["iter"] == "Newton" else CV_FUNCTIONAL) ++ self.cvode_mem= SUNDIALS.CVodeCreate(CV_BDF if self.options["discr"] == "BDF" else CV_ADAMS) + if self.cvode_mem == NULL: + raise CVodeError(CV_MEM_FAIL) + |