diff --git a/.travis.yml b/.travis.yml index 750b4d3..8be320d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -31,53 +31,76 @@ matrix: # needed to work around https://github.com/travis-ci/travis-ci/issues/4794 - python: 2.7 env: - - TOXENV=py27 SUNDIALS_VERSION='2.7.0' - - python: 3.3 - env: - - TOXENV=py33 SUNDIALS_VERSION='2.7.0' + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' - python: 3.4 env: - - TOXENV=py34 SUNDIALS_VERSION='2.7.0' + - TOXENV=py34 SUNDIALS_VERSION='3.1.1' - python: 3.5 env: - - TOXENV=py35 SUNDIALS_VERSION='2.7.0' + - TOXENV=py35 SUNDIALS_VERSION='3.1.1' - python: 3.6 env: - - TOXENV=py36 SUNDIALS_VERSION='2.7.0' + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' + # commented as 3.7 not yet on travis + #- python: 3.7 + # env: + # - TOXENV=py37 SUNDIALS_VERSION='3.1.1' - python: 3.5 env: - - TOXENV=doctr SUNDIALS_VERSION='2.7.0' + - TOXENV=doctr SUNDIALS_VERSION='3.1.1' - python: 3.5 env: - - TOXENV=docs SUNDIALS_VERSION='2.7.0' + - TOXENV=docs SUNDIALS_VERSION='3.1.1' - python: 3.5 env: - - TOXENV=check-manifest SUNDIALS_VERSION='2.7.0' + - TOXENV=check-manifest SUNDIALS_VERSION='3.1.1' - python: 3.5 env: - - TOXENV=checkreadme SUNDIALS_VERSION='2.7.0' + - TOXENV=checkreadme SUNDIALS_VERSION='3.1.1' + + # Additional precisions + - python: 2.7 + env: + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' + - python: 3.6 + env: + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' + - python: 2.7 + env: + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='extended' + - python: 3.6 + env: + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='extended' + + allow_failures: + - python: 2.7 + env: + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' + - python: 3.6 + env: + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' # Additional precisions - python: 2.7 env: - - TOXENV=py27 SUNDIALS_VERSION='2.7.0' SUNDIALS_PRECISION='single' + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' - python: 3.6 env: - - TOXENV=py36 SUNDIALS_VERSION='2.7.0' SUNDIALS_PRECISION='single' + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' - python: 2.7 env: - - TOXENV=py27 SUNDIALS_VERSION='2.7.0' SUNDIALS_PRECISION='extended' + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='extended' - python: 3.6 env: - - TOXENV=py36 SUNDIALS_VERSION='2.7.0' SUNDIALS_PRECISION='extended' + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='extended' allow_failures: - python: 2.7 env: - - TOXENV=py27 SUNDIALS_VERSION='2.7.0' SUNDIALS_PRECISION='single' + - TOXENV=py27 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' - python: 3.6 env: - - TOXENV=py36 SUNDIALS_VERSION='2.7.0' SUNDIALS_PRECISION='single' + - TOXENV=py36 SUNDIALS_VERSION='3.1.1' SUNDIALS_PRECISION='single' install: - source ci_support/ensure_sundials_installed.sh diff --git a/ci_support/ensure_sundials_installed.sh b/ci_support/ensure_sundials_installed.sh index 8c60e83..33a4dbd 100755 --- a/ci_support/ensure_sundials_installed.sh +++ b/ci_support/ensure_sundials_installed.sh @@ -1,9 +1,10 @@ #!/bin/sh -SUNDIALS_DEFAULT_VERSION='2.7.0' +SUNDIALS_DEFAULT_VERSION='3.1.1' SUNDIALS_DEFAULT_PRECISION='double' +SUNDIALS_DEFAULT_INDEX_TYPE='int32_t' -export SUNDIALS_DIR=$HOME/sundials/"${SUNDIALS_VERSION:-$SUNDIALS_DEFAULT_VERSION}"/"${SUNDIALS_PRECISION:-$SUNDIALS_DEFAULT_PRECISION}" +export SUNDIALS_DIR=$HOME/sundials/"${SUNDIALS_VERSION:-$SUNDIALS_DEFAULT_VERSION}"/"${SUNDIALS_PRECISION:-$SUNDIALS_DEFAULT_PRECISION}"/"${SUNDIALS_INDEX_TYPE:-$SUNDIALS_DEFAULT_INDEX_TYPE}" SUNDIALS_LIBDIR=$SUNDIALS_DIR/lib SUNDIALS_INCLUDEDIR=$SUNDIALS_DIR/include diff --git a/ci_support/install_sundials.sh b/ci_support/install_sundials.sh index 9de7c6e..68b39f9 100755 --- a/ci_support/install_sundials.sh +++ b/ci_support/install_sundials.sh @@ -1,10 +1,11 @@ #!/bin/sh set -ex -SUNDIALS=sundials-"${SUNDIALS_VERSION:-2.7.0}" +SUNDIALS=sundials-"${SUNDIALS_VERSION:-3.1.1}" SUNDIALS_FILE=$SUNDIALS.tar.gz SUNDIALS_URL=http://computation.llnl.gov/projects/sundials-suite-nonlinear-differential-algebraic-equation-solvers/download/$SUNDIALS_FILE PRECISION="${SUNDIALS_PRECISION:-double}" +INDEX_TYPE="${SUNDIALS_INDEX_TYPE:-int32_t}" wget "$SUNDIALS_URL" @@ -16,5 +17,5 @@ mkdir sundials_build # ARKODE and examples are disabled as they does not support multiple precisions # (should be fixed in the 3.0 release) cd sundials_build && - cmake -DCMAKE_INSTALL_PREFIX=$SUNDIALS_DIR -DLAPACK_ENABLE=ON -DSUNDIALS_PRECISION="$PRECISION" -DBUILD_ARKODE:BOOL=OFF -DEXAMPLES_ENABLE:BOOL=OFF -DEXAMPLES_INSTALL:BOOL=OFF ../$SUNDIALS && + cmake -DCMAKE_INSTALL_PREFIX=$SUNDIALS_DIR -DLAPACK_ENABLE=ON -DSUNDIALS_INDEX_TYPE="$INDEX_TYPE" -DSUNDIALS_PRECISION="$PRECISION" -DBUILD_ARKODE:BOOL=OFF -DEXAMPLES_ENABLE:BOOL=OFF -DEXAMPLES_INSTALL:BOOL=OFF ../$SUNDIALS && make && make install diff --git a/common.py b/common.py index bbd1af1..c5ebcb6 100644 --- a/common.py +++ b/common.py @@ -21,8 +21,8 @@ BUILD_REQUIRES = ['numpy', 'cython'] # This is need for older pip MAJOR = 2 -MINOR = 3 -MICRO = 3 +MINOR = 4 +MICRO = 0 DEV = True CLASSIFIERS = [ @@ -35,10 +35,10 @@ 'Programming Language :: Python :: 2', 'Programming Language :: Python :: 2.7', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.2', - 'Programming Language :: Python :: 3.3', 'Programming Language :: Python :: 3.4', 'Programming Language :: Python :: 3.5', + 'Programming Language :: Python :: 3.6', + 'Programming Language :: Python :: 3.7', ] def build_verstring(): diff --git a/docs/examples/doublependulum.py b/docs/examples/doublependulum.py index 6a03347..4fd2887 100644 --- a/docs/examples/doublependulum.py +++ b/docs/examples/doublependulum.py @@ -168,7 +168,7 @@ def ddaspk_res(self, tres, yy, yp, res): def ddaspk_jac(self, tres, yy, yp, cj, jac): """the jacobian function as required by ddaspk""" - self.jac.evaluate(tres, yy, yp, cj, jac) + self.jac.evaluate(tres, yy, yp, None, cj, jac) #classes for the equations, as needed for the chosen solution method class resindex2(ida.IDA_RhsFunction): @@ -242,7 +242,7 @@ def set_dblpend(self, dblpend): the data """ self.dblpend = dblpend - def evaluate(self, tres, yy, yp, cj, jac): + def evaluate(self, tres, yy, yp, residual, cj, jac): m1 = self.dblpend.m1 m2 = self.dblpend.m2 diff --git a/docs/examples/ode/simpleoscillator_jac.py b/docs/examples/ode/simpleoscillator_jac.py index 2e71184..28ba7c3 100644 --- a/docs/examples/ode/simpleoscillator_jac.py +++ b/docs/examples/ode/simpleoscillator_jac.py @@ -1,12 +1,12 @@ -# Authors: B. Malengier +# Authors: B. Malengier """ -This example shows the most simple way of using a solver. +This example shows the most simple way of using a solver. We solve free vibration of a simple oscillator:: m \ddot{u} + k u = 0, u(0) = u_0, \dot{u}(0) = \dot{u}_0 using the CVODE solver, which means we use a rhs function of \dot{u}. Solution:: u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m) - + """ from __future__ import print_function from numpy import asarray, cos, sin, sqrt @@ -22,11 +22,11 @@ def rhseqn(t, x, xdot): """ we create rhs equations for the problem""" xdot[0] = x[1] xdot[1] = - k/m * x[0] - -def jaceqn(t, x, jac): + +def jaceqn(t, x, fx, jac): jac[0,1] = 1 jac[1,0] = -k/m - + #instantiate the solver from scikits.odes import ode solver = ode('cvode', rhseqn, jacfn=jaceqn) diff --git a/docs/examples/ode/simpleoscillator_prec.py b/docs/examples/ode/simpleoscillator_prec.py index c3a9cd2..61543ae 100644 --- a/docs/examples/ode/simpleoscillator_prec.py +++ b/docs/examples/ode/simpleoscillator_prec.py @@ -7,11 +7,11 @@ using the CVODE solver. The rhs function is given by \dot{u}. The preconditioning is implemented in prec_solvefn which solves the system P = I - gamma * J, where J is the Jacobian of the rhs. -The jac_times_vecfn calculates the Jacobian times vector product -using the analytic Jacobian. +The jac_times_vecfn calculates the Jacobian times vector product +using the analytic Jacobian. Solution:: u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m) - + """ from __future__ import print_function from numpy import asarray, cos, sin, sqrt @@ -32,17 +32,17 @@ def prec_solvefn(t, y, r, z, gamma, delta, lr): def jac_times_vecfn(v, Jv, t, y, user_data): """ Calculate Jacobian times vector product Jv = J*v""" Jv[0] = v[1] - Jv[1] = -k/m * v[0] + Jv[1] = -k/m * v[0] #define function for the right-hand-side equations which has specific signature def rhseqn(t, x, xdot): """ we create rhs equations for the problem""" xdot[0] = x[1] xdot[1] = - k/m * x[0] - + #instantiate the solver using a left-preconditioned BiCGStab as linear solver from scikits.odes import ode -solver = ode('cvode', rhseqn, linsolver='spbcg', precond_type='left', prec_solvefn=prec_solvefn, jac_times_vecfn=jac_times_vecfn) +solver = ode('cvode', rhseqn, linsolver='spbcgs', precond_type='left', prec_solvefn=prec_solvefn, jac_times_vecfn=jac_times_vecfn) #obtain solution at a required time result = solver.solve([0., 1., 2.], initx) diff --git a/docs/installation.rst b/docs/installation.rst index 7d8d442..9a36164 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -10,17 +10,17 @@ Before building ``odes``, you need to have installed: distributions, ``python-devel`` on Fedora) * C compiler * Fortran compiler (e.g. gfortran) - * `Sundials 2.7.0 `_ + * `Sundials 3.1.1 `_ In addition, if building from a git checkout, you'll also need Cython. It is required that Sundials is built with the BLAS/LAPACK interface enabled, so check the Fortran Settings section. A typical install if sundials download package is -extracted into directory sundials-2.7.0 is on a \*nix system:: +extracted into directory sundials-3.1.1 is on a \*nix system:: - mkdir build-sundials-2.7.0 - cd build-sundials-2.7.0/ - cmake -DLAPACK_ENABLE=ON -DCMAKE_INSTALL_PREFIX= ../sundials-2.7.0/ + mkdir build-sundials-3.1.1 + cd build-sundials-3.1.1/ + cmake -DLAPACK_ENABLE=ON -DSUNDIALS_INDEX_TYPE=int32_t -DCMAKE_INSTALL_PREFIX= ../sundials-3.1.1/ make install .. warning:: @@ -43,6 +43,15 @@ which will download the latest version from PyPI. This will handle the installat If you have installed SUNDIALS in a non-standard path (e.g. ``/usr/`` or ``/usr/local/``), you can set ``$SUNDIALS_INST`` in your environment to the installation prefix of SUNDIALS (i.e. value of ```` mentioned above). + +Testing your version of ``odes`` +................................ +To test the version in python, use in the python shell:: + + >>> import pkg_resources + >>> pkg_resources.get_distribution("scikits.odes").version + + Running the Tests ................. You need nose to run the tests. To install nose, run:: @@ -52,6 +61,10 @@ You need nose to run the tests. To install nose, run:: To run the tests, in the python shell:: >>> import scikits.odes as od; od.test() + +Note that the sundials library must be in your ``LD_LIBRARY_PATH``. So, make sure the directory ``$SUNDIALS_INST/lib`` is included. You can do this for example as follows (assuming sundials was installed in ``/usr/local``:: + + export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH Installation of ODES from git checkout --------------------------------------------- diff --git a/scikits/odes/dae.py b/scikits/odes/dae.py index 48d0a0d..896f56c 100644 --- a/scikits/odes/dae.py +++ b/scikits/odes/dae.py @@ -442,7 +442,8 @@ def __del__(self): """ Clean up what is needed """ - del self._integrator + if hasattr(self, '_integrator'): + del self._integrator #------------------------------------------------------------------------------ # DAE integrators diff --git a/scikits/odes/sundials/c_cvode.pxd b/scikits/odes/sundials/c_cvode.pxd index 98d863a..2a1eade 100644 --- a/scikits/odes/sundials/c_cvode.pxd +++ b/scikits/odes/sundials/c_cvode.pxd @@ -127,29 +127,24 @@ cdef extern from "cvode/cvode_direct.h": enum: CVDLS_JACFUNC_UNRECVR # -5 enum: CVDLS_JACFUNC_RECVR # -6 + enum: CVDLS_SUNMAT_FAIL # -7 - ctypedef int (*CVDlsDenseJacFn)(long int N, realtype t, - N_Vector y, N_Vector fy, - DlsMat Jac, void *user_data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except? -1 - ctypedef int (*CVDlsBandJacFn)(long int N, long int mupper, long int mlower, - 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 jac) - int CVDlsSetBandJacFn(void *cvode_mem, CVDlsBandJacFn jac) + ctypedef int (*CVDlsJacFn)(realtype t, N_Vector y, N_Vector fy, + SUNMatrix Jac, void *user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except? -1 + + int CVDlsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS, + SUNMatrix A) + int CVDlsSetJacFn(void *cvode_mem, CVDlsJacFn jac) int CVDlsGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS) int CVDlsGetNumJacEvals(void *cvode_mem, long int *njevals) int CVDlsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) int CVDlsGetLastFlag(void *cvode_mem, long int *flag) char *CVDlsGetReturnFlagName(long int flag) -cdef extern from "cvode/cvode_band.h": - int CVBand(void *cvode_mem, long int N, long int mupper, long int mlower) - cdef extern from "cvode/cvode_bandpre.h": - int CVBandPrecInit(void *cvode_mem, long int N, long int mu, long int ml) + int CVBandPrecInit(void *cvode_mem, sunindextype N, sunindextype mu, + sunindextype ml); int CVBandPrecGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS) int CVBandPrecGetNumRhsEvals(void *cvode_mem, long int *nfevalsBP) @@ -172,28 +167,21 @@ cdef extern from "cvode/cvode_diag.h": enum: CVDIAG_RHSFUNC_RECVR # -7 cdef extern from "cvode/cvode_bbdpre.h": - ctypedef int (*CVLocalFn)(long int Nlocal, realtype t, N_Vector y, + ctypedef int (*CVLocalFn)(sunindextype Nlocal, realtype t, N_Vector y, N_Vector g, void *user_data) - ctypedef int (*CVCommFn)(long int Nlocal, realtype t, N_Vector y, + ctypedef int (*CVCommFn)(sunindextype Nlocal, realtype t, N_Vector y, void *user_data) - int CVBBDPrecInit(void *cvode_mem, long int Nlocal, - long int mudq, long int mldq, - long int mukeep, long int mlkeep, + int CVBBDPrecInit(void *cvode_mem, sunindextype Nlocal, + sunindextype mudq, sunindextype mldq, + sunindextype mukeep, sunindextype mlkeep, realtype dqrely, CVLocalFn gloc, CVCommFn cfn) - int CVBBDPrecReInit(void *cvode_mem, long int mudq, long int mldq, + int CVBBDPrecReInit(void *cvode_mem, sunindextype mudq, sunindextype mldq, realtype dqrely) int CVBBDPrecGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS) int CVBBDPrecGetNumGfnEvals(void *cvode_mem, long int *ngevalsBBDP) -cdef extern from "cvode/cvode_dense.h": - int CVDense(void *cvode_mem, long int N) - -cdef extern from "cvode/cvode_lapack.h": - int CVLapackDense(void *cvode_mem, int N) - int CVLapackBand(void *cvode_mem, int N, int mupper, int mlower) - cdef extern from "cvode/cvode_spils.h": # CVSPILS return values enum: CVSPILS_SUCCESS # 0 @@ -202,50 +190,41 @@ cdef extern from "cvode/cvode_spils.h": enum: CVSPILS_ILL_INPUT # -3 enum: CVSPILS_MEM_FAIL # -4 enum: CVSPILS_PMEM_NULL # -5 + enum: CVSPILS_SUNLS_FAIL # -6 - enum: CVSPILS_MAXL # 5 enum: CVSPILS_MSBPRE # 50 enum: CVSPILS_DGMAX # RCONST(0.2) enum: CVSPILS_EPLIN # RCONST(0.05) ctypedef int (*CVSpilsPrecSetupFn)(realtype t, N_Vector y, N_Vector fy, booleantype jok, booleantype *jcurPtr, - realtype gamma, void *user_data, - N_Vector tmp1, N_Vector tmp2, - N_Vector tmp3) except? -1 + realtype gamma, void *user_data) except? -1 ctypedef int (*CVSpilsPrecSolveFn)(realtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, realtype gamma, realtype delta, - int lr, void *user_data, N_Vector tmp) except? -1 + int lr, void *user_data) except? -1 + ctypedef int (*CVSpilsJacTimesSetupFn)(realtype t, N_Vector y, + N_Vector fy, void *user_data) except? -1 ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t, N_Vector y, N_Vector fy, void *user_data, N_Vector tmp) except? -1 - int CVSpilsSetPrecType(void *cvode_mem, int pretype) - int CVSpilsSetGSType(void *cvode_mem, int gstype) - int CVSpilsSetMaxl(void *cvode_mem, int maxl) + int CVSpilsSetLinearSolver(void *cvode_mem, SUNLinearSolver LS) int CVSpilsSetEpsLin(void *cvode_mem, realtype eplifac) int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn pset, CVSpilsPrecSolveFn psolve) - int CVSpilsSetJacTimesVecFn(void *cvode_mem, - CVSpilsJacTimesVecFn jtv) + int CVSpilsSetJacTimes(void *cvode_mem, + CVSpilsJacTimesSetupFn jtsetup, + CVSpilsJacTimesVecFn jtimes) int CVSpilsGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS) int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals) int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves) int CVSpilsGetNumLinIters(void *cvode_mem, long int *nliters) int CVSpilsGetNumConvFails(void *cvode_mem, long int *nlcfails) + int CVSpilsGetNumJTSetupEvals(void *cvode_mem, long int *njtsetups) int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) int CVSpilsGetLastFlag(void *cvode_mem, long int *flag) char *CVSpilsGetReturnFlagName(long int flag) - -cdef extern from "cvode/cvode_spgmr.h": - int CVSpgmr(void *cvode_mem, int pretype, int maxl) - -cdef extern from "cvode/cvode_spbcgs.h": - int CVSpbcg(void *cvode_mem, int pretype, int maxl) - -cdef extern from "cvode/cvode_sptfqmr.h": - int CVSptfqmr(void *cvode_mem, int pretype, int maxl) diff --git a/scikits/odes/sundials/c_ida.pxd b/scikits/odes/sundials/c_ida.pxd index 3d6fd09..996ec52 100644 --- a/scikits/odes/sundials/c_ida.pxd +++ b/scikits/odes/sundials/c_ida.pxd @@ -45,6 +45,7 @@ cdef extern from "ida/ida.h": enum: IDA_BAD_K #-25 enum: IDA_BAD_T #-26 enum: IDA_BAD_DKY #-27 + enum: IDA_UNRECOGNISED_ERROR #-99 ctypedef int (*IDAResFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data) @@ -144,74 +145,52 @@ cdef extern from "ida/ida_direct.h": #/* Additional last_flag values */ enum: IDADLS_JACFUNC_UNRECVR enum: IDADLS_JACFUNC_RECVR - - ctypedef int (*IDADlsDenseJacFn)(long int N, realtype t, realtype c_j, - N_Vector y, N_Vector yp, N_Vector r, - DlsMat Jac, void *user_data, N_Vector tmp1, - N_Vector tmp2, N_Vector tmp3) - - ctypedef int (*IDADlsBandJacFn)(long int N, long int mupper, long int mlower, - realtype t, realtype c_j, - N_Vector y, N_Vector yp, N_Vector r, - DlsMat Jac, void *user_data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) - ctypedef int (*IDADlsBandJacFn)(int N, int mupper, int mlower, - realtype t, realtype c_j, - N_Vector y, N_Vector yp, N_Vector r, - DlsMat Jac, void *user_data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) - int IDADlsSetDenseJacFn(void *ida_mem, IDADlsDenseJacFn jac) - int IDADlsSetBandJacFn(void *ida_mem, IDADlsBandJacFn jac) + + ctypedef int (*IDADlsJacFn)(realtype t, realtype c_j, N_Vector y, + N_Vector yp, N_Vector r, SUNMatrix Jac, + void *user_data, N_Vector tmp1, + N_Vector tmp2, N_Vector tmp3) + int IDADlsSetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix A) + int IDADlsSetJacFn(void *ida_mem, IDADlsJacFn jac) + int IDADlsGetWorkSpace(void *ida_mem, long int *lenrwLS, long int *leniwLS) int IDADlsGetNumJacEvals(void *ida_mem, long int *njevals) int IDADlsGetNumResEvals(void *ida_mem, long int *nfevalsLS) int IDADlsGetLastFlag(void *ida_mem, long int *flag) char *IDADlsGetReturnFlagName(long int flag) -cdef extern from "ida/ida_band.h": - int IDABand(void *ida_mem, long int Neq, long int mupper, long int mlower) - -cdef extern from "ida/ida_dense.h": - int IDADense(void *ida_mem, long int Neq) - -cdef extern from "ida/ida_lapack.h": - int IDALapackDense(void *ida_mem, int N) - int IDALapackBand(void *ida_mem, int N, int mupper, int mlower) - cdef extern from "ida/ida_spils.h": - enum: IDASPILS_SUCCESS # 0 enum: IDASPILS_MEM_NULL # -1 enum: IDASPILS_LMEM_NULL # -2 enum: IDASPILS_ILL_INPUT # -3 enum: IDASPILS_MEM_FAIL # -4 enum: IDASPILS_PMEM_NULL # -5 + enum: IDASPILS_SUNLS_FAIL # -6 ctypedef int (*IDASpilsPrecSetupFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, - realtype c_j, void *user_data, - N_Vector tmp1, N_Vector tmp2, - N_Vector tmp3) + realtype c_j, void *user_data) ctypedef int (*IDASpilsPrecSolveFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector rvec, N_Vector zvec, - realtype c_j, realtype delta, void *user_data, - N_Vector tmp) + realtype c_j, realtype delta, void *user_data) + ctypedef int (*IDASpilsJacTimesSetupFn)(realtype tt, N_Vector yy, + N_Vector yp, N_Vector rr, + realtype c_j, void *user_data) ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, realtype c_j, void *user_data, N_Vector tmp1, N_Vector tmp2) + + int IDASpilsSetLinearSolver(void *ida_mem, SUNLinearSolver LS) int IDASpilsSetPreconditioner(void *ida_mem, IDASpilsPrecSetupFn pset, IDASpilsPrecSolveFn psolve) - int IDASpilsSetJacTimesVecFn(void *ida_mem, - IDASpilsJacTimesVecFn jtv) - - int IDASpilsSetGSType(void *ida_mem, int gstype) - int IDASpilsSetMaxRestarts(void *ida_mem, int maxrs) - int IDASpilsSetMaxl(void *ida_mem, int maxl) + int IDASpilsSetJacTimes(void *ida_mem, IDASpilsJacTimesSetupFn jtsetup, + IDASpilsJacTimesVecFn jtimes) int IDASpilsSetEpsLin(void *ida_mem, realtype eplifac) int IDASpilsSetIncrementFactor(void *ida_mem, realtype dqincfac) int IDASpilsGetWorkSpace(void *ida_mem, long int *lenrwLS, long int *leniwLS) @@ -224,30 +203,21 @@ cdef extern from "ida/ida_spils.h": int IDASpilsGetLastFlag(void *ida_mem, long int *flag) char *IDASpilsGetReturnFlagName(long int flag) -cdef extern from "ida/ida_spgmr.h": - int IDASpgmr(void *ida_mem, int maxl) - -cdef extern from "ida/ida_spbcgs.h": - int IDASpbcg(void *ida_mem, int maxl) - -cdef extern from "ida/ida_sptfqmr.h": - int IDASptfqmr(void *ida_mem, int maxl) - cdef extern from "ida/ida_bbdpre.h": - ctypedef int (*IDABBDLocalFn)(long int Nlocal, realtype tt, + ctypedef int (*IDABBDLocalFn)(sunindextype Nlocal, realtype tt, N_Vector yy, N_Vector yp, N_Vector gval, void *user_data) - ctypedef int (*IDABBDCommFn)(long int Nlocal, realtype tt, + ctypedef int (*IDABBDCommFn)(sunindextype Nlocal, realtype tt, N_Vector yy, N_Vector yp, void *user_data) - int IDABBDPrecInit(void *ida_mem, long int Nlocal, - long int mudq, long int mldq, - long int mukeep, long int mlkeep, + int IDABBDPrecInit(void *ida_mem, sunindextype Nlocal, + sunindextype mudq, sunindextype mldq, + sunindextype mukeep, sunindextype mlkeep, realtype dq_rel_yy, IDABBDLocalFn Gres, IDABBDCommFn Gcomm) int IDABBDPrecReInit(void *ida_mem, - long int mudq, long int mldq, + sunindextype mudq, sunindextype mldq, realtype dq_rel_yy) int IDABBDPrecGetWorkSpace(void *ida_mem, long int *lenrwBBDP, long int *leniwBBDP) diff --git a/scikits/odes/sundials/c_nvector_serial.pxd b/scikits/odes/sundials/c_nvector_serial.pxd index e3df75c..40779d2 100644 --- a/scikits/odes/sundials/c_nvector_serial.pxd +++ b/scikits/odes/sundials/c_nvector_serial.pxd @@ -1,24 +1,29 @@ -include "c_sundials.pxd" +from libc.stdio cimport FILE +from .c_sundials cimport * cdef extern from "nvector/nvector_serial.h": cdef struct _N_VectorContent_Serial: - long int length + sunindextype length booleantype own_data realtype *data ctypedef _N_VectorContent_Serial *N_VectorContent_Serial - N_Vector N_VNew_Serial(long int vec_length) - N_Vector N_VNewEmpty_Serial(long int vec_length) - N_Vector N_VMake_Serial(long int vec_length, realtype *v_data) + N_Vector N_VNew_Serial(sunindextype vec_length) + N_Vector N_VNewEmpty_Serial(sunindextype vec_length) + N_Vector N_VMake_Serial(sunindextype 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_VDestroyVectorArray_Serial(N_Vector *vs, int count) + sunindextype N_VGetLength_Serial(N_Vector v) void N_VPrint_Serial(N_Vector v) + void N_VPrintFile_Serial(N_Vector v, FILE *outfile) + + N_Vector_ID N_VGetVectorID_Serial(N_Vector v) N_Vector N_VCloneEmpty_Serial(N_Vector w) N_Vector N_VClone_Serial(N_Vector w) void N_VDestroy_Serial(N_Vector v) - void N_VSpace_Serial(N_Vector v, long int *lrw, long int *liw) + void N_VSpace_Serial(N_Vector v, sunindextype *lrw, sunindextype *liw) realtype *N_VGetArrayPointer_Serial(N_Vector v) void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v) void N_VLinearSum_Serial(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z) @@ -40,4 +45,7 @@ cdef extern from "nvector/nvector_serial.h": booleantype N_VInvTest_Serial(N_Vector x, N_Vector z) booleantype N_VConstrMask_Serial(N_Vector c, N_Vector x, N_Vector m) realtype N_VMinQuotient_Serial(N_Vector num, N_Vector denom) + # Macros + sunindextype NV_LENGTH_S(N_Vector vc_s) + realtype* NV_DATA_S(N_Vector vc_s) diff --git a/scikits/odes/sundials/c_sundials.pxd b/scikits/odes/sundials/c_sundials.pxd index d9dd93f..cd190ac 100644 --- a/scikits/odes/sundials/c_sundials.pxd +++ b/scikits/odes/sundials/c_sundials.pxd @@ -1,8 +1,22 @@ +from libc.stdio cimport FILE + cdef extern from "sundials/sundials_types.h": ctypedef float realtype ctypedef unsigned int booleantype + ctypedef long sunindextype cdef extern from "sundials/sundials_nvector.h": + cdef enum N_Vector_ID: + SUNDIALS_NVEC_SERIAL, + SUNDIALS_NVEC_PARALLEL, + SUNDIALS_NVEC_OPENMP, + SUNDIALS_NVEC_PTHREADS, + SUNDIALS_NVEC_PARHYP, + SUNDIALS_NVEC_PETSC, + SUNDIALS_NVEC_CUDA, + SUNDIALS_NVEC_RAJA, + SUNDIALS_NVEC_CUSTOM + struct _generic_N_Vector_Ops: pass struct _generic_N_Vector: @@ -11,10 +25,11 @@ cdef extern from "sundials/sundials_nvector.h": ctypedef _generic_N_Vector_Ops *N_Vector_Ops struct _generic_N_Vector_Ops: + N_Vector_ID (*nvgetvectorid)(N_Vector) N_Vector (*nvclone)(N_Vector) N_Vector (*nvcloneempty)(N_Vector) void (*nvdestroy)(N_Vector) - void (*nvspace)(N_Vector, long int *, long int *) + void (*nvspace)(N_Vector, sunindextype *, sunindextype *) realtype* (*nvgetarraypointer)(N_Vector) void (*nvsetarraypointer)(realtype *, N_Vector) void (*nvlinearsum)(realtype, N_Vector, realtype, N_Vector, N_Vector) @@ -38,10 +53,11 @@ cdef extern from "sundials/sundials_nvector.h": realtype (*nvminquotient)(N_Vector, N_Vector) # * FUNCTIONS * + N_Vector_ID N_VGetVectorID(N_Vector w) N_Vector N_VClone(N_Vector w) N_Vector N_VCloneEmpty(N_Vector w) void N_VDestroy(N_Vector v) - void N_VSpace(N_Vector v, long int *lrw, long int *liw) + void N_VSpace(N_Vector v, sunindextype *lrw, sunindextype *liw) realtype *N_VGetArrayPointer(N_Vector v) void N_VSetArrayPointer(realtype *v_data, N_Vector v) void N_VLinearSum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z) @@ -68,82 +84,130 @@ cdef extern from "sundials/sundials_nvector.h": N_Vector *N_VCloneVectorArray(int count, N_Vector w) void N_VDestroyVectorArray(N_Vector *vs, int count) -cdef extern from "nvector/nvector_serial.h": - N_Vector N_VNew_Serial(long int vec_length) - N_Vector N_VNewEmpty_Serial(long int vec_length) - 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_VDestroyVectorArray_Serial(N_Vector *vs, int count) - - void N_VPrint_Serial(N_Vector v) - N_Vector N_VCloneEmpty_Serial(N_Vector w) - N_Vector N_VClone_Serial(N_Vector w) - void N_VDestroy_Serial(N_Vector v) - void N_VSpace_Serial(N_Vector v, long int *lrw, long int *liw) - realtype *N_VGetArrayPointer_Serial(N_Vector v) - void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v) - void N_VLinearSum_Serial(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z) - void N_VConst_Serial(realtype c, N_Vector z) - void N_VProd_Serial(N_Vector x, N_Vector y, N_Vector z) - void N_VDiv_Serial(N_Vector x, N_Vector y, N_Vector z) - void N_VScale_Serial(realtype c, N_Vector x, N_Vector z) - void N_VAbs_Serial(N_Vector x, N_Vector z) - void N_VInv_Serial(N_Vector x, N_Vector z) - void N_VAddConst_Serial(N_Vector x, realtype b, N_Vector z) - realtype N_VDotProd_Serial(N_Vector x, N_Vector y) - realtype N_VMaxNorm_Serial(N_Vector x) - realtype N_VWrmsNorm_Serial(N_Vector x, N_Vector w) - realtype N_VWrmsNormMask_Serial(N_Vector x, N_Vector w, N_Vector id) - realtype N_VMin_Serial(N_Vector x) - realtype N_VWL2Norm_Serial(N_Vector x, N_Vector w) - realtype N_VL1Norm_Serial(N_Vector x) - void N_VCompare_Serial(realtype c, N_Vector x, N_Vector z) - booleantype N_VInvTest_Serial(N_Vector x, N_Vector z) - booleantype N_VConstrMask_Serial(N_Vector c, N_Vector x, N_Vector m) - realtype N_VMinQuotient_Serial(N_Vector num, N_Vector denom) - # Macros - long int NV_LENGTH_S(N_Vector vc_s) - realtype* NV_DATA_S(N_Vector vc_s) - -cdef extern from "sundials/sundials_lapack.h": - pass - # void dcopy_(int *n, double *x, int *inc_x, double *y, int *inc_y) - # void dscal_(int *n, double *alpha, double *x, int *inc_x) - - # # Level-2 BLAS - # void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, - # int *lda, double *x, int *inc_x, double *beta, double *y, - # int *inc_y, int len_trans) - # void dtrsv_(char *uplo, char *trans, char *diag, int *n, - # double *a, int *lda, double *x, int *inc_x, - # int len_uplo, int len_trans, int len_diag) - - # # Level-3 BLAS - # void dsyrk_(char *uplo, char *trans, int *n, int *k, - # double *alpha, double *a, int *lda, double *beta, - # double *c, int *ldc, int len_uplo, int len_trans) - - # # LAPACK - # void dgbtrf_(int *m, int *n, int *kl, int *ku, - # double *ab, int *ldab, int *ipiv, int *info) - # void dgbtrs_(char *trans, int *n, int *kl, int *ku, int *nrhs, - # double *ab, int *ldab, int *ipiv, double *b, int *ldb, - # int *info, int len_trans) - # void dgeqp3_(int *m, int *n, double *a, int *lda, int *jpvt, double *tau, - # double *work, int *lwork, int *info) - # void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau, double *work, - # int *lwork, int *info) - # void dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info) - # void dgetrs_(char *trans, int *n, int *nrhs, double *a, int *lda, - # int *ipiv, double *b, int *ldb, int *info, int len_trans) - - # void dormqr_(char *side, char *trans, int *m, int *n, int *k, double *a, - # int *lda, double *tau, double *c, int *ldc, double *work, - # int *lwork, int *info, int len_side, int len_trans) - # void dpotrf_(char *uplo, int *n, double *a, int *lda, int *info, int len_uplo) - # void dpotrs_(char *uplo, int *n, int *nrhs, double *a, int *lda, - # double *b, int *ldb, int * info, int len_uplo) +cdef extern from "sundials/sundials_matrix.h": + cdef enum SUNMatrix_ID: + SUNMATRIX_DENSE, + SUNMATRIX_BAND, + SUNMATRIX_SPARSE, + SUNMATRIX_CUSTOM + + struct _generic_SUNMatrix_Ops: + pass + struct _generic_SUNMatrix: + pass + ctypedef _generic_SUNMatrix *SUNMatrix + ctypedef _generic_SUNMatrix_Ops *SUNMatrix_Ops + + 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*) + + #struct _generic_SUNMatrix: + # void *content + # struct _generic_SUNMatrix_Ops *ops + + # * FUNCTIONS * + SUNMatrix_ID SUNMatGetID(SUNMatrix A) + SUNMatrix SUNMatClone(SUNMatrix A) + void SUNMatDestroy(SUNMatrix A) + int SUNMatZero(SUNMatrix A) + int SUNMatCopy(SUNMatrix A, SUNMatrix B) + int SUNMatScaleAdd(realtype c, SUNMatrix A, SUNMatrix B) + int SUNMatScaleAddI(realtype c, SUNMatrix A) + int SUNMatMatvec(SUNMatrix A, N_Vector x, N_Vector y) + int SUNMatSpace(SUNMatrix A, long int *lenrw, long int *leniw) + +cdef extern from "sundials/sundials_iterative.h": + enum: + PREC_NONE + PREC_LEFT + PREC_RIGHT + PREC_BOTH + enum: + MODIFIED_GS = 1 + CLASSICAL_GS = 2 + ctypedef int (*ATimesFn)(void *A_data, N_Vector v, N_Vector z) + ctypedef int (*PSetupFn)(void *P_data) + ctypedef int (*PSolveFn)(void *P_data, N_Vector r, N_Vector z, + realtype tol, int lr) + + int ModifiedGS(N_Vector *v, realtype **h, int k, int p, + realtype *new_vk_norm) + int ClassicalGS(N_Vector *v, realtype **h, int k, int p, + realtype *new_vk_norm, N_Vector temp, realtype *s) + int QRfact(int n, realtype **h, realtype *q, int job) + int QRsol(int n, realtype **h, realtype *q, realtype *b) + +# We don't support KLU for now +#cdef extern from "sundials/sundials_klu_impl.h": +# cdef struct KLUDataRec: +# klu_symbolic *s_Symbolic +# klu_numeric *s_Numeric +# klu_common s_Common +# int s_ordering +# ctypedef KLUDataRec *KLUData + +cdef extern from "sundials/sundials_linearsolver.h": + + cdef enum SUNLinearSolver_Type: + SUNLINEARSOLVER_DIRECT, + SUNLINEARSOLVER_ITERATIVE, + SUNLINEARSOLVER_CUSTOM + + struct _generic_SUNLinearSolver_Ops: + pass + struct _generic_SUNLinearSolver: + pass + ctypedef _generic_SUNLinearSolver *SUNLinearSolver + ctypedef _generic_SUNLinearSolver_Ops *SUNLinearSolver_Ops + + 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) + + #struct _generic_SUNLinearSolver: + # void *content + # struct _generic_SUNLinearSolver_Ops *ops + + SUNLinearSolver_Type SUNLinSolGetType(SUNLinearSolver S) + int SUNLinSolSetATimes(SUNLinearSolver S, void* A_data, + ATimesFn ATimes) + int SUNLinSolSetPreconditioner(SUNLinearSolver S, void* P_data, + PSetupFn Pset, PSolveFn Psol) + int SUNLinSolSetScalingVectors(SUNLinearSolver S, N_Vector s1, + N_Vector s2) + int SUNLinSolInitialize(SUNLinearSolver S) + int SUNLinSolSetup(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve(SUNLinearSolver S, SUNMatrix A, N_Vector x, + N_Vector b, realtype tol) + int SUNLinSolNumIters(SUNLinearSolver S) + realtype SUNLinSolResNorm(SUNLinearSolver S) + N_Vector SUNLinSolResid(SUNLinearSolver S) + long int SUNLinSolLastFlag(SUNLinearSolver S) + int SUNLinSolSpace(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree(SUNLinearSolver S) cdef extern from "sundials/sundials_direct.h": @@ -152,111 +216,83 @@ cdef extern from "sundials/sundials_direct.h": cdef struct _DlsMat: int type - long int M - long int N - long int ldim - long int mu - long int ml - long int s_mu + sunindextype M + sunindextype N + sunindextype ldim + sunindextype mu + sunindextype ml + sunindextype s_mu realtype *data - long int ldata + sunindextype ldata realtype **cols ctypedef _DlsMat *DlsMat - DlsMat NewDenseMat(long int M, long int N) - DlsMat NewBandMat(long int N, long int mu, long int ml, long int smu) + DlsMat NewDenseMat(sunindextype M, sunindextype N) + DlsMat NewBandMat(sunindextype N, sunindextype mu, sunindextype ml, sunindextype smu) void DestroyMat(DlsMat A) int *NewIntArray(int N) - long int *NewLintArray(long int N) + sunindextype *NewIndexArray(sunindextype N) realtype *NewRealArray(long int N) void DestroyArray(void *p) void AddIdentity(DlsMat A) void SetToZero(DlsMat A) - void PrintMat(DlsMat A) + void PrintMat(DlsMat A, FILE *outfile) # * Exported function prototypes (functions working on realtype**) - realtype **newDenseMat(long int m, long int n) - realtype **newBandMat(long int n, long int smu, long int ml) + realtype **newDenseMat(sunindextype m, sunindextype n) + realtype **newBandMat(sunindextype n, sunindextype smu, sunindextype ml) void destroyMat(realtype **a) int *newIntArray(int n) - long int *newLintArray(long int n) - realtype *newRealArray(long int m) + sunindextype *newIndexArray(sunindextype n) + realtype *newRealArray(sunindextype m) void destroyArray(void *v) cdef extern from "sundials/sundials_band.h": - long int BandGBTRF(DlsMat A, long int *p) - long int bandGBTRF(realtype **a, long int n, long int mu, long int ml, - long int smu, long int *p) - void BandGBTRS(DlsMat A, long int *p, realtype *b) - void bandGBTRS(realtype **a, long int n, long int smu, long int ml, - long int *p, realtype *b) - void BandCopy(DlsMat A, DlsMat B, long int copymu, long int copyml) - void bandCopy(realtype **a, realtype **b, long int n, long int a_smu, - long int b_smu, long int copymu, long int copyml) + sunindextype BandGBTRF(DlsMat A, sunindextype *p) + sunindextype bandGBTRF(realtype **a, sunindextype n, sunindextype mu, sunindextype ml, + sunindextype smu, sunindextype *p) + void BandGBTRS(DlsMat A, sunindextype *p, realtype *b) + void bandGBTRS(realtype **a, sunindextype n, sunindextype smu, sunindextype ml, + sunindextype *p, realtype *b) + void BandCopy(DlsMat A, DlsMat B, sunindextype copymu, sunindextype copyml) + void bandCopy(realtype **a, realtype **b, sunindextype n, sunindextype a_smu, + sunindextype b_smu, sunindextype copymu, sunindextype copyml) void BandScale(realtype c, DlsMat A) - void bandScale(realtype c, realtype **a, long int n, long int mu, - long int ml, long int smu) - void bandAddIdentity(realtype **a, long int n, long int smu) + void bandScale(realtype c, realtype **a, sunindextype n, sunindextype mu, + sunindextype ml, sunindextype smu) + void bandAddIdentity(realtype **a, sunindextype n, sunindextype smu) void BandMatvec(DlsMat A, realtype *x, realtype *y) - void bandMatvec(realtype **a, realtype *x, realtype *y, long int n, - long int mu, long int ml, long int smu) + void bandMatvec(realtype **a, realtype *x, realtype *y, sunindextype n, + sunindextype mu, sunindextype ml, sunindextype smu) cdef extern from "sundials/sundials_dense.h": - long int DenseGETRF(DlsMat A, long int *p) - void DenseGETRS(DlsMat A, long int *p, realtype *b) + sunindextype DenseGETRF(DlsMat A, sunindextype *p) + void DenseGETRS(DlsMat A, sunindextype *p, realtype *b) - long int denseGETRF(realtype **a, long int m, long int n, long int *p) - void denseGETRS(realtype **a, long int n, long int *p, realtype *b) + sunindextype denseGETRF(realtype **a, sunindextype m, sunindextype n, sunindextype *p) + void denseGETRS(realtype **a, sunindextype n, sunindextype *p, realtype *b) - long int DensePOTRF(DlsMat A) + sunindextype DensePOTRF(DlsMat A) void DensePOTRS(DlsMat A, realtype *b) - long int densePOTRF(realtype **a, long int m) - void densePOTRS(realtype **a, long int m, realtype *b) + sunindextype densePOTRF(realtype **a, sunindextype m) + void densePOTRS(realtype **a, sunindextype m, realtype *b) int DenseGEQRF(DlsMat A, realtype *beta, realtype *wrk) int DenseORMQR(DlsMat A, realtype *beta, realtype *vn, realtype *vm, realtype *wrk) - int denseGEQRF(realtype **a, long int m, long int n, realtype *beta, realtype *v) - int denseORMQR(realtype **a, long int m, long int n, realtype *beta, + int denseGEQRF(realtype **a, sunindextype m, sunindextype n, realtype *beta, realtype *v) + int denseORMQR(realtype **a, sunindextype m, sunindextype n, realtype *beta, realtype *v, realtype *w, realtype *wrk) void DenseCopy(DlsMat A, DlsMat B) - void denseCopy(realtype **a, realtype **b, long int m, long int n) + void denseCopy(realtype **a, realtype **b, sunindextype m, sunindextype n) void DenseScale(realtype c, DlsMat A) - void denseScale(realtype c, realtype **a, long int m, long int n) - void denseAddIdentity(realtype **a, long int n) + void denseScale(realtype c, realtype **a, sunindextype m, sunindextype n) + void denseAddIdentity(realtype **a, sunindextype n) void DenseMatvec(DlsMat A, realtype *x, realtype *y) - void denseMatvec(realtype **a, realtype *x, realtype *y, long int m, long int n) - -cdef extern from "sundials/sundials_iterative.h": - enum: - PREC_NONE - PREC_LEFT - PREC_RIGHT - PREC_BOTH - enum: - MODIFIED_GS = 1 - CLASSICAL_GS = 2 - ctypedef int (*ATimesFn)(void *A_data, N_Vector v, N_Vector z) - ctypedef int (*PSolveFn)(void *P_data, N_Vector r, N_Vector z, int lr) - - int ModifiedGS(N_Vector *v, realtype **h, int k, int p, - realtype *new_vk_norm) - int ClassicalGS(N_Vector *v, realtype **h, int k, int p, - realtype *new_vk_norm, N_Vector temp, realtype *s) - int QRfact(int n, realtype **h, realtype *q, int job) - int QRsol(int n, realtype **h, realtype *q, realtype *b) - -# We don't support KLU for now -#cdef extern from "sundials/sundials_klu_impl.h": -# cdef struct KLUDataRec: -# klu_symbolic *s_Symbolic; -# klu_numeric *s_Numeric; -# klu_common s_Common; -# int s_ordering; -# ctypedef KLUDataRec *KLUData + void denseMatvec(realtype **a, realtype *x, realtype *y, sunindextype m, sunindextype n) cdef extern from "sundials/sundials_pcg.h": ctypedef struct _PcgMemRec: @@ -288,26 +324,36 @@ cdef extern from "sundials/sundials_pcg.h": enum: PCG_PSET_FAIL_UNREC # -4 /* pset failed unrecoverably */ cdef extern from "sundials/sundials_sparse.h": - ctypedef struct _SlsMat: + + enum: CSC_MAT #0 + enum: CSR_MAT #1 + + ctypedef struct _SlsMat: int M int N int NNZ + int NP realtype *data - int *rowvals - int *colptrs - ctypedef _SlsMat *SlsMat - - SlsMat NewSparseMat(int M, int N, int NNZ) - SlsMat SlsConvertDls(DlsMat A) - void DestroySparseMat(SlsMat A) - void SlsSetToZero(SlsMat A) - void CopySparseMat(SlsMat A, SlsMat B) - void ScaleSparseMat(realtype b, SlsMat A) - void AddIdentitySparseMat(SlsMat A) - int SlsAddMat(SlsMat A, SlsMat B) - void ReallocSparseMat(SlsMat A) - int SlsMatvec(SlsMat A, realtype *x, realtype *y) - void PrintSparseMat(SlsMat A) + int sparsetype + int *indexvals + int *indexptrs + int **rowvals + int **colptrs + int **colvals + int **rowptrs + ctypedef _SlsMat *SlsMat + + SlsMat SparseNewMat(int M, int N, int NNZ, int sparsetype) + SlsMat SparseFromDenseMat(const DlsMat A, int sparsetype) + int SparseDestroyMat(SlsMat A) + int SparseSetMatToZero(SlsMat A) + int SparseCopyMat(const SlsMat A, SlsMat B) + int SparseScaleMat(realtype b, SlsMat A) + int SparseAddIdentityMat(SlsMat A) + int SparseAddMat(SlsMat A, const SlsMat B) + int SparseReallocMat(SlsMat A) + int SparseMatvec(const SlsMat A, const realtype *x, realtype *y) + void SparsePrintMat(const SlsMat A, FILE* outfile) cdef extern from "sundials/sundials_spgmr.h": cdef struct _SpgmrMemRec: @@ -372,6 +418,7 @@ cdef extern from "sundials/sundials_spbcgs.h": enum: SPBCG_PSOLVE_FAIL_REC #3 /* psolve failed recoverably */ enum: SPBCG_ATIMES_FAIL_REC #4 /* atimes failed recoverably */ enum: SPBCG_PSET_FAIL_REC #5 /* pset faild recoverably */ + enum: SPBCG_MEM_NULL #-1 /* mem argument is NULL */ enum: SPBCG_ATIMES_FAIL_UNREC #-2 /* atimes returned failure flag */ enum: SPBCG_PSOLVE_FAIL_UNREC #-3 /* psolve failed unrecoverably */ @@ -452,7 +499,7 @@ cdef extern from "sundials/sundials_sptfqmr.h": void SptfqmrFree(SptfqmrMem mem) -# We don't use SuperLUMT +# We don't use SuperLUMT - "slu_mt_ddefs.h" required #cdef extern from "sundials/sundials_superlumt_impl.h": # ctypedef struct SLUMTDataRec: # SuperMatrix *s_A, *s_AC, *s_L, *s_U, *s_B diff --git a/scikits/odes/sundials/c_sunlinsol.pxd b/scikits/odes/sundials/c_sunlinsol.pxd new file mode 100644 index 0000000..fcb43bc --- /dev/null +++ b/scikits/odes/sundials/c_sunlinsol.pxd @@ -0,0 +1,391 @@ +from c_sundials cimport * +from c_sunmatrix cimport * + +include "sundials_config.pxi" + +cdef extern from "sunlinsol/sunlinsol_dense.h": + + struct _SUNLinearSolverContent_Dense: + sunindextype N + sunindextype *pivots + long int last_flag + + ctypedef _SUNLinearSolverContent_Dense *SUNLinearSolverContent_Dense + + SUNLinearSolver SUNDenseLinearSolver(N_Vector y, SUNMatrix A) + + SUNLinearSolver_Type SUNLinSolGetType_Dense(SUNLinearSolver S) + int SUNLinSolInitialize_Dense(SUNLinearSolver S) + int SUNLinSolSetup_Dense(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_Dense(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + long int SUNLinSolLastFlag_Dense(SUNLinearSolver S) + int SUNLinSolSpace_Dense(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_Dense(SUNLinearSolver S) + +cdef extern from "sunlinsol/sunlinsol_band.h": + + struct _SUNLinearSolverContent_Band: + sunindextype N + sunindextype *pivots + long int last_flag + + ctypedef _SUNLinearSolverContent_Band *SUNLinearSolverContent_Band + + SUNLinearSolver SUNBandLinearSolver(N_Vector y, SUNMatrix A) + + SUNLinearSolver_Type SUNLinSolGetType_Band(SUNLinearSolver S) + int SUNLinSolInitialize_Band(SUNLinearSolver S) + int SUNLinSolSetup_Band(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_Band(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + long int SUNLinSolLastFlag_Band(SUNLinearSolver S) + int SUNLinSolSpace_Band(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_Band(SUNLinearSolver S) + +# We don't support KLU for now +#cdef extern from "sunlinsol/sunlinsol_klu.h": + +IF SUNDIALS_BLAS_LAPACK: + cdef extern from "sunlinsol/sunlinsol_lapackdense.h": + + struct _SUNLinearSolverContent_LapackDense: + sunindextype N + sunindextype *pivots + long int last_flag + + ctypedef _SUNLinearSolverContent_LapackDense *SUNLinearSolverContent_LapackDense + + SUNLinearSolver SUNLapackDense(N_Vector y, SUNMatrix A) + + SUNLinearSolver_Type SUNLinSolGetType_LapackDense(SUNLinearSolver S) + int SUNLinSolInitialize_LapackDense(SUNLinearSolver S) + int SUNLinSolSetup_LapackDense(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_LapackDense(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + long int SUNLinSolLastFlag_LapackDense(SUNLinearSolver S) + int SUNLinSolSpace_LapackDense(SUNLinearSolver S, + long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_LapackDense(SUNLinearSolver S) + + + cdef extern from "sunlinsol/sunlinsol_lapackband.h": + + struct _SUNLinearSolverContent_LapackBand: + sunindextype N + sunindextype *pivots + long int last_flag + + ctypedef _SUNLinearSolverContent_LapackBand *SUNLinearSolverContent_LapackBand + + SUNLinearSolver SUNLapackBand(N_Vector y, SUNMatrix A) + + SUNLinearSolver_Type SUNLinSolGetType_LapackBand(SUNLinearSolver S) + int SUNLinSolInitialize_LapackBand(SUNLinearSolver S) + int SUNLinSolSetup_LapackBand(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_LapackBand(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + long int SUNLinSolLastFlag_LapackBand(SUNLinearSolver S) + int SUNLinSolSpace_LapackBand(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_LapackBand(SUNLinearSolver S) + + +cdef extern from "sunlinsol/sunlinsol_pcg.h": + + struct _SUNLinearSolverContent_PCG: + int maxl + int pretype + int numiters + realtype resnorm + long int last_flag + + ATimesFn ATimes + void* ATData + PSetupFn Psetup + PSolveFn Psolve + void* PData + + N_Vector s + N_Vector r + N_Vector p + N_Vector z + N_Vector Ap + + ctypedef _SUNLinearSolverContent_PCG *SUNLinearSolverContent_PCG + + SUNLinearSolver SUNPCG(N_Vector y, int pretype, int maxl) + int SUNPCGSetPrecType(SUNLinearSolver S, int pretype) + int SUNPCGSetMaxl(SUNLinearSolver S, int maxl) + + SUNLinearSolver_Type SUNLinSolGetType_PCG(SUNLinearSolver S) + int SUNLinSolInitialize_PCG(SUNLinearSolver S) + int SUNLinSolSetATimes_PCG(SUNLinearSolver S, void* A_data, ATimesFn ATimes) + int SUNLinSolSetPreconditioner_PCG(SUNLinearSolver S, void* P_data, + PSetupFn Pset, PSolveFn Psol) + int SUNLinSolSetScalingVectors_PCG(SUNLinearSolver S, N_Vector s, + N_Vector nul) + int SUNLinSolSetup_PCG(SUNLinearSolver S, SUNMatrix nul) + int SUNLinSolSolve_PCG(SUNLinearSolver S, SUNMatrix nul, + N_Vector x, N_Vector b, realtype tol) + int SUNLinSolNumIters_PCG(SUNLinearSolver S) + realtype SUNLinSolResNorm_PCG(SUNLinearSolver S) + N_Vector SUNLinSolResid_PCG(SUNLinearSolver S) + long int SUNLinSolLastFlag_PCG(SUNLinearSolver S) + int SUNLinSolSpace_PCG(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_PCG(SUNLinearSolver S) + + +cdef extern from "sunlinsol/sunlinsol_spbcgs.h": + + struct _SUNLinearSolverContent_SPBCGS: + int maxl + int pretype + int numiters + realtype resnorm + long int last_flag + + ATimesFn ATimes + void* ATData + PSetupFn Psetup + PSolveFn Psolve + void* PData + + N_Vector s1 + N_Vector s2 + N_Vector r + N_Vector r_star + N_Vector p + N_Vector q + N_Vector u + N_Vector Ap + N_Vector vtemp + + ctypedef _SUNLinearSolverContent_SPBCGS *SUNLinearSolverContent_SPBCGS + + SUNLinearSolver SUNSPBCGS(N_Vector y, int pretype, int maxl) + int SUNSPBCGSSetPrecType(SUNLinearSolver S, int pretype) + int SUNSPBCGSSetMaxl(SUNLinearSolver S, int maxl) + + SUNLinearSolver_Type SUNLinSolGetType_SPBCGS(SUNLinearSolver S) + int SUNLinSolInitialize_SPBCGS(SUNLinearSolver S) + int SUNLinSolSetATimes_SPBCGS(SUNLinearSolver S, void* A_data, + ATimesFn ATimes) + int SUNLinSolSetPreconditioner_SPBCGS(SUNLinearSolver S, + void* P_data, + PSetupFn Pset, + PSolveFn Psol) + int SUNLinSolSetScalingVectors_SPBCGS(SUNLinearSolver S, + N_Vector s1, + N_Vector s2) + int SUNLinSolSetup_SPBCGS(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_SPBCGS(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + int SUNLinSolNumIters_SPBCGS(SUNLinearSolver S) + realtype SUNLinSolResNorm_SPBCGS(SUNLinearSolver S) + N_Vector SUNLinSolResid_SPBCGS(SUNLinearSolver S) + long int SUNLinSolLastFlag_SPBCGS(SUNLinearSolver S) + int SUNLinSolSpace_SPBCGS(SUNLinearSolver S, + long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_SPBCGS(SUNLinearSolver S) + + +cdef extern from "sunlinsol/sunlinsol_spfgmr.h": + + struct _SUNLinearSolverContent_SPFGMR: + int maxl + int pretype + int gstype + int max_restarts + int numiters + realtype resnorm + long int last_flag + + ATimesFn ATimes + void* ATData + PSetupFn Psetup + PSolveFn Psolve + void* PData + + N_Vector s1 + N_Vector s2 + N_Vector *V + N_Vector *Z + realtype **Hes + realtype *givens + N_Vector xcor + realtype *yg + N_Vector vtemp + + ctypedef _SUNLinearSolverContent_SPFGMR *SUNLinearSolverContent_SPFGMR + + SUNLinearSolver SUNSPFGMR(N_Vector y, int pretype, int maxl) + int SUNSPFGMRSetPrecType(SUNLinearSolver S, int pretype) + int SUNSPFGMRSetGSType(SUNLinearSolver S, int gstype) + int SUNSPFGMRSetMaxRestarts(SUNLinearSolver S, int maxrs) + + SUNLinearSolver_Type SUNLinSolGetType_SPFGMR(SUNLinearSolver S) + int SUNLinSolInitialize_SPFGMR(SUNLinearSolver S) + int SUNLinSolSetATimes_SPFGMR(SUNLinearSolver S, void* A_data, + ATimesFn ATimes) + int SUNLinSolSetPreconditioner_SPFGMR(SUNLinearSolver S, + void* P_data, + PSetupFn Pset, + PSolveFn Psol) + int SUNLinSolSetScalingVectors_SPFGMR(SUNLinearSolver S, + N_Vector s1, + N_Vector s2) + int SUNLinSolSetup_SPFGMR(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_SPFGMR(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + int SUNLinSolNumIters_SPFGMR(SUNLinearSolver S) + realtype SUNLinSolResNorm_SPFGMR(SUNLinearSolver S) + N_Vector SUNLinSolResid_SPFGMR(SUNLinearSolver S) + long int SUNLinSolLastFlag_SPFGMR(SUNLinearSolver S) + int SUNLinSolSpace_SPFGMR(SUNLinearSolver S, + long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_SPFGMR(SUNLinearSolver S) + + +cdef extern from "sunlinsol/sunlinsol_spgmr.h": + + struct _SUNLinearSolverContent_SPGMR: + int maxl + int pretype + int gstype + int max_restarts + int numiters + realtype resnorm + long int last_flag + + ATimesFn ATimes + void* ATData + PSetupFn Psetup + PSolveFn Psolve + void* PData + + N_Vector s1 + N_Vector s2 + N_Vector *V + realtype **Hes + realtype *givens + N_Vector xcor + realtype *yg + N_Vector vtemp + + ctypedef _SUNLinearSolverContent_SPGMR *SUNLinearSolverContent_SPGMR + + SUNLinearSolver SUNSPGMR(N_Vector y, int pretype, int maxl) + int SUNSPGMRSetPrecType(SUNLinearSolver S, int pretype) + int SUNSPGMRSetGSType(SUNLinearSolver S, int gstype) + int SUNSPGMRSetMaxRestarts(SUNLinearSolver S, int maxrs) + + SUNLinearSolver_Type SUNLinSolGetType_SPGMR(SUNLinearSolver S) + int SUNLinSolInitialize_SPGMR(SUNLinearSolver S) + int SUNLinSolSetATimes_SPGMR(SUNLinearSolver S, void* A_data, + ATimesFn ATimes) + int SUNLinSolSetPreconditioner_SPGMR(SUNLinearSolver S, void* P_data, + PSetupFn Pset, PSolveFn Psol) + int SUNLinSolSetScalingVectors_SPGMR(SUNLinearSolver S, N_Vector s1, + N_Vector s2) + int SUNLinSolSetup_SPGMR(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_SPGMR(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + int SUNLinSolNumIters_SPGMR(SUNLinearSolver S) + realtype SUNLinSolResNorm_SPGMR(SUNLinearSolver S) + N_Vector SUNLinSolResid_SPGMR(SUNLinearSolver S) + long int SUNLinSolLastFlag_SPGMR(SUNLinearSolver S) + int SUNLinSolSpace_SPGMR(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_SPGMR(SUNLinearSolver S) + + +cdef extern from "sunlinsol/sunlinsol_sptfqmr.h": + + struct _SUNLinearSolverContent_SPTFQMR: + int maxl + int pretype + int numiters + realtype resnorm + long int last_flag + + ATimesFn ATimes + void* ATData + PSetupFn Psetup + PSolveFn Psolve + void* PData + + N_Vector s1 + N_Vector s2 + N_Vector r_star + N_Vector q + N_Vector d + N_Vector v + N_Vector p + N_Vector *r + N_Vector u + N_Vector vtemp1 + N_Vector vtemp2 + N_Vector vtemp3 + + ctypedef _SUNLinearSolverContent_SPTFQMR *SUNLinearSolverContent_SPTFQMR + + SUNLinearSolver SUNSPTFQMR(N_Vector y, int pretype, int maxl) + int SUNSPTFQMRSetPrecType(SUNLinearSolver S, int pretype) + int SUNSPTFQMRSetMaxl(SUNLinearSolver S, int maxl) + + SUNLinearSolver_Type SUNLinSolGetType_SPTFQMR(SUNLinearSolver S) + int SUNLinSolInitialize_SPTFQMR(SUNLinearSolver S) + int SUNLinSolSetATimes_SPTFQMR(SUNLinearSolver S, void* A_data, + ATimesFn ATimes) + int SUNLinSolSetPreconditioner_SPTFQMR(SUNLinearSolver S, void* P_data, + PSetupFn Pset, PSolveFn Psol) + int SUNLinSolSetScalingVectors_SPTFQMR(SUNLinearSolver S, N_Vector s1, + N_Vector s2) + int SUNLinSolSetup_SPTFQMR(SUNLinearSolver S, SUNMatrix A) + int SUNLinSolSolve_SPTFQMR(SUNLinearSolver S, SUNMatrix A, + N_Vector x, N_Vector b, realtype tol) + int SUNLinSolNumIters_SPTFQMR(SUNLinearSolver S) + realtype SUNLinSolResNorm_SPTFQMR(SUNLinearSolver S) + N_Vector SUNLinSolResid_SPTFQMR(SUNLinearSolver S) + long int SUNLinSolLastFlag_SPTFQMR(SUNLinearSolver S) + int SUNLinSolSpace_SPTFQMR(SUNLinearSolver S, long int *lenrwLS, + long int *leniwLS) + int SUNLinSolFree_SPTFQMR(SUNLinearSolver S) + + +# We don't use SuperLUMT - "slu_mt_ddefs.h" required +#cdef extern from "sunlinsol/sunlinsol_superlumt.h": +# +# struct _SUNLinearSolverContent_SuperLUMT: +# long int last_flag +# int first_factorize +# SuperMatrix *A, *AC, *L, *U, *B +# Gstat_t *Gstat +# sunindextype *perm_r, *perm_c +# sunindextype N +# int num_threads +# realtype diag_pivot_thresh +# int ordering +# superlumt_options_t *options +# +# ctypedef _SUNLinearSolverContent_SuperLUMT *SUNLinearSolverContent_SuperLUMT +# +# SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads) +# +# int SUNSuperLUMTSetOrdering(SUNLinearSolver S, int ordering_choice) +# SUNLinearSolver_Type SUNLinSolGetType_SuperLUMT(SUNLinearSolver S) +# int SUNLinSolInitialize_SuperLUMT(SUNLinearSolver S) +# int SUNLinSolSetup_SuperLUMT(SUNLinearSolver S, SUNMatrix A) +# int SUNLinSolSolve_SuperLUMT(SUNLinearSolver S, SUNMatrix A, +# N_Vector x, N_Vector b, realtype tol) +# long int SUNLinSolLastFlag_SuperLUMT(SUNLinearSolver S) +# int SUNLinSolSpace_SuperLUMT(SUNLinearSolver S, long int *lenrwLS, +# long int *leniwLS) +# int SUNLinSolFree_SuperLUMT(SUNLinearSolver S) diff --git a/scikits/odes/sundials/c_sunmatrix.pxd b/scikits/odes/sundials/c_sunmatrix.pxd new file mode 100644 index 0000000..2153bbd --- /dev/null +++ b/scikits/odes/sundials/c_sunmatrix.pxd @@ -0,0 +1,122 @@ +from libc.stdio cimport FILE +from c_sundials cimport * + +cdef extern from "sunmatrix/sunmatrix_dense.h": + + cdef struct _SUNMatrixContent_Dense: + sunindextype M + sunindextype N + realtype *data + sunindextype ldata + realtype **cols + + ctypedef _SUNMatrixContent_Dense *SUNMatrixContent_Dense + + SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N) + void SUNDenseMatrix_Print(SUNMatrix A, FILE* outfile) + sunindextype SUNDenseMatrix_Rows(SUNMatrix A) + sunindextype SUNDenseMatrix_Columns(SUNMatrix A) + sunindextype SUNDenseMatrix_LData(SUNMatrix A) + realtype* SUNDenseMatrix_Data(SUNMatrix A) + realtype** SUNDenseMatrix_Cols(SUNMatrix A) + realtype* SUNDenseMatrix_Column(SUNMatrix A, sunindextype j) + SUNMatrix_ID SUNMatGetID_Dense(SUNMatrix A) + SUNMatrix SUNMatClone_Dense(SUNMatrix A) + void SUNMatDestroy_Dense(SUNMatrix A) + int SUNMatZero_Dense(SUNMatrix A) + int SUNMatCopy_Dense(SUNMatrix A, SUNMatrix B) + int SUNMatScaleAdd_Dense(realtype c, SUNMatrix A, SUNMatrix B) + int SUNMatScaleAddI_Dense(realtype c, SUNMatrix A) + int SUNMatMatvec_Dense(SUNMatrix A, N_Vector x, N_Vector y) + int SUNMatSpace_Dense(SUNMatrix A, long int *lenrw, long int *leniw) + + +cdef extern from "sunmatrix/sunmatrix_band.h": + + cdef struct _SUNMatrixContent_Band: + sunindextype M + sunindextype N + sunindextype ldim + sunindextype mu + sunindextype ml + sunindextype s_mu + realtype *data + sunindextype ldata + realtype **cols + + ctypedef _SUNMatrixContent_Band *SUNMatrixContent_Band + + SUNMatrix SUNBandMatrix(sunindextype N, sunindextype mu, + sunindextype ml, sunindextype smu) + + void SUNBandMatrix_Print(SUNMatrix A, FILE* outfile) + + sunindextype SUNBandMatrix_Rows(SUNMatrix A) + sunindextype SUNBandMatrix_Columns(SUNMatrix A) + sunindextype SUNBandMatrix_LowerBandwidth(SUNMatrix A) + sunindextype SUNBandMatrix_UpperBandwidth(SUNMatrix A) + sunindextype SUNBandMatrix_StoredUpperBandwidth(SUNMatrix A) + sunindextype SUNBandMatrix_LDim(SUNMatrix A) + realtype* SUNBandMatrix_Data(SUNMatrix A) + realtype** SUNBandMatrix_Cols(SUNMatrix A) + realtype* SUNBandMatrix_Column(SUNMatrix A, sunindextype j) + + SUNMatrix_ID SUNMatGetID_Band(SUNMatrix A) + SUNMatrix SUNMatClone_Band(SUNMatrix A) + void SUNMatDestroy_Band(SUNMatrix A) + int SUNMatZero_Band(SUNMatrix A) + int SUNMatCopy_Band(SUNMatrix A, SUNMatrix B) + int SUNMatScaleAdd_Band(realtype c, SUNMatrix A, SUNMatrix B) + int SUNMatScaleAddI_Band(realtype c, SUNMatrix A) + int SUNMatMatvec_Band(SUNMatrix A, N_Vector x, N_Vector y) + int SUNMatSpace_Band(SUNMatrix A, long int *lenrw, long int *leniw) + + +cdef extern from "sunmatrix/sunmatrix_sparse.h": + + cdef struct _SUNMatrixContent_Sparse: + sunindextype M + sunindextype N + sunindextype NNZ + sunindextype NP + realtype *data + int sparsetype + sunindextype *indexvals + sunindextype *indexptrs + # CSC indices + sunindextype **rowvals + sunindextype **colptrs + # CSR indices + sunindextype **colvals + sunindextype **rowptrs + + ctypedef _SUNMatrixContent_Sparse *SUNMatrixContent_Sparse + + + SUNMatrix SUNSparseMatrix(sunindextype M, sunindextype N, + sunindextype NNZ, int sparsetype) + SUNMatrix SUNSparseFromDenseMatrix(SUNMatrix A, realtype droptol, + int sparsetype) + SUNMatrix SUNSparseFromBandMatrix(SUNMatrix A, realtype droptol, + int sparsetype) + int SUNSparseMatrix_Realloc(SUNMatrix A) + void SUNSparseMatrix_Print(SUNMatrix A, FILE* outfile) + + sunindextype SUNSparseMatrix_Rows(SUNMatrix A) + sunindextype SUNSparseMatrix_Columns(SUNMatrix A) + sunindextype SUNSparseMatrix_NNZ(SUNMatrix A) + sunindextype SUNSparseMatrix_NP(SUNMatrix A) + int SUNSparseMatrix_SparseType(SUNMatrix A) + realtype* SUNSparseMatrix_Data(SUNMatrix A) + sunindextype* SUNSparseMatrix_IndexValues(SUNMatrix A) + sunindextype* SUNSparseMatrix_IndexPointers(SUNMatrix A) + + SUNMatrix_ID SUNMatGetID_Sparse(SUNMatrix A) + SUNMatrix SUNMatClone_Sparse(SUNMatrix A) + void SUNMatDestroy_Sparse(SUNMatrix A) + int SUNMatZero_Sparse(SUNMatrix A) + int SUNMatCopy_Sparse(SUNMatrix A, SUNMatrix B) + int SUNMatScaleAdd_Sparse(realtype c, SUNMatrix A, SUNMatrix B) + int SUNMatScaleAddI_Sparse(realtype c, SUNMatrix A) + int SUNMatMatvec_Sparse(SUNMatrix A, N_Vector x, N_Vector y) + int SUNMatSpace_Sparse(SUNMatrix A, long int *lenrw, long int *leniw) diff --git a/scikits/odes/sundials/common_defs.pxd b/scikits/odes/sundials/common_defs.pxd index 8eaeb05..9ab0479 100644 --- a/scikits/odes/sundials/common_defs.pxd +++ b/scikits/odes/sundials/common_defs.pxd @@ -1,5 +1,5 @@ cimport numpy as np -from .c_sundials cimport N_Vector, DlsMat +from .c_sundials cimport N_Vector, DlsMat, SUNMatrix include "sundials_config.pxi" @@ -13,10 +13,17 @@ ELSE: # fall back to double ctypedef np.double_t DTYPE_t -cdef inline int nv_s2ndarray(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a) -cdef inline int ndarray2nv_s(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a) +IF SUNDIALS_INDEX_TYPE == "int32": + ctypedef np.int32_t INDEX_TYPE_t +ELIF SUNDIALS_INDEX_TYPE == "int64": + ctypedef np.int64_t INDEX_TYPE_t +ELSE: + ctypedef np.int64_t INDEX_TYPE_t + +cdef inline int nv_s2ndarray(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a) except? -1 +cdef inline int ndarray2nv_s(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a) except? -1 -cdef inline int DlsMatd2ndarray(DlsMat m, np.ndarray a) -cdef inline int ndarray2DlsMatd(DlsMat m, np.ndarray a) +cdef inline int SUNMatrix2ndarray(SUNMatrix m, np.ndarray a) except? -1 +cdef inline int ndarray2SUNMatrix(SUNMatrix m, np.ndarray a) except? -1 cdef ensure_numpy_float_array(object value) diff --git a/scikits/odes/sundials/common_defs.pyx b/scikits/odes/sundials/common_defs.pyx index d9aff6f..e63fa2d 100644 --- a/scikits/odes/sundials/common_defs.pyx +++ b/scikits/odes/sundials/common_defs.pyx @@ -4,10 +4,22 @@ import numpy as np cimport numpy as np import inspect from .c_sundials cimport ( - realtype, N_Vector, DlsMat, booleantype, SpfgmrMem, - NV_LENGTH_S as nv_length_s, NV_DATA_S as nv_data_s + realtype, sunindextype, N_Vector, DlsMat, booleantype, SpfgmrMem, + SUNMatrix, SUNMatGetID, SUNMATRIX_DENSE, SUNMATRIX_BAND, SUNMATRIX_SPARSE, + SUNMATRIX_CUSTOM, +) +from .c_nvector_serial cimport ( + N_VGetLength_Serial as nv_length_s, # use function not macro + NV_DATA_S as nv_data_s, +) +from .c_sunmatrix cimport ( + SUNDenseMatrix_Rows, SUNDenseMatrix_Columns, SUNDenseMatrix_Column, + SUNBandMatrix_Columns, SUNBandMatrix_UpperBandwidth, + SUNBandMatrix_LowerBandwidth, SUNBandMatrix_Column, ) +from libc.stdio cimport stderr + include "sundials_config.pxi" precision = SUNDIALS_FLOAT_TYPE @@ -21,6 +33,16 @@ ELSE: # fall back to double from numpy import double as DTYPE +index_precision = SUNDIALS_INDEX_TYPE +IF SUNDIALS_INDEX_TYPE == "int32": + from numpy import int32 as INDEX_TYPE +ELIF SUNDIALS_INDEX_TYPE == "int64": + from numpy import int64 as INDEX_TYPE +ELSE: + from numpy import int64 as INDEX_TYPE + +has_lapack = SUNDIALS_BLAS_LAPACK + ctypedef realtype *DlsMat_col ctypedef realtype *nv_content_data_s @@ -96,16 +118,16 @@ cdef inline N_Vector spfgmr_vtemp(SpfgmrMem mem): return mem.vtemp # Public functions -cdef inline int nv_s2ndarray(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a): - """ copy a serial N_Vector v to a nympy array a """ - cdef unsigned int N, i +cdef inline int nv_s2ndarray(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a) except? -1: + """ copy a serial N_Vector v to a numpy array a """ + cdef sunindextype N, i N = nv_length_s(v) cdef nv_content_data_s v_data = nv_data_s(v) for i in range(N): a[i] = get_nv_ith_s(v_data, i) -cdef inline int ndarray2nv_s(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a): +cdef inline int ndarray2nv_s(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a) except? -1: """ copy a numpy array a to a serial N_Vector v t""" cdef unsigned int N, i N = nv_length_s(v) @@ -114,28 +136,59 @@ cdef inline int ndarray2nv_s(N_Vector v, np.ndarray[DTYPE_t, ndim=1] a): for i in range(N): set_nv_ith_s(v_data, i, a[i]) -cdef inline int DlsMatd2ndarray(DlsMat m, np.ndarray a): - """ copy a Dense DlsMat m to a nympy array a """ - cdef unsigned int N, i, j +cdef inline int SUNMatrix2ndarray(SUNMatrix m, np.ndarray a) except? -1: + """ copy a SUNMatrix m to a numpy array a """ + cdef sunindextype N, M, i, j, ml, mu, stride cdef nv_content_data_s v_col - N = get_dense_N(m) + if SUNMatGetID(m) == SUNMATRIX_DENSE: + N = SUNDenseMatrix_Columns(m) + M = SUNDenseMatrix_Rows(m) - for i in range(N): - v_col = get_dense_col(m, i) for j in range(N): - a[i,j] = get_nv_ith_s(v_col, j) + v_col = SUNDenseMatrix_Column(m, j) + for i in range(M): + a[i,j] = get_nv_ith_s(v_col, i) + + elif SUNMatGetID(m) == SUNMATRIX_BAND: + N = SUNBandMatrix_Columns(m) + ml = SUNBandMatrix_LowerBandwidth(m) + mu = SUNBandMatrix_UpperBandwidth(m) + stride = ml + mu + 1 + for j in range(N): + v_col = SUNBandMatrix_Column(m, j) + for i in range(stride): + a[i,j] = v_col[i - mu] + + else: + raise NotImplementedError("SUNMatrix type not supported") -cdef inline int ndarray2DlsMatd(DlsMat m, np.ndarray a): - """ copy a nympy array a to a Dense DlsMat m""" - cdef unsigned int N, i, j +cdef inline int ndarray2SUNMatrix(SUNMatrix m, np.ndarray a) except? -1: + """ copy a numpy array a to a SUNMatrix m""" + cdef sunindextype N, M, i, j, ml, mu, stride cdef nv_content_data_s v_col - N = get_dense_N(m) + if SUNMatGetID(m) == SUNMATRIX_DENSE: + N = SUNDenseMatrix_Columns(m) + M = SUNDenseMatrix_Rows(m) - for i in range(N): for j in range(N): - set_dense_element(m, i, j, a[i,j]) + v_col = SUNDenseMatrix_Column(m, j) + for i in range(M): + v_col[i] = a[i,j] + + elif SUNMatGetID(m) == SUNMATRIX_BAND: + N = SUNBandMatrix_Columns(m) + ml = SUNBandMatrix_LowerBandwidth(m) + mu = SUNBandMatrix_UpperBandwidth(m) + stride = ml + mu + 1 + for j in range(N): + v_col = SUNBandMatrix_Column(m, j) + for i in range(stride): + v_col[i - mu] = a[i,j] + + else: + raise NotImplementedError("SUNMatrix type not supported") cdef ensure_numpy_float_array(object value): try: diff --git a/scikits/odes/sundials/cvode.pxd b/scikits/odes/sundials/cvode.pxd index 40092f9..0033dad 100644 --- a/scikits/odes/sundials/cvode.pxd +++ b/scikits/odes/sundials/cvode.pxd @@ -2,7 +2,7 @@ cimport numpy as np from .c_sundials cimport N_Vector, realtype -from .common_defs cimport DTYPE_t +from .common_defs cimport DTYPE_t, INDEX_TYPE_t cdef class CV_RhsFunction: cpdef int evaluate(self, DTYPE_t t, @@ -29,6 +29,7 @@ cdef class CV_WrapRootFunction(CV_RootFunction): cdef class CV_JacRhsFunction: cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] fy, np.ndarray[DTYPE_t, ndim=2] J) except? -1 cdef class CV_WrapJacRhsFunction(CV_JacRhsFunction): @@ -78,6 +79,18 @@ cdef class CV_WrapJacTimesVecFunction(CV_JacTimesVecFunction): cdef int with_userdata cpdef set_jac_times_vecfn(self, object jac_times_vecfn) +cdef class CV_JacTimesSetupFunction: + cpdef int evaluate(self, + DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] fy, + object userdata = *) except? -1 + +cdef class CV_WrapJacTimesSetupFunction(CV_JacTimesSetupFunction): + cpdef object _jac_times_setupfn + cdef int with_userdata + cpdef set_jac_times_setupfn(self, object jac_times_setupfn) + cdef class CV_ContinuationFunction: cpdef object _fn cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, @@ -105,6 +118,7 @@ cdef class CV_data: cdef CV_PrecSolveFunction prec_solvefn cdef CV_PrecSetupFunction prec_setupfn cdef CV_JacTimesVecFunction jac_times_vecfn + cdef CV_JacTimesSetupFunction jac_times_setupfn cdef bint parallel_implementation cdef object user_data cdef CV_ErrHandler err_handler @@ -117,7 +131,7 @@ cdef class CVODE: cdef bint parallel_implementation, initialized, _old_api, _step_compute, _validate_flags cdef CV_data aux_data - cdef long int N #problem size, i.e. len(y0) = N + cdef INDEX_TYPE_t N #problem size, i.e. len(y0) = N cdef N_Vector y0, y, yp # for 'step' method cdef list t_roots cdef list y_roots diff --git a/scikits/odes/sundials/cvode.pyx b/scikits/odes/sundials/cvode.pyx index 92318b6..3789564 100644 --- a/scikits/odes/sundials/cvode.pyx +++ b/scikits/odes/sundials/cvode.pyx @@ -5,6 +5,8 @@ from enum import IntEnum import inspect from warnings import warn +include "sundials_config.pxi" + import numpy as np cimport numpy as np @@ -14,11 +16,16 @@ from . import ( ) from .c_sundials cimport realtype, N_Vector +from .c_nvector_serial cimport * +from .c_sunmatrix cimport * +from .c_sunlinsol cimport * + from .c_cvode cimport * from .common_defs cimport ( - nv_s2ndarray, ndarray2nv_s, ndarray2DlsMatd, DTYPE_t, + nv_s2ndarray, ndarray2nv_s, ndarray2SUNMatrix, DTYPE_t, INDEX_TYPE_t, ) -from .common_defs import DTYPE # this is needed because we want DTYPE to be +from .common_defs import DTYPE, INDEX_TYPE +# this is needed because we want DTYPE and INDEX_TYPE to be # accessible from python (not only in cython) @@ -234,7 +241,8 @@ cdef class CV_JacRhsFunction: """ cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, - np.ndarray J) except? -1: + np.ndarray[DTYPE_t, ndim=1] fy, + np.ndarray[DTYPE_t, ndim=2] J) except? -1: """ Returns the Jacobi matrix of the right hand side function, as d(rhs)/d y @@ -256,6 +264,7 @@ cdef class CV_WrapJacRhsFunction(CV_JacRhsFunction): cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] fy, np.ndarray J) except? -1: """ Returns the Jacobi matrix (for dense the full matrix, for band only @@ -266,18 +275,20 @@ cdef class CV_WrapJacRhsFunction(CV_JacRhsFunction): ## self._jacfn(t, y, ydot, cj, J, userdata) ## else: ## self._jacfn(t, y, ydot, cj, J) - user_flag = self._jacfn(t, y, J) + user_flag = self._jacfn(t, y, fy, J) if user_flag is None: user_flag = 0 return user_flag -cdef int _jacdense(long int Neq, realtype tt, - N_Vector yy, N_Vector ff, DlsMat Jac, +cdef int _jacdense(realtype tt, + N_Vector yy, N_Vector ff, SUNMatrix Jac, void *auxiliary_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except? -1: - """function with the signature of CVDlsDenseJacFn that calls python Jac""" - cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp - cdef np.ndarray jac_tmp + """function with the signature of CVDlsJacFn that calls python Jac + Note: signature of Jac is SUNMatrix + """ + cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp, ff_tmp + cdef np.ndarray[DTYPE_t, ndim=2] jac_tmp aux_data = auxiliary_data cdef bint parallel_implementation = aux_data.parallel_implementation @@ -285,19 +296,18 @@ cdef int _jacdense(long int Neq, realtype tt, raise NotImplemented else: yy_tmp = aux_data.yy_tmp - if aux_data.jac_tmp is None: - N = np.alen(yy_tmp) - aux_data.jac_tmp = np.empty((N,N), DTYPE) jac_tmp = aux_data.jac_tmp nv_s2ndarray(yy, yy_tmp) - user_flag = aux_data.jac.evaluate(tt, yy_tmp, jac_tmp) + ff_tmp = aux_data.z_tmp + nv_s2ndarray(ff, ff_tmp) + + user_flag = aux_data.jac.evaluate(tt, yy_tmp, ff_tmp, jac_tmp) if parallel_implementation: raise NotImplemented else: - #we convert the python jac_tmp array to DslMat of sundials - ndarray2DlsMatd(Jac, jac_tmp) + ndarray2SUNMatrix(Jac, jac_tmp) return user_flag @@ -359,8 +369,9 @@ class MutableBool(object): def __init__(self, value): self.value = value -cdef int _prec_setupfn(realtype tt, N_Vector yy, N_Vector ff, booleantype jok, booleantype *jcurPtr, - realtype gamma, void *auxiliary_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except? -1: +cdef int _prec_setupfn(realtype tt, N_Vector yy, N_Vector ff, booleantype jok, + booleantype *jcurPtr, realtype gamma, + void *auxiliary_data) except? -1: """ function with the signature of CVSpilsPrecSetupFn, that calls python function """ cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp @@ -440,7 +451,7 @@ cdef class CV_WrapPrecSolveFunction(CV_PrecSolveFunction): return user_flag cdef int _prec_solvefn(realtype tt, N_Vector yy, N_Vector ff, N_Vector r, N_Vector z, - realtype gamma, realtype delta, int lr, void *auxiliary_data, N_Vector tmp) except? -1: + realtype gamma, realtype delta, int lr, void *auxiliary_data) except? -1: """ function with the signature of CVSpilsPrecSolveFn, that calls python function """ cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp, r_tmp, z_tmp @@ -452,14 +463,6 @@ cdef int _prec_solvefn(realtype tt, N_Vector yy, N_Vector ff, N_Vector r, N_Vect else: yy_tmp = aux_data.yy_tmp - if aux_data.r_tmp is None: - N = np.alen(yy_tmp) - aux_data.r_tmp = np.empty(N, DTYPE) - - if aux_data.z_tmp is None: - N = np.alen(yy_tmp) - aux_data.z_tmp = np.empty(N, DTYPE) - r_tmp = aux_data.r_tmp z_tmp = aux_data.z_tmp @@ -543,14 +546,6 @@ cdef int _jac_times_vecfn(N_Vector v, N_Vector Jv, realtype t, N_Vector y, else: y_tmp = aux_data.yy_tmp - if aux_data.r_tmp is None: - N = np.alen(y_tmp) - aux_data.r_tmp = np.empty(N, DTYPE) - - if aux_data.z_tmp is None: - N = np.alen(y_tmp) - aux_data.z_tmp = np.empty(N, DTYPE) - v_tmp = aux_data.r_tmp Jv_tmp = aux_data.z_tmp @@ -566,6 +561,86 @@ cdef int _jac_times_vecfn(N_Vector v, N_Vector Jv, realtype t, N_Vector y, return user_flag +# JacTimesVec function +cdef class CV_JacTimesSetupFunction: + """ + Prototype for jacobian times setup function. + + Note that evaluate must return a integer, 0 for success, non-zero for error + (as per CVODE documentation), with >0 a recoverable error (step is retried). + """ + cpdef int evaluate(self, + DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] fy, + object userdata = None) except? -1: + """ + This function calculates the product of the Jacobian with a given vector v. + Use the userdata object to expose Jacobian related data to the solve function. + + This is a generic class, you should subclass it for the problem specific + purposes. + """ + return 0 + +cdef class CV_WrapJacTimesSetupFunction(CV_JacTimesSetupFunction): + cpdef set_jac_times_setupfn(self, object jac_times_setupfn): + """ + Set some CV_JacTimesSetupFn executable class. + """ + """ + set a jacobian-times-vector method setup as a CV_JacTimesSetupFunction + executable class + """ + self.with_userdata = 0 + nrarg = _get_num_args(jac_times_setupfn) + if nrarg > 4: + #hopefully a class method, self gives 5 arg! + self.with_userdata = 1 + elif nrarg == 4 and inspect.isfunction(jac_times_setupfn): + self.with_userdata = 1 + self._jac_times_setupfn = jac_times_setupfn + + cpdef int evaluate(self, + DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] fy, + object userdata = None) except? -1: + if self.with_userdata == 1: + user_flag = self._jac_times_setupfn(t, y, fy, userdata) + else: + user_flag = self._jac_times_setupfn(t, y, fy) + if user_flag is None: + user_flag = 0 + return user_flag + +cdef int _jac_times_setupfn(realtype t, N_Vector y, N_Vector fy, + void *user_data) except? -1: + """ function with the signature of CVSpilsJacTimesSetupFn, that calls python function """ + cdef np.ndarray[DTYPE_t, ndim=1] y_tmp, fy_tmp + + aux_data = user_data + cdef bint parallel_implementation = aux_data.parallel_implementation + + if parallel_implementation: + raise NotImplemented + else: + y_tmp = aux_data.yy_tmp + + fy_tmp = aux_data.z_tmp + + nv_s2ndarray(y, y_tmp) + nv_s2ndarray(fy, fy_tmp) + + user_flag = aux_data.jac_times_setupfn.evaluate(t, y_tmp, fy_tmp, aux_data.user_data) + + #if parallel_implementation: + # raise NotImplemented + #else: + # ndarray2nv_s(fy, fy_tmp) + + return user_flag + cdef class CV_ContinuationFunction: """ Simple wrapper for functions called when ROOT or TSTOP are returned. @@ -639,8 +714,8 @@ cdef class CV_data: self.yp_tmp = np.empty(N, DTYPE) self.jac_tmp = None self.g_tmp = None - self.r_tmp = None - self.z_tmp = None + self.r_tmp = np.empty(N, DTYPE) + self.z_tmp = np.empty(N, DTYPE) cdef class CVODE: @@ -683,6 +758,7 @@ cdef class CVODE: 'prec_setupfn': None, 'prec_solvefn': None, 'jac_times_vecfn': None, + 'jac_times_setupfn': None, 'err_handler': None, 'err_user_data': None, 'old_api': None, @@ -768,9 +844,11 @@ cdef class CVODE: Defines the jacobian function and has to be a subclass of CV_JacRhsFunction class or python function. This function takes as input arguments current time t, current value of y, + current value of f(t,y), and a 2D numpy array of returned jacobian and optional userdata. Return value is 0 if successfull. Jacobian is only used for dense or lapackdense linear solver + TODO: cvode supports Jacobian for band also, this is not supported. 'rtol': Values: float, 1e-6 = default Description: @@ -810,13 +888,14 @@ cdef class CVODE: of 0.0 uses the solver's internal default value. 'linsolver': Values: 'dense' (= default), 'lapackdense', 'band', - 'lapackband', 'spgmr', 'spbcg', 'sptfqmr' + 'lapackband', 'spgmr', 'spbcgs', 'sptfqmr' Description: Specifies used linear solver. Limitations: Linear solvers for dense and band matrices can be used only for serial implementation. For parallel implementation use_relaxation use lapackdense or lapackband respectively. + TODO: to add new solvers: pcg, spfgmr, superlumt, klu 'lband', 'uband': Values: non-negative integer, 0 = default Description: @@ -831,7 +910,7 @@ cdef class CVODE: Values: 0 (= default), 1, 2, 3, 4, 5 Description: Dimension of the number of used Krylov subspaces - (used only by 'spgmr', 'spbcg', 'sptfqmr' linsolvers) + (used only by 'spgmr', 'spbcgs', 'sptfqmr' linsolvers) 'tstop': Values: float, 0.0 = default Description: @@ -883,6 +962,11 @@ cdef class CVODE: This function takes as input arguments the vector v, result vector Jv, current time t, current value of y, and optional userdata. + 'jac_times_setupfn': + Values: function of class CV_JacTimesSetupFunction + Description: + Optional. Default is to internal finite difference with no + extra setup. 'bdf_stability_detection': default = False, only used if lmm_type == 'bdf 'max_conv_fails': @@ -1177,8 +1261,8 @@ cdef class CVODE: (opts['implementation'].lower() == 'parallel') if self.parallel_implementation: raise ValueError('Error: Parallel implementation not implemented !') - cdef long int N - N = np.alen(y0) + cdef INDEX_TYPE_t N + N = np.alen(y0) if opts['rfn'] is None: raise ValueError('The right-hand-side function rfn not assigned ' @@ -1304,6 +1388,14 @@ cdef class CVODE: opts['jac_times_vecfn'] = tmpfun self.aux_data.jac_times_vecfn = jac_times_vecfn + jac_times_setupfn = opts['jac_times_setupfn'] + if jac_times_setupfn is not None and not isinstance(jac_times_setupfn, CV_JacTimesSetupFunction): + tmpfun = CV_WrapJacTimesSetupFunction() + tmpfun.set_jac_times_setupfn(jac_times_setupfn) + jac_times_setupfn = tmpfun + opts['jac_times_setupfn'] = tmpfun + self.aux_data.jac_times_setupfn = jac_times_setupfn + self.aux_data.user_data = opts['user_data'] # As cv_mem is now initialized, set also options changeable at runtime @@ -1338,48 +1430,37 @@ cdef class CVODE: if iter_type == 'newton': if linsolver == 'dense': - if self.parallel_implementation: - raise ValueError('Linear solver for dense matrices can be' - 'used only for serial implementation. For ' - 'parallel implementation use ''lapackdense''.') - else: - flag = CVDense(cv_mem, N) - if flag == CVDLS_ILL_INPUT: - raise ValueError('CVDense solver is not compatible with' - ' the current nvector implementation.') - elif flag == CVDLS_MEM_FAIL: - raise MemoryError('CVDense memory allocation error.') - elif linsolver == 'lapackdense': - flag = CVLapackDense(cv_mem, N) + A = SUNDenseMatrix(N, N) + LS = SUNDenseLinearSolver(self.y0, A) + # check if memory was allocated + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + # attach matrix and linear solver to cvode + flag = CVDlsSetLinearSolver(cv_mem, LS, A) if flag == CVDLS_ILL_INPUT: - raise ValueError('CVLapackDense solver is not compatible ' - 'with the current nvector implementation.') + raise ValueError('CVDense linear solver setting failed, ' + 'arguments incompatible') elif flag == CVDLS_MEM_FAIL: - raise MemoryError('CVLapackDense memory allocation error.') + raise MemoryError('CVDense linear solver memory allocation error.') + elif flag != CVDLS_SUCCESS: + raise ValueError('CVDlsSetLinearSolver failed with code {}' + .format(flag)) elif linsolver == 'band': - if self.parallel_implementation: - raise ValueError('Linear solver for band matrices can be ' - 'used only for serial implementation. ' - 'Use ''lapackband'' instead for parallel ' - 'implementation.') - else: - flag = CVBand(cv_mem, N, opts['uband'], - opts['lband']) - if flag == CVDLS_ILL_INPUT: - raise ValueError('CVBand solver is not compatible with ' - 'the current nvector implementation ' - 'or bandwith outside range.') - elif flag == CVDLS_MEM_FAIL: - raise MemoryError('CVBand memory allocation error.') - elif linsolver == 'lapackband': - flag = CVLapackBand(cv_mem, N, opts['uband'], - opts['lband']) + A = SUNBandMatrix(N, opts['uband'], opts['lband'], + opts['uband'] + opts['lband']); + LS = SUNBandLinearSolver(self.y0, A); + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + flag = CVDlsSetLinearSolver(cv_mem, LS, A) + if flag == CVDLS_ILL_INPUT: - raise ValueError('CVLapackBand solver is not compatible' - ' with the current nvector implementation' - ' or bandwith outside range.') + raise ValueError('CVBand linear solver setting failed, ' + 'arguments incompatible') elif flag == CVDLS_MEM_FAIL: - raise MemoryError('CVLapackBand memory allocation error.') + raise MemoryError('CVBand linear solver memory allocation error.') + elif flag != CVDLS_SUCCESS: + raise ValueError('CVDlsSetLinearSolver failed with code {}' + .format(flag)) elif linsolver == 'diag': flag = CVDiag(cv_mem) if flag == CVDIAG_ILL_INPUT: @@ -1387,7 +1468,10 @@ cdef class CVODE: ' the current nvector implementation.') elif flag == CVDIAG_MEM_FAIL: raise MemoryError('CVDiag memory allocation error.') - elif ((linsolver == 'spgmr') or (linsolver == 'spbcg') + elif flag != CVDIAG_SUCCESS: + raise ValueError('CVDiag failed with code {}' + .format(flag)) + elif ((linsolver == 'spgmr') or (linsolver == 'spbcgs') or (linsolver == 'sptfqmr')): precond_type = opts['precond_type'].lower() if precond_type == 'none': @@ -1403,16 +1487,32 @@ cdef class CVODE: % opts['precond_type']) if linsolver == 'spgmr': - flag = CVSpgmr(cv_mem, pretype, opts['maxl']) - elif linsolver == 'spbcg': - flag = CVSpbcg(cv_mem, pretype, opts['maxl']) + LS = SUNSPGMR(self.y0, pretype, opts['maxl']); + if LS == NULL: + raise ValueError('Could not allocate linear solver') + elif linsolver == 'spbcgs': + LS = SUNSPBCGS(self.y0, pretype, opts['maxl']); + if LS == NULL: + raise ValueError('Could not allocate linear solver') + elif linsolver == 'sptfqmr': + LS = SUNSPTFQMR(self.y0, pretype, opts['maxl']); + if LS == NULL: + raise ValueError('Could not allocate linear solver') else: - flag = CVSptfqmr(cv_mem, pretype, opts['maxl']) - + raise ValueError('Given linsolver {} not implemented in odes'.format(linsolver)) + + flag = CVSpilsSetLinearSolver(cv_mem, LS); if flag == CVSPILS_MEM_FAIL: raise MemoryError('LinSolver:CVSpils memory allocation ' 'error.') - + elif flag != CVSPILS_SUCCESS: + raise ValueError('CVSpilsSetLinearSolver failed with code {}' + .format(flag)) + # TODO: make option for the Gram-Schmidt orthogonalization + #flag = SUNSPGMRSetGSType(LS, gstype); + + # TODO make option + #flag = CVSpilsSetEpsLin(cvode_mem, DELT); if self.aux_data.prec_solvefn: if self.aux_data.prec_setupfn: flag = CVSpilsSetPreconditioner(cv_mem, _prec_setupfn, _prec_solvefn) @@ -1425,19 +1525,71 @@ cdef class CVODE: 'not been initialized.') if self.aux_data.jac_times_vecfn: - flag = CVSpilsSetJacTimesVecFn(cv_mem, _jac_times_vecfn) + if self.aux_data.jac_times_setupfn: + flag = CVSpilsSetJacTimes(cv_mem, _jac_times_setupfn, _jac_times_vecfn) + else: + flag = CVSpilsSetJacTimes(cv_mem, NULL, _jac_times_vecfn) if flag == CVSPILS_MEM_NULL: raise ValueError('LinSolver: The cvode mem pointer is NULL.') elif flag == CVSPILS_LMEM_NULL: raise ValueError('LinSolver: The cvspils linear solver has ' 'not been initialized.') - + elif flag != CVSPILS_SUCCESS: + raise ValueError('CVSpilsSetJacTimes failed with code {}' + .format(flag)) else: - raise ValueError('LinSolver: Unknown solver type: %s' - % opts['linsolver']) - - if (linsolver in ['dense', 'lapackdense']) and self.aux_data.jac: - CVDlsSetDenseJacFn(cv_mem, _jacdense) + IF SUNDIALS_BLAS_LAPACK: + if linsolver == 'lapackdense': + A = SUNDenseMatrix(N, N) + LS = SUNLapackDense(self.y0, A) + # check if memory was allocated + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + # attach matrix and linear solver to cvode + flag = CVDlsSetLinearSolver(cv_mem, LS, A) + if flag == CVDLS_ILL_INPUT: + raise ValueError('CVDense lapack linear solver setting failed, ' + 'arguments incompatible') + elif flag == CVDLS_MEM_FAIL: + raise MemoryError('CVDense lapack linear solver memory allocation error.') + elif flag != CVDLS_SUCCESS: + raise ValueError('CVDlsSetLinearSolver failed with code {}' + .format(flag)) + elif linsolver == 'lapackband': + A = SUNBandMatrix(N, opts['uband'], opts['lband'], + opts['uband'] + opts['lband']) + LS = SUNLapackBand(self.y0, A) + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + flag = CVDlsSetLinearSolver(cv_mem, LS, A) + if flag == CVDLS_ILL_INPUT: + raise ValueError('CVLapackBand linear solver setting failed, ' + 'arguments incompatible') + elif flag == CVDLS_MEM_FAIL: + raise MemoryError('CVLapackBand linear solver memory allocation error.') + elif flag != CVDLS_SUCCESS: + raise ValueError('CVDlsSetLinearSolver failed with code {}' + .format(flag)) + else: + raise ValueError('LinSolver: Unknown solver type: %s' + % opts['linsolver']) + ELSE: + raise ValueError('LinSolver: Unknown solver type: %s' + % opts['linsolver']) + + if (linsolver in ['dense', 'lapackdense', 'lapackband', 'band'] + and self.aux_data.jac): + # we need to create the correct shape for jacobian output, here is + # the best place + if linsolver == 'lapackband' or linsolver == 'band': + self.aux_data.jac_tmp = np.empty(( + opts['uband'] + opts['lband'] + 1, + np.alen(y0) + ), DTYPE + ) + else: + self.aux_data.jac_tmp = np.empty((np.alen(y0), np.alen(y0)), DTYPE) + CVDlsSetJacFn(cv_mem, _jacdense) #we test if jac don't give errors due to bad coding, as #cvode will ignore errors, it only checks return value (0 or 1 for error) @@ -1447,7 +1599,8 @@ cdef class CVODE: DTYPE) else: _test = np.empty((np.alen(y0), np.alen(y0)), DTYPE) - jac._jacfn(t0, y0, _test) + _fy_test = np.zeros(np.alen(y0), DTYPE) + jac._jacfn(t0, y0, _fy_test, _test) _test = None #now we initialize storage which is persistent over steps @@ -1513,8 +1666,8 @@ cdef class CVODE: self._init_step(t0, y0) return - cdef long int N - N = np.alen(y0) + cdef INDEX_TYPE_t N + N = np.alen(y0) if N == self.N: self.y0 = N_VMake_Serial(N, y0.data) else: @@ -1822,7 +1975,7 @@ cdef class CVODE: cdef int qlast, qcur cdef realtype hinused, hlast, hcur, tcur - # for extra output from SPILS modules (SPGMR, SPBCG, SPTFQMR) + # for extra output from SPILS modules (SPGMR, SPBCGS, SPTFQMR) cdef long int npevals, npsolves, njvevals, nliters, nfevalsLS flagCV = CVodeGetIntegratorStats(self._cv_mem, &nsteps, &nfevals, @@ -1838,7 +1991,7 @@ cdef class CVODE: 'LastStep': hlast, 'CurrentStep': hcur, 'CurrentStep': tcur} linsolver = self.options['linsolver'].lower() - if linsolver == 'spgmr' or linsolver == 'spbcg' or linsolver == 'sptfqmr': + if linsolver == 'spgmr' or linsolver == 'spbcgs' or linsolver == 'sptfqmr': flagCV = CVSpilsGetNumPrecEvals(self._cv_mem, &npevals) flagCV = CVSpilsGetNumPrecSolves(self._cv_mem, &npsolves) flagCV = CVSpilsGetNumJtimesEvals(self._cv_mem, &njvevals) diff --git a/scikits/odes/sundials/ida.pxd b/scikits/odes/sundials/ida.pxd index 6f25c61..998bca0 100644 --- a/scikits/odes/sundials/ida.pxd +++ b/scikits/odes/sundials/ida.pxd @@ -30,13 +30,43 @@ cdef class IDA_JacRhsFunction: cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, np.ndarray[DTYPE_t, ndim=1] ydot, + np.ndarray[DTYPE_t, ndim=1] residual, DTYPE_t cj, np.ndarray[DTYPE_t, ndim=2] J) except? -1 cdef class IDA_WrapJacRhsFunction(IDA_JacRhsFunction): cpdef object _jacfn cdef int with_userdata - cpdef set_jacfn(self, object jacfn) + cpdef set_jacfn(self, object jacfn) + +cdef class IDA_PrecSetupFunction: + cpdef int evaluate(self, DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] yp, + np.ndarray[DTYPE_t, ndim=1] rr, + DTYPE_t cj, + object userdata = *) except? -1 + +cdef class IDA_WrapPrecSetupFunction(IDA_PrecSetupFunction): + cpdef object _prec_setupfn + cdef int with_userdata + cpdef set_prec_setupfn(self, object prec_setupfn) + +cdef class IDA_PrecSolveFunction: + cpdef int evaluate(self, DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] yp, + np.ndarray[DTYPE_t, ndim=1] r, + np.ndarray[DTYPE_t, ndim=1] rvec, + np.ndarray[DTYPE_t, ndim=1] z, + DTYPE_t cj, + DTYPE_t delta, + object userdata = *) except? -1 + +cdef class IDA_WrapPrecSolveFunction(IDA_PrecSolveFunction): + cpdef object _prec_solvefn + cdef int with_userdata + cpdef set_prec_solvefn(self, object prec_solvefn) cdef class IDA_ContinuationFunction: cpdef object _fn @@ -61,10 +91,12 @@ cdef class IDA_WrapErrHandler(IDA_ErrHandler): cdef class IDA_data: - cdef np.ndarray yy_tmp, yp_tmp, residual_tmp, jac_tmp, g_tmp + cdef np.ndarray yy_tmp, yp_tmp, residual_tmp, jac_tmp, g_tmp, z_tmp, rvec_tmp cdef IDA_RhsFunction res cdef IDA_JacRhsFunction jac cdef IDA_RootFunction rootfn + cdef IDA_PrecSetupFunction prec_setupfn + cdef IDA_PrecSolveFunction prec_solvefn cdef bint parallel_implementation cdef object user_data cdef IDA_ErrHandler err_handler diff --git a/scikits/odes/sundials/ida.pyx b/scikits/odes/sundials/ida.pyx index c41e9e7..314f1c7 100644 --- a/scikits/odes/sundials/ida.pyx +++ b/scikits/odes/sundials/ida.pyx @@ -5,20 +5,29 @@ from enum import IntEnum import inspect from warnings import warn +include "sundials_config.pxi" + import numpy as np cimport numpy as np + from .c_sundials cimport realtype, N_Vector +from .c_nvector_serial cimport * +from .c_sunmatrix cimport * +from .c_sunlinsol cimport * + from .c_ida cimport * from .common_defs cimport ( - nv_s2ndarray, ndarray2nv_s, ndarray2DlsMatd, DTYPE_t, + nv_s2ndarray, ndarray2nv_s, ndarray2SUNMatrix, DTYPE_t, INDEX_TYPE_t, ) -from .common_defs import DTYPE # this is needed because we want DTYPE to be +from .common_defs import DTYPE, INDEX_TYPE +# this is needed because we want DTYPE and INDEX_TYPE to be # accessible from python (not only in cython) from . import ( IDASolveFailed, IDASolveFoundRoot, IDASolveReachedTSTOP, _get_num_args, ) + # TODO: parallel implementation: N_VectorParallel # TODO: linsolvers: check the output value for errors # TODO: optimize code for compiler @@ -65,6 +74,8 @@ class StatusEnumIDA(IntEnum): BAD_T = IDA_BAD_T BAD_DKY = IDA_BAD_DKY + UNRECOGNISED_ERROR = IDA_UNRECOGNISED_ERROR + STATUS_MESSAGE = { StatusEnumIDA.SUCCESS: "Successful function return.", StatusEnumIDA.TSTOP_RETURN: "Reached specified stopping point", @@ -92,6 +103,7 @@ STATUS_MESSAGE = { StatusEnumIDA.BAD_K: "Illegal value for k. If the requested k is not in the range 0,1,...,order used ", StatusEnumIDA.BAD_T: "Illegal value for t. If t is not within the interval of the last step taken.", StatusEnumIDA.BAD_DKY: "The dky vector is NULL", + StatusEnumIDA.UNRECOGNISED_ERROR: "Unrecognised error", } @@ -250,6 +262,7 @@ cdef class IDA_JacRhsFunction: cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, np.ndarray[DTYPE_t, ndim=1] ydot, + np.ndarray[DTYPE_t, ndim=1] residual, DTYPE_t cj, np.ndarray J) except? -1: """ @@ -274,6 +287,7 @@ cdef class IDA_WrapJacRhsFunction(IDA_JacRhsFunction): cpdef int evaluate(self, DTYPE_t t, np.ndarray[DTYPE_t, ndim=1] y, np.ndarray[DTYPE_t, ndim=1] ydot, + np.ndarray[DTYPE_t, ndim=1] residual, DTYPE_t cj, np.ndarray J) except? -1: """ @@ -285,16 +299,16 @@ cdef class IDA_WrapJacRhsFunction(IDA_JacRhsFunction): ## self._jacfn(t, y, ydot, cj, J, userdata) ## else: ## self._jacfn(t, y, ydot, cj, J) - user_flag = self._jacfn(t, y, ydot, cj, J) + user_flag = self._jacfn(t, y, ydot, residual, cj, J) if user_flag is None: user_flag = 0 return user_flag -cdef int _jacdense(long int Neq, realtype tt, realtype cj, - N_Vector yy, N_Vector yp, N_Vector rr, DlsMat Jac, +cdef int _jacdense(realtype tt, realtype cj, + N_Vector yy, N_Vector yp, N_Vector rr, SUNMatrix Jac, void *auxiliary_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3): - """function with the signature of IDADlsDenseJacFn """ - cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp, yp_tmp + """function with the signature of IDADlsJacFn """ + cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp, yp_tmp, residual_tmp cdef np.ndarray jac_tmp aux_data = auxiliary_data @@ -304,6 +318,7 @@ cdef int _jacdense(long int Neq, realtype tt, realtype cj, else: yy_tmp = aux_data.yy_tmp yp_tmp = aux_data.yp_tmp + residual_tmp = aux_data.residual_tmp if aux_data.jac_tmp is None: N = np.alen(yy_tmp) aux_data.jac_tmp = np.empty((N,N), DTYPE) @@ -311,13 +326,199 @@ cdef int _jacdense(long int Neq, realtype tt, realtype cj, nv_s2ndarray(yy, yy_tmp) nv_s2ndarray(yp, yp_tmp) - user_flag = aux_data.jac.evaluate(tt, yy_tmp, yp_tmp, cj, jac_tmp) + nv_s2ndarray(rr, residual_tmp) + user_flag = aux_data.jac.evaluate(tt, yy_tmp, yp_tmp, residual_tmp, cj, jac_tmp) + + if parallel_implementation: + raise NotImplemented + else: + ndarray2SUNMatrix(Jac, jac_tmp) + + return user_flag + + +# Precondioner setup funtion +cdef class IDA_PrecSetupFunction: + """ + Prototype for preconditioning setup function. + + Note that evaluate must return a integer, 0 for success, positive for + recoverable failure, negative for unrecoverable failure (as per CVODE + documentation). + """ + cpdef int evaluate(self, DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] yp, + np.ndarray[DTYPE_t, ndim=1] rr, + DTYPE_t cj, + object userdata = None) except? -1: + """ + This function preprocesses and/or evaluates Jacobian-related data + needed by the preconditioner. Use the userdata object to expose + the preprocessed data to the solve function. + + This is a generic class, you should subclass it for the problem specific + purposes. + """ + return 0 + +cdef class IDA_WrapPrecSetupFunction(IDA_PrecSetupFunction): + cpdef set_prec_setupfn(self, object prec_setupfn): + """ + set a precondititioning setup method as a IDA_PrecSetupFunction + executable class + """ + self.with_userdata = 0 + nrarg = _get_num_args(prec_setupfn) + if nrarg > 5: + #hopefully a class method, self gives 6 arg! + self.with_userdata = 1 + elif nrarg == 5 and inspect.isfunction(prec_setupfn): + self.with_userdata = 1 + self._prec_setupfn = prec_setupfn + + cpdef int evaluate(self, DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] yp, + np.ndarray[DTYPE_t, ndim=1] rr, + DTYPE_t cj, + object userdata = None) except? -1: + if self.with_userdata == 1: + user_flag = self._prec_setupfn(t, y, yp, rr, cj, userdata) + else: + user_flag = self._prec_setupfn(t, y, yp, rr, cj) + if user_flag is None: + user_flag = 0 + return user_flag + +cdef int _prec_setupfn(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, + realtype cj, + void *auxiliary_data): + """ function with the signature of IDASpilsPrecSetupFn, that calls + the python function """ + cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp, rp_tmp, residual_tmp + + aux_data = auxiliary_data + cdef bint parallel_implementation = aux_data.parallel_implementation + + if parallel_implementation: + raise NotImplemented + else: + yy_tmp = aux_data.yy_tmp + yp_tmp = aux_data.yp_tmp + residual_tmp = aux_data.residual_tmp + nv_s2ndarray(yy, yy_tmp) + nv_s2ndarray(yp, yp_tmp) + nv_s2ndarray(rr, residual_tmp) + + user_flag = aux_data.prec_setupfn.evaluate(tt, yy_tmp, yp_tmp, + residual_tmp, cj, aux_data.user_data) + return user_flag + +# Precondioner solve funtion +cdef class IDA_PrecSolveFunction: + """ + Prototype for precondititioning solution function. + + Note that evaluate must return a integer, 0 for success, positive for + recoverable failure, negative for unrecoverable failure (as per CVODE + documentation). + """ + cpdef int evaluate(self, DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] yp, + np.ndarray[DTYPE_t, ndim=1] r, + np.ndarray[DTYPE_t, ndim=1] rvec, + np.ndarray[DTYPE_t, ndim=1] z, + DTYPE_t cj, + DTYPE_t delta, + object userdata = None) except? -1: + """ + This function solves the preconditioned system P*z = r, where P may be + either a left or right preconditioner matrix. Here P should approximate + (at least crudely) the Newton matrix M = I − gamma*J, where J is the + Jacobian of the system. If preconditioning is done on both sides, + the product of the two preconditioner matrices should approximate M. + + This is a generic class, you should subclass it for the problem specific + purposes. + """ + return 0 + +cdef class IDA_WrapPrecSolveFunction(IDA_PrecSolveFunction): + cpdef set_prec_solvefn(self, object prec_solvefn): + """ + set a precondititioning solve method as a IDA_PrecSolveFunction + executable class + """ + self.with_userdata = 0 + nrarg = _get_num_args(prec_solvefn) + if nrarg > 9: + #hopefully a class method, self gives 10 arg! + self.with_userdata = 1 + elif nrarg == 9 and inspect.isfunction(prec_solvefn): + self.with_userdata = 1 + self._prec_solvefn = prec_solvefn + + cpdef int evaluate(self, DTYPE_t t, + np.ndarray[DTYPE_t, ndim=1] y, + np.ndarray[DTYPE_t, ndim=1] yp, + np.ndarray[DTYPE_t, ndim=1] r, + np.ndarray[DTYPE_t, ndim=1] rvec, + np.ndarray[DTYPE_t, ndim=1] z, + DTYPE_t cj, + DTYPE_t delta, + object userdata = None) except? -1: + if self.with_userdata == 1: + user_flag = self._prec_solvefn(t, y, yp, r, rvec, z, cj, delta, userdata) + else: + user_flag = self._prec_solvefn(t, y, yp, r, rvec, z, cj, delta) + + if user_flag is None: + user_flag = 0 + return user_flag + +cdef int _prec_solvefn(realtype tt, N_Vector yy, N_Vector yp, N_Vector r, + N_Vector rvec, N_Vector z, realtype cj, + realtype delta, void *auxiliary_data): + """ function with the signature of CVSpilsPrecSolveFn, that calls python function """ + cdef np.ndarray[DTYPE_t, ndim=1] yy_tmp, r_tmp, z_tmp + + aux_data = auxiliary_data + cdef bint parallel_implementation = aux_data.parallel_implementation + + if parallel_implementation: + raise NotImplemented + else: + yy_tmp = aux_data.yy_tmp + yp_tmp = aux_data.yp_tmp + residual_tmp = aux_data.residual_tmp + + if aux_data.r_vec is None: + N = np.alen(yy_tmp) + aux_data.rvec_tmp = np.empty(N, DTYPE) + + if aux_data.z_tmp is None: + N = np.alen(yy_tmp) + aux_data.z_tmp = np.empty(N, DTYPE) + + rvec_tmp = aux_data.rvec_tmp + z_tmp = aux_data.z_tmp + + nv_s2ndarray(yy, yy_tmp) + nv_s2ndarray(yp, yp_tmp) + nv_s2ndarray(r, residual_tmp) + nv_s2ndarray(rvec, rvec_tmp) + nv_s2ndarray(z, z_tmp) + + user_flag = aux_data.prec_solvefn.evaluate(tt, yy_tmp, residual_tmp, + rvec_tmp, z_tmp, cj, delta, + aux_data.user_data) if parallel_implementation: raise NotImplemented else: - #we convert the python jac_tmp array to DslMat of sundials - ndarray2DlsMatd(Jac, jac_tmp) + ndarray2nv_s(z, z_tmp) return user_flag @@ -396,6 +597,8 @@ cdef class IDA_data: self.residual_tmp = np.empty(N, DTYPE) self.jac_tmp = None self.g_tmp = None + self.z_tmp = None + self.rvec_tmp = None cdef class IDA: @@ -434,6 +637,9 @@ cdef class IDA: 'rootfn': None, 'nr_rootfns': 0, 'jacfn': None, + 'precond_type': 'NONE', + 'prec_setupfn': None, + 'prec_solvefn': None, 'err_handler': None, 'err_user_data': None, 'old_api': None, @@ -640,6 +846,27 @@ cdef class IDA: 2.0 - variable has to be positive (i.e. > 0) -1.0 - variable has to be non-positive (i.e. <= 0) -2.0 - variable has to be negative (i.e. < 0) + 'precond_type': + default = None + 'prec_setupfn': + Values: function of class IDA_PrecSetupFunction + Description: + Defines a function that setups the preconditioner on change + of the Jacobian. This function takes as input arguments + current time t, current value of y, flag jok that indicates + whether Jacobian-related data has to be updated, flag jcurPtr + that should be set to True (jcurPtr.value = True) if Jacobian + data was recomputed, parameter gamma and optional userdata. + 'prec_solvefn': + Values: function of class IDA_PrecSolveFunction + Description: + Defines a function that solves the preconditioning problem + P*z = r where P may be a left or right preconditioner + matrix. This function takes as input arguments current time + t, current value of y, right-hand side r, result vector z, + parameters gamma and delta, input flag lr that determines + the flavour of the preconditioner (left = 1, right = 2) and + optional userdata. 'err_handler': Values: function of class IDA_ErrHandler, default = None Description: @@ -961,8 +1188,8 @@ cdef class IDA: self.parallel_implementation = (opts['implementation'].lower() == 'parallel') if self.parallel_implementation: raise ValueError('Error: Parallel implementation not implemented !') - cdef long int N - N = np.alen(y0) + cdef INDEX_TYPE_t N + N = np.alen(y0) if opts['rfn'] is None: raise ValueError('The residual function rfn not assigned ' @@ -1068,6 +1295,22 @@ cdef class IDA: self.aux_data.jac = jac self.aux_data.user_data = opts['user_data'] + prec_setupfn = opts['prec_setupfn'] + if prec_setupfn is not None and not isinstance(prec_setupfn, IDA_PrecSetupFunction): + tmpfun = IDA_WrapPrecSetupFunction() + tmpfun.set_prec_setupfn(prec_setupfn) + prec_setupfn = tmpfun + opts['prec_setupfn'] = tmpfun + self.aux_data.prec_setupfn = prec_setupfn + + prec_solvefn = opts['prec_solvefn'] + if prec_solvefn is not None and not isinstance(prec_solvefn, IDA_PrecSolveFunction): + tmpfun = IDA_WrapPrecSolveFunction() + tmpfun.set_prec_solvefn(prec_solvefn) + prec_solvefn = tmpfun + opts['prec_solvefn'] = tmpfun + self.aux_data.prec_solvefn = prec_solvefn + self._set_runtime_changeable_options(opts, supress_supported_check=True) if flag == IDA_ILL_INPUT: @@ -1097,63 +1340,148 @@ cdef class IDA: linsolver = opts['linsolver'].lower() if linsolver == 'dense': - if self.parallel_implementation: - raise ValueError('Linear solver for dense matrices can be used ' - 'only for serial implementation. For parallel' - ' implementation use ''lapackdense'' instead.') - else: - flag = IDADense(ida_mem, N) - if flag == IDADLS_ILL_INPUT: - raise ValueError('IDADense solver is not compatible with' - ' the current nvector implementation.') - elif flag == IDADLS_MEM_FAIL: - raise MemoryError('IDADense memory allocation error.') - elif linsolver == 'lapackdense': - flag = IDALapackDense(ida_mem, N) + A = SUNDenseMatrix(N, N) + LS = SUNDenseLinearSolver(self.y0, A) + # check if memory was allocated + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + # attach matrix and linear solver to cvode + flag = IDADlsSetLinearSolver(ida_mem, LS, A) if flag == IDADLS_ILL_INPUT: - raise ValueError('IDALapackDense solver is not compatible with' - ' the current nvector implementation.') - elif flag == IDADLS_MEM_FAIL: - raise MemoryError('IDALapackDense memory allocation error.') + raise ValueError('IDADense linear solver setting failed, ' + 'arguments incompatible') + elif flag == IDADLS_MEM_NULL: + raise MemoryError('IDADense linear solver memory allocation error.') + elif flag != IDADLS_SUCCESS: + raise ValueError('IDADlsSetLinearSolver failed with code {}' + .format(flag)) elif linsolver == 'band': - if self.parallel_implementation: - raise ValueError('Linear solver for band matrices can be used' - 'only for serial implementation. For parallel' - ' implementation use ''lapackband'' instead.') - else: - flag = IDABand(ida_mem, N, opts['uband'], - opts['lband']) - if flag == IDADLS_ILL_INPUT: - raise ValueError('IDABand solver is not compatible' - ' with the current nvector implementation' - ' or bandwith outside range.') - elif flag == IDADLS_MEM_FAIL: - raise MemoryError('IDABand memory allocation error.') - elif linsolver == 'lapackband': - flag = IDALapackBand(ida_mem, N, opts['uband'], - opts['lband']) + A = SUNBandMatrix(N, opts['uband'], opts['lband'], + opts['uband'] + opts['lband']); + LS = SUNBandLinearSolver(self.y0, A); + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + flag = IDADlsSetLinearSolver(ida_mem, LS, A) if flag == IDADLS_ILL_INPUT: - raise ValueError('IDALapackBand solver is not compatible' - ' with the current nvector implementation' - ' or bandwith outside range.') - elif flag == IDADLS_MEM_FAIL: - raise MemoryError('IDALapackBand memory allocation error.') + raise ValueError('IDABand linear solver setting failed, ' + 'arguments incompatible') + elif flag == IDADLS_MEM_NULL: + raise MemoryError('IDABand linear solver memory allocation error.') + elif flag != IDADLS_SUCCESS: + raise ValueError('IDADlsSetLinearSolver failed with code {}' + .format(flag)) elif ((linsolver == 'spgmr') or (linsolver == 'spbcg') or (linsolver == 'sptfqmr')): maxl = opts['maxl'] + precond_type = opts['precond_type'].lower() + if precond_type == 'none': + pretype = PREC_NONE + elif precond_type == 'left': + pretype = PREC_LEFT + elif precond_type == 'right': + pretype = PREC_RIGHT + elif precond_type == 'both': + pretype = PREC_BOTH + else: + raise ValueError('LinSolver::Precondition: Unknown type: %s' + % opts['precond_type']) + + if linsolver == 'spgmr': - flag = IDASpgmr(ida_mem, maxl) - elif linsolver == 'spbcg': - flag = IDASpbcg(ida_mem, maxl) + LS = SUNSPGMR(self.y0, pretype, maxl); + if LS == NULL: + raise ValueError('Could not allocate linear solver') + elif linsolver == 'spbcgs': + LS = SUNSPBCGS(self.y0, pretype, maxl); + if LS == NULL: + raise ValueError('Could not allocate linear solver') + elif linsolver == 'sptfqmr': + LS = SUNSPTFQMR(self.y0, pretype, maxl); + if LS == NULL: + raise ValueError('Could not allocate linear solver') else: - flag = IDASptfqmr(ida_mem, maxl) + raise ValueError('Given linsolver {} not implemented in odes'.format(linsolver)) + + flag = IDASpilsSetLinearSolver(ida_mem, LS); + if flag == IDASPILS_MEM_NULL: + raise MemoryError('IDA memory was NULL') + elif flag == IDASPILS_ILL_INPUT: + raise MemoryError('linear solver memory was NULL') + elif flag != IDASPILS_SUCCESS: + raise ValueError('CVSpilsSetLinearSolver failed with code {}' + .format(flag)) + # TODO: make option for the Gram-Schmidt orthogonalization + #flag = SUNSPGMRSetGSType(LS, gstype); + + # TODO make option + #flag = IDASpilsSetEpsLin(cvode_mem, DELT); + + if self.aux_data.prec_solvefn: + if self.aux_data.prec_setupfn: + flag = IDASpilsSetPreconditioner(ida_mem, _prec_setupfn, + _prec_solvefn) + else: + flag = IDASpilsSetPreconditioner(ida_mem, NULL, _prec_solvefn) + if flag == IDASPILS_MEM_NULL: + raise ValueError('LinSolver: The cvode mem pointer is NULL.') + elif flag == IDASPILS_LMEM_NULL: + raise ValueError('LinSolver: The cvspils linear solver has ' + 'not been initialized.') + elif flag != IDASPILS_SUCCESS: + raise ValueError('IDASpilsSetPreconditioner failed with code {}' + .format(flag)) + else: + IF SUNDIALS_BLAS_LAPACK: + if linsolver == 'lapackdense': + A = SUNDenseMatrix(N, N) + LS = SUNLapackDense(self.y0, A) + # check if memory was allocated + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + # attach matrix and linear solver to cvode + flag = IDADlsSetLinearSolver(ida_mem, LS, A) + if flag == IDADLS_ILL_INPUT: + raise ValueError('IDADense linear solver setting failed, ' + 'arguments incompatible') + elif flag == IDADLS_MEM_NULL: + raise MemoryError('IDADense linear solver memory allocation error.') + elif flag != IDADLS_SUCCESS: + raise ValueError('IDADlsSetLinearSolver failed with code {}' + .format(flag)) + elif linsolver == 'lapackband': + A = SUNBandMatrix(N, opts['uband'], opts['lband'], + opts['uband'] + opts['lband']); + LS = SUNLapackBand(self.y0, A) + if (A == NULL or LS == NULL): + raise ValueError('Could not allocate matrix or linear solver') + flag = IDADlsSetLinearSolver(ida_mem, LS, A) + if flag == IDADLS_ILL_INPUT: + raise ValueError('IDABand linear solver setting failed, ' + 'arguments incompatible') + elif flag == IDADLS_MEM_NULL: + raise MemoryError('IDABand linear solver memory allocation error.') + elif flag != IDADLS_SUCCESS: + raise ValueError('IDADlsSetLinearSolver failed with code {}' + .format(flag)) + else: + raise ValueError('LinSolver: Unknown solver type: %s' + % opts['linsolver']) + ELSE: + raise ValueError('LinSolver: Unknown solver type: %s' + % opts['linsolver']) - if flag == IDASPILS_MEM_FAIL: - raise MemoryError('LinSolver:IDASpils memory allocation error.') if (linsolver in ['dense', 'lapackdense']) and self.aux_data.jac: - IDADlsSetDenseJacFn(ida_mem, _jacdense) + self.aux_data.jac_tmp = np.empty((np.alen(y0), np.alen(y0)), DTYPE) + flag = IDADlsSetJacFn(ida_mem, _jacdense) + if flag == IDADLS_MEM_NULL: + raise MemoryError('IDA Memory NULL.') + if flag == IDADLS_LMEM_NULL: + raise ValueError('IDA linear solver memory NULL') + elif flag != IDADLS_SUCCESS: + raise ValueError('IDADlsSetJacFn failed with code {}' + .format(flag)) # Constraints constraints_idx = opts['constraints_idx'] @@ -1335,9 +1663,9 @@ cdef class IDA: self._init_step(t0, y0, yp0) return - cdef long int N - N = np.alen(y0) - Np = np.alen(yp0) + cdef INDEX_TYPE_t N + N = np.alen(y0) + Np = np.alen(yp0) if N == self.N and Np == N: self.y0 = N_VMake_Serial(N, y0.data) self.yp0 = N_VMake_Serial(N, yp0.data) diff --git a/scikits/odes/tests/test_dae.py b/scikits/odes/tests/test_dae.py index 52b65c1..fb02ed3 100644 --- a/scikits/odes/tests/test_dae.py +++ b/scikits/odes/tests/test_dae.py @@ -48,7 +48,7 @@ def _do_problem(self, problem, integrator, old_api=True, **integrator_params): def test_ddaspk(self): """Check the ddaspk solver""" - for problem_cls in PROBLEMS: + for problem_cls in PROBLEMS_DDASPK: problem = problem_cls() self._do_problem(problem, 'ddaspk', **problem.ddaspk_pars) @@ -60,13 +60,13 @@ def test_lsodi(self): def test_ida_old_api(self): """Check the ida solver""" - for problem_cls in PROBLEMS: + for problem_cls in PROBLEMS_IDA: problem = problem_cls() self._do_problem(problem, 'ida', old_api=True, **problem.ida_pars) def test_ida(self): """Check the ida solver""" - for problem_cls in PROBLEMS: + for problem_cls in PROBLEMS_IDA: problem = problem_cls() self._do_problem(problem, 'ida', old_api=False, **problem.ida_pars) @@ -152,6 +152,14 @@ def jac(self, t, y, yp, cj, jac): jac[0][0] = self.m*cj_in ;jac[0][1] = self.k jac[1][0] = -1 ;jac[1][1] = cj_in; +class SimpleOscillatorJacIDA(SimpleOscillator): + def jac(self, t, y, yp, residual, cj, jac): + """Jacobian[i,j] is dRES(i)/dY(j) + CJ*dRES(i)/dYPRIME(j)""" + jc = zeros((len(y), len(y)), DTYPE) + cj_in = cj + jac[0][0] = self.m*cj_in ;jac[0][1] = self.k + jac[1][0] = -1 ;jac[1][1] = cj_in; + class StiffVODECompare(DAE): r""" We create a stiff problem, obtain the vode solution, and compare with @@ -183,7 +191,7 @@ def jac_vode(self, t, y): def f_cvode(self, t, y, ydot): ydot[:] = self.f_vode(t, y) - def jac_cvode(self, t, y, J): + def jac_cvode(self, t, y, fy, J): J[:,:] = self.jac_vode(t, y) def __init__(self): @@ -248,8 +256,10 @@ def verify(self, y, yp, t): return ( allclose(self.sol, y, atol=self.atol, rtol=self.rtol) and allclose(self.sol2, y, atol=self.atol, rtol=self.rtol) ) -PROBLEMS = [SimpleOscillator, StiffVODECompare, +PROBLEMS_DDASPK = [SimpleOscillator, StiffVODECompare, SimpleOscillatorJac ] +PROBLEMS_IDA = [SimpleOscillator, StiffVODECompare, + SimpleOscillatorJacIDA ] PROBLEMS_LSODI = [SimpleOscillator, StiffVODECompare] #------------------------------------------------------------------------------ diff --git a/scikits/odes/tests/test_odeint.py b/scikits/odes/tests/test_odeint.py index 05c5362..04174f6 100644 --- a/scikits/odes/tests/test_odeint.py +++ b/scikits/odes/tests/test_odeint.py @@ -106,7 +106,7 @@ def f(self, t, z, rhs): rhs[:] = 1j*z return 0 - def jac(self, t, z, J): + def jac(self, t, z, fz, J): J[:,:] = 1j*eye(5) return 0 @@ -151,7 +151,7 @@ def f(self, t, z, rhs): -lmbd[2]*z[2] + lmbd[1]*z[1]]) return 0 - def jac(self, t, z, J): + def jac(self, t, z, fz, J): # The full Jacobian is # # [-lmbd[0] 0 0 ] @@ -253,14 +253,14 @@ def test_odeint_trivial_time(): c = array([[-205, 0.01, 0.00, 0.0], [0.1, -2.50, 0.02, 0.0], - [1e-3, 0.01, -2.0, 0.01], + [0.00, 0.01, -2.0, 0.01], [0.00, 0.00, 0.1, -1.0]]) cband1 = array([ - [0. , 0., 0.00, 0.0 ], [0. , 0.01, 0.02, 0.01], [-205, -2.50, -2.0, -1.0 ], - [0.1, 0.01, 0.1, 0.0 ]]) + [0.1, 0.01, 0.1, 0.0 ], + [0. , 0., 0.00, 0.0 ]]) cband2 = array([[0, 0.01, 0.02, 0.01], [-250, -2.50, -2., -1.], [0.1, 0.01, 0.1, 0.]]) @@ -273,18 +273,18 @@ def func(t, y, rhs): rhs[:] = c.dot(y) return 0 - def jac(t, y, J): + def jac(t, y, fy, J): J[:,:] = c return 0 - def bjac1_rows(t, y, J): + def bjac1_rows(t, y, fy, J): J[:,:] = cband1 - return jac + return 0 - def bjac2_rows(t, y, J): + def bjac2_rows(t, y, fy, J): # define BAND_ELEM(A,i,j) ((A->cols)[j][(i)-(j)+(A->s_mu)]) J[:,:] = cband2 - return jac + return 0 y0 = np.ones(4) t = np.array([0, 5, 10, 100]) @@ -354,7 +354,7 @@ def bjac2_rows(t, y, J): linsolver='spgmr') sol11 = odeint(func, t, y0, atol=1e-13, rtol=1e-11, max_steps=10000, - linsolver='spbcg') + linsolver='spbcgs') sol12 = odeint(func, t, y0, atol=1e-13, rtol=1e-11, max_steps=10000, linsolver='sptfqmr') @@ -390,11 +390,11 @@ def bad2(t, x, rhs): rhs[:] = "foo" return 0 - def bad_jac1(t, x, J): + def bad_jac1(t, x, fx, J): J[:,:] = 1.0/0 return 0 - def bad_jac2(t, x, J): + def bad_jac2(t, x, fx, J): J[:,:] = [["foo"]] return 0 @@ -402,7 +402,7 @@ def sys2d(t, x, rhs): rhs[:] = [-100*x[0], -0.1*x[1]] return 0 - def sys2d_bad_jac(t, x, J): + def sys2d_bad_jac(t, x, fx, J): J[:,:] = [[1.0/0, 0], [0, -0.1]] return 0 @@ -427,7 +427,7 @@ def sys1(t, x, rhs): rhs[:] = -100*x return 0 - def badjac(t, x, J): + def badjac(t, x, fx, J): J[:,:] = [[0, 0, 0]] return 0 diff --git a/scikits/odes/tests/test_user_return_vals_cvode.py b/scikits/odes/tests/test_user_return_vals_cvode.py index 1d15b49..8af4b0e 100644 --- a/scikits/odes/tests/test_user_return_vals_cvode.py +++ b/scikits/odes/tests/test_user_return_vals_cvode.py @@ -53,27 +53,27 @@ def root_error_late(t, y, g): def root_error_immediate(t, y, g): return 1 -def normal_jac(t, y, J): +def normal_jac(t, y, fy, J): J[0][0] = 0 -def jac_with_return(t, y, J): +def jac_with_return(t, y, fy, J): J[0][0] = 0 return 0 -def jac_problem_late(t, y, J): +def jac_problem_late(t, y, fy, J): J[0][0] = 1 if t > 0: return 1 -def jac_problem_immediate(t, y, J): +def jac_problem_immediate(t, y, fy, J): return 1 -def jac_error_late(t, y, J): +def jac_error_late(t, y, fy, J): J[0][0] = 1 if t > 0: return -1 -def jac_error_immediate(t, y, J): +def jac_error_immediate(t, y, fy, J): return -1 diff --git a/scikits/odes/tests/test_user_return_vals_ida.py b/scikits/odes/tests/test_user_return_vals_ida.py index 3c0b006..3fc9fd8 100644 --- a/scikits/odes/tests/test_user_return_vals_ida.py +++ b/scikits/odes/tests/test_user_return_vals_ida.py @@ -53,27 +53,27 @@ def root_error_late(t, y, y_dot, g): def root_error_immediate(t, y, y_dot, g): return 1 -def normal_jac(t, y, ydot, cj, J): +def normal_jac(t, y, ydot, residual, cj, J): J[0][0] = cj -def jac_with_return(t, y, ydot, cj, J): +def jac_with_return(t, y, ydot, residual, cj, J): J[0][0] = cj return 0 -def jac_problem_late(t, y, ydot, cj, J): +def jac_problem_late(t, y, ydot, residual, cj, J): J[0][0] = cj + 1 if t > 0: return 1 -def jac_problem_immediate(t, y, ydot, cj, J): +def jac_problem_immediate(t, y, ydot, residual, cj, J): return 1 -def jac_error_late(t, y, ydot, cj, J): +def jac_error_late(t, y, ydot, residual, cj, J): J[0][0] = cj + 1 if t > 0: return -1 -def jac_error_immediate(t, y, ydot, cj, J): +def jac_error_immediate(t, y, ydot, residual, cj, J): return -1 class TestIdaReturn(TestCase): diff --git a/setup_build.py b/setup_build.py index 59207cc..60576c0 100644 --- a/setup_build.py +++ b/setup_build.py @@ -21,6 +21,26 @@ def write_pxi(filename, definitions): return filename +def check_macro_def(cmd, symbol, headers=None, include_dirs=None): + """ + Based on numpy.distutils.command.config:config.check_macro_true, checks if + macro is defined or not + """ + cmd._check_compiler() + body = """ +int main(void) +{ +#ifdef %s +#else +#error undefined macro +#endif + ; + return 0; +}""" % (symbol,) + + return cmd.try_compile(body, headers, include_dirs) + + def get_sundials_config_pxi(include_dirs, dist): """ Create pxi file containing some of sundials build config @@ -33,6 +53,8 @@ def get_sundials_config_pxi(include_dirs, dist): BASE_PATH = join('scikits', 'odes', 'sundials') config_cmd = dist.get_command_obj("config") + + # Get float type if config_cmd.check_macro_true( "SUNDIALS_DOUBLE_PRECISION", headers=[SUNDIALS_CONFIG_H], include_dirs=include_dirs @@ -56,10 +78,46 @@ def get_sundials_config_pxi(include_dirs, dist): SUNDIALS_FLOAT_TYPE = '"double"' info("Failed to find sundials precision, falling back to double...") - return write_pxi(join(BASE_PATH, "sundials_config.pxi"), - dict(SUNDIALS_FLOAT_TYPE=SUNDIALS_FLOAT_TYPE), + # Get index (int) type + if config_cmd.check_macro_true( + "SUNDIALS_INT32_T", headers=[SUNDIALS_CONFIG_H], + include_dirs=include_dirs + ): + SUNDIALS_INDEX_TYPE = '"int32"' + info("Found sundials built with int32.") + elif config_cmd.check_macro_true( + "SUNDIALS_INT64_T", headers=[SUNDIALS_CONFIG_H], + include_dirs=include_dirs + ): + SUNDIALS_INDEX_TYPE = '"int64"' + info("Found sundials built with int64.") + else: + # fall back to int64 + SUNDIALS_INDEX_TYPE = '"int64"' + info("Failed to find sundials index type, falling back to int64...") + + # Check for blas/lapack + if check_macro_def( + config_cmd, + "SUNDIALS_BLAS_LAPACK", headers=[SUNDIALS_CONFIG_H], + include_dirs=include_dirs + ): + has_lapack = True + else: + has_lapack = False + + cfg = dict( + float_type = SUNDIALS_FLOAT_TYPE, + index_type = SUNDIALS_INDEX_TYPE, + has_lapack = has_lapack, ) + return write_pxi(join(BASE_PATH, "sundials_config.pxi"), dict( + SUNDIALS_FLOAT_TYPE=SUNDIALS_FLOAT_TYPE, + SUNDIALS_INDEX_TYPE=SUNDIALS_INDEX_TYPE, + SUNDIALS_BLAS_LAPACK=str(has_lapack), + )), cfg + class build_ext(_build_ext): """ @@ -131,8 +189,40 @@ def _get_cython_ext(self): except ImportError: info("pkgconfig module not found, using preset paths") + sundials_pxi, cfg = get_sundials_config_pxi(SUNDIALS_INCLUDE_DIRS, + self.distribution) + + has_lapack = cfg['has_lapack'] + if not SUNDIALS_LIBRARIES: + # This is where to put N_vector codes (currently only serial is + # supported) SUNDIALS_LIBRARIES.append('sundials_nvecserial') + # SUNDIALS_LIBRARIES.append('sundials_nvecopenmp') + # SUNDIALS_LIBRARIES.append('sundials_nvecparallel') + # SUNDIALS_LIBRARIES.append('sundials_nvecparhyp') + # SUNDIALS_LIBRARIES.append('sundials_nvecpetsc') + # SUNDIALS_LIBRARIES.append('sundials_nvecpthreads') + + # This is where to put SUNLinearSolver codes (klu not supported + # yet) + if has_lapack: + SUNDIALS_LIBRARIES.append('sundials_sunlinsollapackband') + SUNDIALS_LIBRARIES.append('sundials_sunlinsollapackdense') + + SUNDIALS_LIBRARIES.append('sundials_sunlinsolband') + SUNDIALS_LIBRARIES.append('sundials_sunlinsoldense') + SUNDIALS_LIBRARIES.append('sundials_sunlinsolpcg') + SUNDIALS_LIBRARIES.append('sundials_sunlinsolspbcgs') + SUNDIALS_LIBRARIES.append('sundials_sunlinsolspfgmr') + SUNDIALS_LIBRARIES.append('sundials_sunlinsolspgmr') + SUNDIALS_LIBRARIES.append('sundials_sunlinsolsptfqmr') + # SUNDIALS_LIBRARIES.append('sundials_sunlinsolklu') + + # This is where to put SUNMatrix codes + SUNDIALS_LIBRARIES.append('sundials_sunmatrixband') + SUNDIALS_LIBRARIES.append('sundials_sunmatrixdense') + SUNDIALS_LIBRARIES.append('sundials_sunmatrixsparse') if not IDA_LIBRARIES: IDA_LIBRARIES.append('sundials_ida') @@ -141,20 +231,18 @@ def _get_cython_ext(self): CVODE_LIBRARIES.append('sundials_cvode') - try: + if has_lapack: lapack_opt = get_info('lapack_opt', notfound_action=2) if lapack_opt: SUNDIALS_INCLUDE_DIRS.extend(lapack_opt.get('include_dirs',[])) SUNDIALS_LIBRARY_DIRS.extend(lapack_opt.get('library_dirs',[])) SUNDIALS_LIBRARIES.extend(lapack_opt.get('libraries',[])) - use_lapack = True + info('Found LAPACK paths via lapack_opt ...') else: - raise ValueError - info('Found LAPACK paths via lapack_opt ...') - except: - info('LAPACK was not detected, disabling sundials solvers') - return [] + info('LAPACK was not found, but SUNDIALS compiled against ' + 'lapack, check your numpy installation' + ) CVODE_LIBRARIES.extend(SUNDIALS_LIBRARIES) IDA_LIBRARIES.extend(SUNDIALS_LIBRARIES) @@ -163,14 +251,13 @@ def _get_cython_ext(self): CVODE_LIBRARY_DIRS.extend(SUNDIALS_LIBRARY_DIRS) IDA_LIBRARY_DIRS.extend(SUNDIALS_LIBRARY_DIRS) - sundials_pxi = get_sundials_config_pxi(SUNDIALS_INCLUDE_DIRS, - self.distribution) - return [ Extension( base_module + '.' + "common_defs", sources = [join(base_path, 'common_defs.pyx')], include_dirs=SUNDIALS_INCLUDE_DIRS, + library_dirs=SUNDIALS_LIBRARY_DIRS, + libraries=SUNDIALS_LIBRARIES, ), Extension( base_module + '.' + "cvode", diff --git a/tox.ini b/tox.ini index 8c5c807..9b3efa9 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist = py27,py33,py34,py35,py36,check-manifest,checkreadme,docs +envlist = py27,py34,py35,py36,py37,check-manifest,checkreadme,docs #skipsdist=True [testenv] @@ -10,13 +10,12 @@ passenv= LIBRARY_PATH CPATH PIP_VERBOSE + PYTHONFAULTHANDLER deps = - py33: numpy<1.11 - py{27,34,35,36}: numpy + numpy cython nose - py33: pytest<3.3 - py{27,34,35,36}: pytest + pytest wheel commands = env