diff --git a/build-aux/subst-config-vals.in.sh b/build-aux/subst-config-vals.in.sh --- a/build-aux/subst-config-vals.in.sh +++ b/build-aux/subst-config-vals.in.sh @@ -95,9 +95,6 @@ HDF5_LDFLAGS="@HDF5_LDFLAGS@" HDF5_LIBS="@HDF5_LIBS@" includedir="@includedir@" -KLU_CPPFLAGS="@KLU_CPPFLAGS@" -KLU_LDFLAGS="@KLU_LDFLAGS@" -KLU_LIBS="@KLU_LIBS@" LAPACK_LIBS="@LAPACK_LIBS@" LD_CXX="@LD_CXX@" LDFLAGS="@LDFLAGS@" @@ -156,12 +153,6 @@ SH_LD="@SH_LD@" SH_LDFLAGS="@SH_LDFLAGS@" STATIC_LIBS="@STATIC_LIBS@" -SUNDIALS_IDA_CPPFLAGS="@SUNDIALS_IDA_CPPFLAGS@" -SUNDIALS_IDA_LDFLAGS="@SUNDIALS_IDA_LDFLAGS@" -SUNDIALS_IDA_LIBS="@SUNDIALS_IDA_LIBS@" -SUNDIALS_NVECSERIAL_CPPFLAGS="@SUNDIALS_NVECSERIAL_CPPFLAGS@" -SUNDIALS_NVECSERIAL_LDFLAGS="@SUNDIALS_NVECSERIAL_LDFLAGS@" -SUNDIALS_NVECSERIAL_LIBS="@SUNDIALS_NVECSERIAL_LIBS@" SUITESPARSE_CONFIG_LIBS="@SUITESPARSE_CONFIG_LIBS@" TERM_LIBS="@TERM_LIBS@" UMFPACK_CPPFLAGS="@UMFPACK_CPPFLAGS@" @@ -258,9 +249,6 @@ $SED -e "s|%OCTAVE_CONF_HDF5_LDFLAGS%|\"${HDF5_LDFLAGS}\"|" \ -e "s|%OCTAVE_CONF_HDF5_LIBS%|\"${HDF5_LIBS}\"|" \ -e "s|%OCTAVE_CONF_INCLUDEDIR%|\"${includedir}\"|" \ - -e "s|%OCTAVE_CONF_KLU_CPPFLAGS%|\"${KLU_CPPFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_KLU_LDFLAGS%|\"${KLU_LDFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_KLU_LIBS%|\"${KLU_LIBS}\"|" \ -e "s|%OCTAVE_CONF_LAPACK_LIBS%|\"${LAPACK_LIBS}\"|" \ -e "s|%OCTAVE_CONF_LD_CXX%|\"${LD_CXX}\"|" \ -e "s|%OCTAVE_CONF_LDFLAGS%|\"${LDFLAGS}\"|" \ @@ -319,12 +307,6 @@ -e "s|%OCTAVE_CONF_SH_LD%|\"${SH_LD}\"|" \ -e "s|%OCTAVE_CONF_SH_LDFLAGS%|\"${SH_LDFLAGS}\"|" \ -e "s|%OCTAVE_CONF_STATIC_LIBS%|\"${STATIC_LIBS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_IDA_CPPFLAGS%|\"${SUNDIALS_IDA_CPPFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_IDA_LDFLAGS%|\"${SUNDIALS_IDA_LDFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_IDA_LIBS%|\"${SUNDIALS_IDA_LIBS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_NVECSERIAL_CPPFLAGS%|\"${SUNDIALS_NVECSERIAL_CPPFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_NVECSERIAL_LDFLAGS%|\"${SUNDIALS_NVECSERIAL_LDFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_NVECSERIAL_LIBS%|\"${SUNDIALS_NVECSERIAL_LIBS}\"|" \ -e "s|%OCTAVE_CONF_SUITESPARSE_CONFIG_LIBS%|\"${SUITESPARSE_CONFIG_LIBS}\"|" \ -e "s|%OCTAVE_CONF_TERM_LIBS%|\"${TERM_LIBS}\"|" \ -e "s|%OCTAVE_CONF_UMFPACK_CPPFLAGS%|\"${UMFPACK_CPPFLAGS}\"|" \ diff --git a/build-aux/subst-cross-config-vals.in.sh b/build-aux/subst-cross-config-vals.in.sh --- a/build-aux/subst-cross-config-vals.in.sh +++ b/build-aux/subst-cross-config-vals.in.sh @@ -107,9 +107,6 @@ HDF5_LDFLAGS="@HDF5_LDFLAGS@" HDF5_LIBS="@HDF5_LIBS@" includedir="@includedir@" -KLU_CPPFLAGS="@KLU_CPPFLAGS@" -KLU_LDFLAGS="@KLU_LDFLAGS@" -KLU_LIBS="@KLU_LIBS@" LAPACK_LIBS="@LAPACK_LIBS@" LD_CXX="@LD_CXX@" LDFLAGS="@LDFLAGS@" @@ -160,12 +157,6 @@ SH_LD="@SH_LD@" SH_LDFLAGS="@SH_LDFLAGS@" STATIC_LIBS="@STATIC_LIBS@" -SUNDIALS_IDA_CPPFLAGS="@SUNDIALS_IDA_CPPFLAGS@" -SUNDIALS_IDA_LDFLAGS="@SUNDIALS_IDA_LDFLAGS@" -SUNDIALS_IDA_LIBS="@SUNDIALS_IDA_LIBS@" -SUNDIALS_NVECSERIAL_CPPFLAGS="@SUNDIALS_NVECSERIAL_CPPFLAGS@" -SUNDIALS_NVECSERIAL_LDFLAGS="@SUNDIALS_NVECSERIAL_LDFLAGS@" -SUNDIALS_NVECSERIAL_LIBS="@SUNDIALS_NVECSERIAL_LIBS@" SUITESPARSE_CONFIG_LIBS="@SUITESPARSE_CONFIG_LIBS@" TERM_LIBS="@TERM_LIBS@" UMFPACK_CPPFLAGS="@UMFPACK_CPPFLAGS@" @@ -262,9 +253,6 @@ $SED -e "s|%OCTAVE_CONF_HDF5_LDFLAGS%|\"${HDF5_LDFLAGS}\"|" \ -e "s|%OCTAVE_CONF_HDF5_LIBS%|\"${HDF5_LIBS}\"|" \ -e "s|%OCTAVE_CONF_INCLUDEDIR%|\"${includedir}\"|" \ - -e "s|%OCTAVE_CONF_KLU_CPPFLAGS%|\"${KLU_CPPFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_KLU_LDFLAGS%|\"${KLU_LDFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_KLU_LIBS%|\"${KLU_LIBS}\"|" \ -e "s|%OCTAVE_CONF_LAPACK_LIBS%|\"${LAPACK_LIBS}\"|" \ -e "s|%OCTAVE_CONF_LD_CXX%|\"${LD_CXX}\"|" \ -e "s|%OCTAVE_CONF_LDFLAGS%|\"${LDFLAGS}\"|" \ @@ -323,12 +311,6 @@ -e "s|%OCTAVE_CONF_SH_LD%|\"${SH_LD}\"|" \ -e "s|%OCTAVE_CONF_SH_LDFLAGS%|\"${SH_LDFLAGS}\"|" \ -e "s|%OCTAVE_CONF_STATIC_LIBS%|\"${STATIC_LIBS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_IDA_CPPFLAGS%|\"${SUNDIALS_IDA_CPPFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_IDA_LDFLAGS%|\"${SUNDIALS_IDA_LDFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_IDA_LIBS%|\"${SUNDIALS_IDA_LIBS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_NVECSERIAL_CPPFLAGS%|\"${SUNDIALS_NVECSERIAL_CPPFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_NVECSERIAL_LDFLAGS%|\"${SUNDIALS_NVECSERIAL_LDFLAGS}\"|" \ - -e "s|%OCTAVE_CONF_SUNDIALS_NVECSERIAL_LIBS%|\"${SUNDIALS_NVECSERIAL_LIBS}\"|" \ -e "s|%OCTAVE_CONF_SUITESPARSE_CONFIG_LIBS%|\"${SUITESPARSE_CONFIG_LIBS}\"|" \ -e "s|%OCTAVE_CONF_TERM_LIBS%|\"${TERM_LIBS}\"|" \ -e "s|%OCTAVE_CONF_UMFPACK_CPPFLAGS%|\"${UMFPACK_CPPFLAGS}\"|" \ diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -1478,34 +1478,6 @@ FFT_DIR="fftpack" AC_SUBST(FFT_DIR) -### Check for KLU library and header. -OCTAVE_CHECK_LIB(klu, KLU, - [KLU library not found. This will result in some lack of functionality for sparse matrices], - [suitesparse/klu.h ufsparse/klu.h klu/klu.h klu.h], [klu_solve], - [], [don't use KLU library, disable some sparse matrix functionality]) - -### Check for SUNDIALS_IDA library and header. -save_CPPFLAGS="$CPPFLAGS" -save_LIBS="$LIBS" -CPPFLAGS="$KLU_CPPFLAGS $CPPFLAGS" -LIBS="$KLU_LDFLAGS $KLU_LIBS $LIBS" -OCTAVE_CHECK_LIB(sundials_ida, SUNDIALS_IDA, - [SUNDIALS_IDA library configured with klu not found. Solvers ode15i and ode15s will be disabled.], - [ida/ida.h ida.h], [IDAKLU], [], [], - [warn_sundials_ida= - OCTAVE_CHECK_SIZEOF_SUNDIALS_IDA_REALTYPE]) -LIBS="$save_LIBS" -CPPFLAGS="$save_CPPFLAGS" - -## Additional check for sundials_nvecserial library -OCTAVE_CHECK_LIB(sundials_nvecserial, SUNDIALS_NVECSERIAL, - [SUNDIALS_NVECSERIAL library not found. Solvers ode15i and ode15s will be disabled.], - [nvector/nvector_serial.h nvector_serial.h ], [N_VNew_Serial],[],[],[]) - -if test -n "$SUNDIALS_IDA_LIBS" && test -n "$SUNDIALS_NVECSERIAL_LIBS" && test $octave_cv_sizeof_ida_double = yes; then - AC_DEFINE(HAVE_SUNDIALS, 1, [Define if you have both sundials_ida and sundials_nvecserial libraries.]) -fi - ### Check for GLPK library and header. save_CPPFLAGS="$CPPFLAGS" @@ -2979,12 +2951,6 @@ FFTW_XLIBS="$FFTW3_LIBS $FFTW3F_LIBS" -SUNDIALS_XCPPFLAGS="$SUNDIALS_IDA_CPPFLAGS $SUNDIALS_NVECSERIAL_CPPFLAGS $KLU_CPPFLAGS" - -SUNDIALS_XLDFLAGS="$SUNDIALS_IDA_LDFLAGS $SUNDIALS_NVECSERIAL_LDFLAGS $KLU_LDFLAGS" - -SUNDIALS_XLIBS="$SUNDIALS_IDA_LIBS $SUNDIALS_NVECSERIAL_LIBS $KLU_LIBS" - SPARSE_XCPPFLAGS="$CHOLMOD_CPPFLAGS $UMFPACK_CPPFLAGS $AMD_CPPFLAGS $CAMD_CPPFLAGS $COLAMD_CPPFLAGS $CCOLAMD_CPPFLAGS $CXSPARSE_CPPFLAGS" SPARSE_XLDFLAGS="$CHOLMOD_LDFLAGS $UMFPACK_LDFLAGS $AMD_LDFLAGS $CAMD_LDFLAGS $COLAMD_LDFLAGS $CCOLAMD_LDFLAGS $CXSPARSE_LDFLAGS" @@ -2995,10 +2961,6 @@ AC_SUBST(FFTW_XLDFLAGS) AC_SUBST(FFTW_XLIBS) -AC_SUBST(SUNDIALS_XCPPFLAGS) -AC_SUBST(SUNDIALS_XLDFLAGS) -AC_SUBST(SUNDIALS_XLIBS) - AC_SUBST(SPARSE_XCPPFLAGS) AC_SUBST(SPARSE_XLDFLAGS) AC_SUBST(SPARSE_XLIBS) @@ -3007,9 +2969,9 @@ AC_SUBST(GNULIB_LINK_DEPS) -LIBOCTAVE_LINK_DEPS="$GNULIB_LINK_DEPS $CURL_LIBS $SPARSE_XLIBS $ARPACK_LIBS $QRUPDATE_LIBS $FFTW_XLIBS $SUNDIALS_XLIBS $LAPACK_LIBS $BLAS_LIBS $READLINE_LIBS $TERM_LIBS $LIBGLOB $PCRE_LIBS $DL_LIBS $PTHREAD_LIBS $FLIBS $LIBS" - -LIBOCTAVE_LINK_OPTS="$CURL_LDFLAGS $SPARSE_XLDFLAGS $ARPACK_LDFLAGS $QRUPDATE_LDFLAGS $FFTW_XLDFLAGS $SUNDIALS_XLDFLAGS" +LIBOCTAVE_LINK_DEPS="$GNULIB_LINK_DEPS $CURL_LIBS $SPARSE_XLIBS $ARPACK_LIBS $QRUPDATE_LIBS $FFTW_XLIBS $LAPACK_LIBS $BLAS_LIBS $READLINE_LIBS $TERM_LIBS $LIBGLOB $PCRE_LIBS $DL_LIBS $PTHREAD_LIBS $FLIBS $LIBS" + +LIBOCTAVE_LINK_OPTS="$CURL_LDFLAGS $SPARSE_XLDFLAGS $ARPACK_LDFLAGS $QRUPDATE_LDFLAGS $FFTW_XLDFLAGS" AC_SUBST(LIBOCTAVE_LINK_DEPS) AC_SUBST(LIBOCTAVE_LINK_OPTS) @@ -3020,9 +2982,9 @@ LIBOCTINTERP_LINK_DEPS="$DLDFCN_LIBS" fi -LIBOCTINTERP_LINK_DEPS="$LIBOCTINTERP_LINK_DEPS $FT2_LIBS $HDF5_LIBS $MAGICK_LIBS $Z_LIBS $FFTW_XLIBS $SUNDIALS_XLIBS $OPENGL_LIBS $FONTCONFIG_LIBS $FREETYPE_LIBS $X11_LIBS $CARBON_LIBS $GL2PS_LIBS $LLVM_LIBS $JAVA_LIBS $LAPACK_LIBS" - -LIBOCTINTERP_LINK_OPTS="$FT2_LDFLAGS $HDF5_LDFLAGS $MAGICK_LDFLAGS $Z_LDFLAGS $FFTW_XLDFLAGS $SUNDIALS_XLDFLAGS $LLVM_LDFLAGS" +LIBOCTINTERP_LINK_DEPS="$LIBOCTINTERP_LINK_DEPS $FT2_LIBS $HDF5_LIBS $MAGICK_LIBS $Z_LIBS $FFTW_XLIBS $OPENGL_LIBS $FONTCONFIG_LIBS $FREETYPE_LIBS $X11_LIBS $CARBON_LIBS $GL2PS_LIBS $LLVM_LIBS $JAVA_LIBS $LAPACK_LIBS" + +LIBOCTINTERP_LINK_OPTS="$FT2_LDFLAGS $HDF5_LDFLAGS $MAGICK_LDFLAGS $Z_LDFLAGS $FFTW_XLDFLAGS $LLVM_LDFLAGS" OCTAVE_LINK_DEPS="" OCTAVE_LINK_OPTS="" @@ -3293,9 +3255,6 @@ Java JVM path: $JAVA_LDPATH Java CPPFLAGS: $JAVA_CPPFLAGS Java libraries: $JAVA_LIBS - KLU CPPFLAGS: $KLU_CPPFLAGS - KLU LDFLAGS: $KLU_LDFLAGS - KLU libraries: $KLU_LIBS LAPACK libraries: $LAPACK_LIBS LLVM CPPFLAGS: $LLVM_CPPFLAGS LLVM LDFLAGS: $LLVM_LDFLAGS @@ -3331,12 +3290,6 @@ Sndfile CPPFLAGS: $SNDFILE_CPPFLAGS Sndfile LDFLAGS: $SNDFILE_LDFLAGS Sndfile libraries: $SNDFILE_LIBS - SUNDIALS_IDA CPPFLAGS: $SUNDIALS_IDA_CPPFLAGS - SUNDIALS_IDA LDFLAGS: $SUNDIALS_IDA_LDFLAGS - SUNDIALS_IDA libraries: $SUNDIALS_IDA_LIBS - SUNDIALS_NVECSERIAL CPPFLAGS: $SUNDIALS_NVECSERIAL_CPPFLAGS - SUNDIALS_NVECSERIAL LDFLAGS: $SUNDIALS_NVECSERIAL_LDFLAGS - SUNDIALS_NVECSERIAL libraries: $SUNDIALS_NVECSERIAL_LIBS SuiteSparse config libraries: $SUITESPARSE_CONFIG_LIBS TERM libraries: $TERM_LIBS UMFPACK CPPFLAGS: $UMFPACK_CPPFLAGS diff --git a/libinterp/build-env.h b/libinterp/build-env.h --- a/libinterp/build-env.h +++ b/libinterp/build-env.h @@ -104,9 +104,6 @@ extern const char *HDF5_LDFLAGS; extern const char *HDF5_LIBS; extern const char *INCLUDEDIR; - extern const char *KLU_CPPFLAGS; - extern const char *KLU_LDFLAGS; - extern const char *KLU_LIBS; extern const char *LAPACK_LIBS; extern const char *LDFLAGS; extern const char *LD_CXX; @@ -157,12 +154,6 @@ extern const char *SH_LD; extern const char *SH_LDFLAGS; extern const char *STATIC_LIBS; - extern const char *SUNDIALS_IDA_CPPFLAGS; - extern const char *SUNDIALS_IDA_LDFLAGS; - extern const char *SUNDIALS_IDA_LIBS; - extern const char *SUNDIALS_NVECSERIAL_CPPFLAGS; - extern const char *SUNDIALS_NVECSERIAL_LDFLAGS; - extern const char *SUNDIALS_NVECSERIAL_LIBS; extern const char *SUITESPARSE_CONFIG_LIBS; extern const char *TERM_LIBS; extern const char *UMFPACK_CPPFLAGS; diff --git a/libinterp/build-env.in.cc b/libinterp/build-env.in.cc --- a/libinterp/build-env.in.cc +++ b/libinterp/build-env.in.cc @@ -173,12 +173,6 @@ const char *INCLUDEDIR = %OCTAVE_CONF_INCLUDEDIR%; - const char *KLU_CPPFLAGS = %OCTAVE_CONF_KLU_CPPFLAGS%; - - const char *KLU_LDFLAGS = %OCTAVE_CONF_KLU_LDFLAGS%; - - const char *KLU_LIBS = %OCTAVE_CONF_KLU_LIBS%; - const char *LAPACK_LIBS = %OCTAVE_CONF_LAPACK_LIBS%; const char *LDFLAGS = %OCTAVE_CONF_LDFLAGS%; @@ -281,18 +275,6 @@ const char *SUITESPARSE_CONFIG_LIBS = %OCTAVE_CONF_SUITESPARSE_CONFIG_LIBS%; - const char *SUNDIALS_IDA_CPPFLAGS = %OCTAVE_CONF_SUNDIALS_IDA_CPPFLAGS%; - - const char *SUNDIALS_IDA_LDFLAGS = %OCTAVE_CONF_SUNDIALS_IDA_LDFLAGS%; - - const char *SUNDIALS_IDA_LIBS = %OCTAVE_CONF_SUNDIALS_IDA_LIBS%; - - const char *SUNDIALS_NVECSERIAL_CPPFLAGS = %OCTAVE_CONF_SUNDIALS_NVECSERIAL_CPPFLAGS%; - - const char *SUNDIALS_NVECSERIAL_LDFLAGS = %OCTAVE_CONF_SUNDIALS_NVECSERIAL_LDFLAGS%; - - const char *SUNDIALS_NVECSERIAL_LIBS = %OCTAVE_CONF_SUNDIALS_NVECSERIAL_LIBS%; - const char *TERM_LIBS = %OCTAVE_CONF_TERM_LIBS%; const char *UMFPACK_CPPFLAGS = %OCTAVE_CONF_UMFPACK_CPPFLAGS%; diff --git a/libinterp/dldfcn/__ode15__.cc b/libinterp/dldfcn/__ode15__.cc deleted file mode 100644 --- a/libinterp/dldfcn/__ode15__.cc +++ /dev/null @@ -1,1323 +0,0 @@ -/* - -Copyright (C) 2016 Francesco Faccio - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by -the Free Software Foundation; either version 3 of the License, or (at -your option) any later version. - -Octave is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -General Public License for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -. - -*/ - -#if defined (HAVE_CONFIG_H) -# include "config.h" -#endif - -#if defined (HAVE_IDA_IDA_H) -#include -#include -#include -#include -#endif - -#if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H) -#include -#endif - - -#include -#include -#include - -#include "defun-dld.h" -#include "iostream" -#include "error.h" -#include "ov.h" -#include "ov-struct.h" -#include "errwarn.h" - -#include "algorithm" -#include "iterator" - - -#if defined (HAVE_SUNDIALS) - -namespace octave -{ - class IDA - { - public: - - typedef - ColumnVector (*DAERHSFuncIDA) (const ColumnVector& x, - const ColumnVector& xdot, - realtype t, octave_function *idaf); - - typedef - Matrix (*DAEJacFuncDense) (const ColumnVector& x, - const ColumnVector& xdot, realtype t, - realtype cj, octave_function *idaj); - - typedef - SparseMatrix (*DAEJacFuncSparse) (const ColumnVector& x, - const ColumnVector& xdot, - realtype t, realtype cj, - octave_function *idaj); - - typedef - Matrix (*DAEJacCellDense) (Matrix *dfdy, Matrix *dfdyp, - realtype cj); - - typedef - SparseMatrix (*DAEJacCellSparse) (SparseMatrix *dfdy, - SparseMatrix *dfdyp, realtype cj); - - //Default - IDA (void) - : t0 (0.0), y0 (), yp0 (), havejac (false), havejacfun (false), - havejacsparse (false), mem (NULL), num (), ida_fun (NULL), - ida_jac (NULL), dfdy (NULL), dfdyp (NULL), spdfdy (NULL), - spdfdyp (NULL), fun (0), jacfun (NULL), jacspfun (0), - jacdcell (0), jacspcell (0) - { }; - - - IDA (realtype t, ColumnVector y, ColumnVector yp, - octave_function *ida_fcn, DAERHSFuncIDA daefun) - : t0 (t), y0 (y), yp0 (yp), havejac (false), havejacfun (false), - havejacsparse (false), mem (NULL), num (), ida_fun (ida_fcn), - ida_jac (NULL), dfdy (NULL), dfdyp (NULL), spdfdy (NULL), - spdfdyp (NULL), fun (daefun), jacfun (NULL), jacspfun (0), - jacdcell (0), jacspcell (0) - { }; - - - ~IDA (void) { IDAFree (&mem); } - - IDA& - set_jacobian (octave_function *jac, DAEJacFuncDense j) - { - jacfun = j; - ida_jac = jac; - havejac = true; - havejacfun = true; - havejacsparse = false; - return *this; - } - - IDA& - set_jacobian (octave_function *jac, DAEJacFuncSparse j) - { - jacspfun = j; - ida_jac = jac; - havejac = true; - havejacfun = true; - havejacsparse = true; - return *this; - } - - IDA& - set_jacobian (Matrix *dy, Matrix *dyp, DAEJacCellDense j) - { - jacdcell = j; - dfdy = dy; - dfdyp = dyp; - havejac = true; - havejacfun = false; - havejacsparse = false; - return *this; - } - - IDA& - set_jacobian (SparseMatrix *dy, SparseMatrix *dyp, - DAEJacCellSparse j) - { - jacspcell = j; - spdfdy = dy; - spdfdyp = dyp; - havejac = true; - havejacfun = false; - havejacsparse = true; - return *this; - } - - void - set_userdata (void); - - void - initialize (void); - - static ColumnVector - NVecToCol (N_Vector& v, long int n); - - static N_Vector - ColToNVec (const ColumnVector& data, long int n); - - void - set_up (void); - - void - set_tolerance (ColumnVector& abstol, realtype reltol); - - void - set_tolerance (realtype abstol, realtype reltol); - - static int - resfun (realtype t, N_Vector yy, N_Vector yyp, - N_Vector rr, void *user_data); - - void - resfun_impl (realtype t, N_Vector& yy, - N_Vector& yyp, N_Vector& rr); - - static int - jacdense (long int Neq, realtype t, realtype cj, - N_Vector yy, N_Vector yyp, N_Vector resvec, - DlsMat JJ, void *user_data, N_Vector tempv1, - N_Vector tempv2, N_Vector tempv3) - { - IDA *self = - static_cast (user_data); - - self -> jacdense_impl (Neq, t, cj, yy, yyp, JJ); - return 0; - } - - void - jacdense_impl (long int Neq, realtype t, realtype cj, - N_Vector& yy, N_Vector& yyp, DlsMat& JJ); - - static int - jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp, - N_Vector r, SlsMat Jac, void *user_data, N_Vector tmp1, - N_Vector tmp2, N_Vector tmp3) - { - IDA *self = - static_cast (user_data); - - self -> jacsparse_impl (t, cj, yy, yyp, Jac); - return 0; - } - - void - jacsparse_impl (realtype t, realtype cj, N_Vector& yy, - N_Vector& yyp, SlsMat& Jac); - - void - set_maxstep (realtype maxstep); - - void - set_initialstep (realtype initialstep); - - bool - interpolate (int& cont, Matrix& output, ColumnVector& tout, - int refine, realtype tend, bool haveoutputfcn, - bool haveoutputsel, octave_function *output_fcn, - ColumnVector& outputsel, bool haveeventfunction, - octave_function * event_fcn, ColumnVector& te, - Matrix& ye, ColumnVector& ie, ColumnVector& oldval, - ColumnVector& oldisterminal, ColumnVector& olddir, - int& temp, ColumnVector& yold); - - bool - outputfun (octave_function *output_fcn, bool haveoutputsel, - const ColumnVector& output, realtype tout, realtype tend, - ColumnVector& outputsel, std::string flag); - - - bool - event (octave_function *event_fcn, - ColumnVector& te, Matrix& ye, ColumnVector& ie, - realtype tsol, const ColumnVector& y, std::string flag, - const ColumnVector& yp, ColumnVector& oldval, - ColumnVector& oldisterminal, ColumnVector& olddir, - int cont, int& temp, realtype told, ColumnVector& yold); - - void - set_maxorder (int maxorder); - - octave_value_list - integrate (const int numt, const ColumnVector& tt, - const ColumnVector& y0, const ColumnVector& yp0, - const int refine, bool haverefine, bool haveoutputfcn, - octave_function *output_fcn, bool haveoutputsel, - ColumnVector& outputsel, bool haveeventfunction, - octave_function *event_fcn); - - void - print_stat (void); - - private: - - realtype t0; - ColumnVector y0; - ColumnVector yp0; - bool havejac; - bool havejacfun; - bool havejacsparse; - void *mem; - int num; - octave_function *ida_fun; - octave_function *ida_jac; - Matrix *dfdy; - Matrix *dfdyp; - SparseMatrix *spdfdy; - SparseMatrix *spdfdyp; - DAERHSFuncIDA fun; - DAEJacFuncDense jacfun; - DAEJacFuncSparse jacspfun; - DAEJacCellDense jacdcell; - DAEJacCellSparse jacspcell; - - }; - - int - IDA::resfun (realtype t, N_Vector yy, N_Vector yyp, N_Vector rr, - void *user_data) - { - IDA *self = - static_cast (user_data); - - self -> resfun_impl (t, yy, yyp, rr); - return 0; - } - - void - IDA::resfun_impl (realtype t, N_Vector& yy, - N_Vector& yyp, N_Vector& rr) - { - BEGIN_INTERRUPT_WITH_EXCEPTIONS; - - ColumnVector y = - IDA::NVecToCol (yy, num); - - ColumnVector yp = - IDA::NVecToCol (yyp, num); - - ColumnVector res = - (*fun) (y, yp, t, ida_fun); - - realtype *puntrr = - NV_DATA_S (rr); - - for (octave_idx_type i = 0; i < num; i++) - puntrr [i] = res (i); - - END_INTERRUPT_WITH_EXCEPTIONS; - } - - void - IDA::set_up (void) - { - if (havejacsparse) - { - int flag = - IDAKLU (mem, num, num*num, CSC_MAT); - if (flag != 0) - { - error ("IDAKLU solver not initialized"); - } - flag = IDASlsSetSparseJacFn (mem, IDA::jacsparse); - } - else - { - int flag = - IDADense (mem, num); - if (flag != 0) - { - error ("IDADense solver not initialized"); - } - if (havejac) - { - flag = IDADlsSetDenseJacFn (mem, IDA::jacdense); - if (flag != 0) - { - error ("Dense Jacobian not set"); - } - } - } - } - - void - IDA::jacdense_impl (long int Neq, realtype t, realtype cj, - N_Vector& yy, N_Vector& yyp, DlsMat& JJ) - - { - BEGIN_INTERRUPT_WITH_EXCEPTIONS; - - ColumnVector y = - NVecToCol (yy, Neq); - - ColumnVector yp = - NVecToCol (yyp, Neq); - - Matrix jac; - - if (havejacfun) - jac = (*jacfun) (y, yp, t, cj, ida_jac); - else - jac = (*jacdcell) (dfdy, dfdyp, cj); - - std::copy (jac.fortran_vec (), - jac.fortran_vec () + jac.numel (), - JJ -> data); - - END_INTERRUPT_WITH_EXCEPTIONS; - } - - void - IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, - SlsMat& Jac) - - { - BEGIN_INTERRUPT_WITH_EXCEPTIONS; - - ColumnVector y = - NVecToCol (yy, num); - - ColumnVector yp = - NVecToCol (yyp, num); - - SparseMatrix jac; - - if (havejacfun) - jac = (*jacspfun) (y, yp, t, cj, ida_jac); - else - jac = (*jacspcell) (spdfdy, spdfdyp, cj); - - SparseSetMatToZero (Jac); - int *colptrs = *(Jac -> colptrs); - int *rowvals = *(Jac -> rowvals); - - for (int i = 0; i < num + 1; i++) - colptrs[i] = jac.cidx(i); - - for (int i = 0; i < jac.nnz (); i++) - { - rowvals[i] = jac.ridx (i); - Jac -> data[i] = jac.data (i); - } - - END_INTERRUPT_WITH_EXCEPTIONS; - } - - ColumnVector - IDA::NVecToCol (N_Vector& v, long int n) - { - ColumnVector data (n); - realtype *punt; - punt = NV_DATA_S (v); - - for (octave_idx_type i = 0; i < n; i++) - data (i) = punt [i]; - - return data; - } - - N_Vector - IDA::ColToNVec (const ColumnVector &data, long int n) - { - N_Vector v = - N_VNew_Serial (n); - - realtype * punt; - punt = NV_DATA_S (v); - - for (octave_idx_type i = 0; i < n; i++) - punt [i] = data (i); - - return v; - } - - void - IDA::set_userdata (void) - { - void * userdata = this; - - int flag = - IDASetUserData(mem, userdata); - if (flag != 0) - { - error ("User data not set"); - } - } - - void - IDA::initialize (void) - { - num = y0.numel(); - mem = IDACreate (); - - N_Vector yy = - ColToNVec(y0, num); - - N_Vector yyp = - ColToNVec(yp0, num); - - IDA::set_userdata (); - - int flag = - IDAInit (mem, IDA::resfun, t0, yy, yyp); - if (flag != 0) - { - error ("IDA not initialized"); - } - } - - void - IDA::set_tolerance (ColumnVector& abstol, realtype reltol) - { - N_Vector abs_tol = - ColToNVec (abstol, num); - - int flag = - IDASVtolerances (mem, reltol, abs_tol); - if (flag != 0) - { - error ("IDA: Tolerance not set"); - } - N_VDestroy_Serial (abs_tol); - } - - void - IDA::set_tolerance (realtype abstol, realtype reltol) - { - int flag = - IDASStolerances (mem, reltol, abstol); - if (flag != 0) - { - error ("IDA: Tolerance not set"); - } - } - - octave_value_list - IDA::integrate (const int numt, const ColumnVector& tspan, - const ColumnVector& y, const ColumnVector& yp, - const int refine, bool haverefine, bool haveoutputfcn, - octave_function *output_fcn, bool haveoutputsel, - ColumnVector& outputsel, bool haveeventfunction, - octave_function *event_fcn) - { - // Set up output - ColumnVector tout, yout (num), ypout (num), ysel (outputsel.numel ()); - ColumnVector ie, te, oldval, oldisterminal, olddir; - Matrix output, ye; - int cont = 0, flag = 0, temp = 0; - bool status = 0; - std::string string = ""; - ColumnVector yold = y; - octave_value_list retval; - - realtype tsol = - tspan (0); - realtype tend = - tspan (numt - 1); - - N_Vector yyp = - ColToNVec (yp, num); - - N_Vector yy = - ColToNVec (y, num); - - // Initialize OutputFcn - if (haveoutputfcn) - status = IDA::outputfun (output_fcn, haveoutputsel, y, - tsol, tend, outputsel, "init"); - - // Initialize Events - if (haveeventfunction) - status = IDA::event (event_fcn, te, ye, ie, tsol, y, - "init", yp, oldval, oldisterminal, - olddir, cont, temp, tsol, yold); - - if (numt > 2) - { - // First output value - tout.resize (numt); - tout (0) = tsol; - output.resize (numt, num); - - for (octave_idx_type i = 0; i < num; i++) - output.elem (0, i) = y.elem (i); - - //Main loop - for (octave_idx_type j = 1; j < numt && status == 0; j++) - { - // IDANORMAL already interpolates tspan(j) - flag = IDASolve (mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL); - if (flag != 0) - { - error ("IDASolve failed"); - } - - yout = NVecToCol (yy, num); - ypout = NVecToCol (yyp, num); - tout (j) = tsol; - - for (octave_idx_type i = 0; i < num; i++) - output.elem (j, i) = yout.elem (i); - - if (haveoutputfcn) - status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol, - tend, outputsel, string); - - if (haveeventfunction) - status = IDA::event (event_fcn, te, ye, ie, tout (j), yout, - string, ypout, oldval, oldisterminal, - olddir, j, temp, tout (j - 1), yold); - - // If integration is stopped, return only the reached steps - if (status == 1) - { - output.resize (j + 1, num); - tout.resize (j + 1); - } - - } - } - else // numel (tspan) == 2 - { - // First output value - tout.resize (1); - tout (0) = tsol; - output.resize (1, num); - - for (octave_idx_type i = 0; i < num; i++) - output.elem (0, i) = y.elem (i); - - bool posdirection = - (tend > tsol); - - //main loop - while (((posdirection == 1 && tsol < tend) - || (posdirection == 0 && tsol > tend)) - && status == 0) - { - flag = IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP); - if (flag != 0) - { - error ("IDASolve failed"); - } - - if (haverefine) - status = IDA::interpolate (cont, output, tout, refine, tend, - haveoutputfcn, haveoutputsel, - output_fcn, outputsel, - haveeventfunction, event_fcn, te, - ye, ie, oldval, oldisterminal, - olddir, temp, yold); - - ypout = NVecToCol (yyp, num); - cont += 1; - output.resize (cont + 1, num); // This may be not efficient - tout.resize (cont + 1); - tout (cont) = tsol; - yout = NVecToCol (yy, num); - - for (octave_idx_type i = 0; i < num; i++) - output.elem (cont, i) = yout.elem (i); - - if (haveoutputfcn && ! haverefine && tout (cont) < tend) - status = IDA::outputfun (output_fcn, haveoutputsel, yout, tsol, - tend, outputsel, string); - - if (haveeventfunction && ! haverefine && tout (cont) < tend) - status = IDA::event (event_fcn, te, ye, ie, tout (cont), yout, - string, ypout, oldval, oldisterminal, - olddir, cont, temp, tout (cont - 1), yold); - } - if (status == 0) - { - // Interpolate in tend - N_Vector dky = - N_VNew_Serial (num); - - flag = IDAGetDky (mem, tend, 0, dky); - if (flag != 0) - { - error ("IDA failed to interpolate y"); - } - - tout (cont) = tend; - yout = NVecToCol (dky, num); - - for (octave_idx_type i = 0; i < num; i++) - output.elem (cont, i) = yout.elem (i); - - // Plot final value - if (haveoutputfcn) - { - status = IDA::outputfun (output_fcn, haveoutputsel, yout, - tend, tend, outputsel, string); - - // Events during last step - if (haveeventfunction) - status = IDA::event (event_fcn, te, ye, ie, tend, yout, - string, ypout, oldval, oldisterminal, - olddir, cont, temp, tout (cont - 1), - yold); - } - - N_VDestroy_Serial (dky); - } - // Cleanup plotter - status = IDA::outputfun (output_fcn, haveoutputsel, yout, tend, tend, - outputsel, "done"); - - } - - retval (0) = tout; - retval (1) = output; - retval (2) = te; - retval (3) = ye; - retval (4) = ie; - - return retval; - } - - bool - IDA::event (octave_function *event_fcn, - ColumnVector& te, Matrix& ye, ColumnVector& ie, - realtype tsol, const ColumnVector& y, std::string flag, - const ColumnVector& yp, ColumnVector& oldval, - ColumnVector& oldisterminal, ColumnVector& olddir, int cont, - int& temp, realtype told, ColumnVector& yold) - { - bool status = 0; - octave_value_list args (3); - octave_value_list output (3); - ColumnVector val, isterminal, dir; - args (0) = tsol; - args (1) = y; - args (2) = yp; - - // cont is the number of steps reached by the solver - // temp is the number of events registered - - if (flag == "init") - { - output = feval (event_fcn, args, 3); - oldval = output(0).vector_value (); - oldisterminal = output(1).vector_value (); - olddir = output(2).vector_value (); - } - else if (flag == "") - { - ColumnVector index (0); - output = feval (event_fcn, args, 3); - val = output(0).vector_value (); - isterminal = output(1).vector_value (); - dir = output(2).vector_value (); - - // Get the index of the changed values - for (octave_idx_type i = 0; i < val.numel (); i++) - { - if ((val(i) > 0 && oldval(i) < 0 && dir(i) != -1) // increasing - || (val(i) < 0 && oldval(i) > 0 && dir(i) != 1)) // decreasing - { - index.resize (index.numel () + 1); - index (index.numel () - 1) = i; - } - } - - if (cont == 1 && index.numel () > 0) // Events in first step - { - temp = 1; // register only the first event - te.resize (1); - ye.resize (1, num); - ie.resize (1); - - // Linear interpolation - ie (0) = index (0); - te (0) = tsol - val (index (0)) * (tsol - told) - / (val (index (0)) - oldval (index (0))); - - ColumnVector ytemp = - y - ((tsol - te (0)) * (y - yold ) / (tsol - told)); - - for (octave_idx_type i = 0; i < num; i++) - ye.elem (0, i) = ytemp.elem (i); - - } - else if (index.numel () > 0) - // Not first step: register all events and test if stop integration or not - { - te.resize (temp + index.numel ()); - ye.resize (temp + index.numel (), num); - ie.resize (temp + index.numel ()); - - for (octave_idx_type i = 0; i < index.numel (); i++) - { - - if (isterminal (index (i)) == 1) - status = 1; // Stop integration - - // Linear interpolation - ie (temp + i) = index (i); - te (temp + i) = tsol - - val (index (i)) * (tsol - told) / - (val (index (i)) - oldval (index (i))); - - ColumnVector ytemp = - y - (tsol - te (temp + i)) * (y - yold) / (tsol - told); - - for (octave_idx_type j = 0; j < num; j++) - ye.elem (temp + i, j) = ytemp.elem (j); - - } - - temp += index.numel (); - } - - // Update variables - yold = y; - told = tsol; - olddir = dir; - oldval = val; - oldisterminal = isterminal; - } - - return status; - } - - bool - IDA::interpolate (int& cont, Matrix& output, ColumnVector& tout, - int refine, realtype tend, bool haveoutputfcn, - bool haveoutputsel, octave_function *output_fcn, - ColumnVector& outputsel, bool haveeventfunction, - octave_function * event_fcn, ColumnVector& te, - Matrix& ye, ColumnVector& ie, ColumnVector& oldval, - ColumnVector& oldisterminal, ColumnVector& olddir, - int& temp, ColumnVector& yold) - { - realtype h = 0, tcur = 0; - bool status = 0; - - N_Vector dky = - N_VNew_Serial (num); - - N_Vector dkyp = - N_VNew_Serial (num); - - ColumnVector yout (num); - ColumnVector ypout (num); - std::string string = ""; - - int flag = - IDAGetLastStep (mem, &h); - - if (flag != 0) - { - error ("IDA failed to return last step"); - } - flag = IDAGetCurrentTime (mem, &tcur); - if (flag != 0) - { - error ("IDA failed to return the current time"); - } - - realtype tin = - tcur - h; - - realtype step = - h / refine; - - for (octave_idx_type i = 1; i < refine && tin + step * i < tend - && status == 0; i++) - { - flag = IDAGetDky (mem, tin + step*i, 0, dky); - if (flag != 0) - { - error ("IDA failed to interpolate y"); - } - flag = IDAGetDky (mem, tin + step*i, 1, dkyp); - if (flag != 0) - { - error ("IDA failed to interpolate yp"); - } - cont += 1; - output.resize (cont + 1, num); - tout.resize (cont + 1); - - tout (cont) = tin + step * i; - yout = NVecToCol (dky, num); - ypout = NVecToCol (dkyp, num); - - for (octave_idx_type j = 0; j < num; j++) - output.elem (cont, j) = yout.elem (j); - - if (haveoutputfcn) - status = IDA::outputfun (output_fcn, haveoutputsel, yout, - tout (cont), tend, outputsel, ""); - - if (haveeventfunction) - status = IDA::event (event_fcn, te, ye, ie, tout (cont), - yout, string, ypout, oldval, - oldisterminal, olddir, cont, temp, - tout (cont - 1), yold); - } - N_VDestroy_Serial (dky); - return status; - } - - bool - IDA::outputfun (octave_function *output_fcn, bool haveoutputsel, - const ColumnVector& yout, realtype tsol, - realtype tend, ColumnVector& outputsel, - std::string flag) - { - bool status = 0; - octave_value_list output, val; - ColumnVector ysel (outputsel.numel ()); - - if (haveoutputsel) - { - for (octave_idx_type i = 0; i < outputsel.numel (); i++) - ysel (i) = yout (outputsel (i)); - - output (1) = ysel; - } - else - output (1) = yout; - - output (2) = flag; - - if (flag == "init") - { - ColumnVector toutput (2); - toutput (0) = tsol; - toutput (1) = tend; - output (0) = toutput; - - feval (output_fcn, output, 0); - } - else if (flag == "") - { - output (0) = tsol; - val = feval (output_fcn, output, 1); - status = val(0).bool_value (); - } - else - { // Cleanup plotter - output (0) = tend; - feval (output_fcn, output, 0); - } - - return status; - } - - void - IDA::set_maxstep (realtype maxstep) - { - int flag = - IDASetMaxStep (mem, maxstep); - - if (flag != 0) - { - error ("IDA: Max Step not set"); - } - } - - void - IDA::set_initialstep (realtype initialstep) - { - int flag = - IDASetInitStep (mem, initialstep); - - if (flag != 0) - { - error ("IDA: Initial Step not set"); - } - } - - void - IDA::set_maxorder (int maxorder) - { - int flag = - IDASetMaxOrd (mem, maxorder); - - if (flag != 0) - { - error ("IDA: Max Order not set"); - } - } - - void - IDA::print_stat (void) - { - long int nsteps = 0, netfails = 0, nrevals = 0; - int flag = - IDAGetNumSteps(mem, &nsteps); - - if (flag != 0) - { - error ("IDA failed to return the number of internal steps"); - } - flag = IDAGetNumErrTestFails(mem, &netfails); - if (flag != 0) - { - error ("IDA failed to return the number of internal errors"); - } - flag = IDAGetNumResEvals(mem, &nrevals); - if (flag != 0) - { - error ("IDA failed to return the number of residual evaluations"); - } - - std::cout< do_multi_index_op (1, args); - } - catch (octave::execution_exception& e) - { - err_user_supplied_eval (e, "__ode15__"); - } - - retval = tmp(0).vector_value (); - - return retval; - } - - Matrix - ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot, - double t, double cj, octave_function *ida_jc) - { - Matrix retval; - octave_value_list newargs (3); - - newargs (0) = t; - newargs (1) = x; - newargs (2) = xdot; - - octave_value_list tmp; - - try - { - tmp = ida_jc -> do_multi_index_op (2, newargs); - } - catch (octave::execution_exception& e) - { - err_user_supplied_eval (e, "__ode15__"); - } - - retval = tmp(0).matrix_value () + cj * tmp(1).matrix_value (); - - return retval; - } - - SparseMatrix - ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot, - double t, double cj, octave_function *ida_jc) - { - SparseMatrix retval; - octave_value_list newargs (3); - - newargs (0) = t; - newargs (1) = x; - newargs (2) = xdot; - - octave_value_list tmp; - - try - { - tmp = ida_jc -> do_multi_index_op (2, newargs); - } - catch (octave::execution_exception& e) - { - err_user_supplied_eval (e, "__ode15__"); - } - - retval = tmp(0).sparse_matrix_value () + - cj * tmp(1).sparse_matrix_value (); - - return retval; - } - - Matrix - ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj) - { - Matrix retval; - retval = (*dfdy) + cj * (*dfdyp); - return retval; - } - - SparseMatrix - ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp, - double cj) - { - SparseMatrix retval; - retval = (*spdfdy) + cj * (*spdfdyp); - return retval; - } - - octave_value_list - do_ode15 (octave_function *ida_fcn, - const ColumnVector &tspan, - const int numt, - const realtype t0, - const ColumnVector &y0, - const ColumnVector &yp0, - const octave_scalar_map &options) - { - octave_value_list retval; - - // Create object - IDA dae (t0, y0, yp0, ida_fcn, ida_user_function); - - // Set Jacobian - bool havejac = - options.getfield ("havejac").bool_value (); - - bool havejacsparse = - options.getfield ("havejacsparse").bool_value (); - - bool havejacfun = - options.getfield ("havejacfun").bool_value (); - - Matrix ida_dfdy, ida_dfdyp, *dfdy, *dfdyp; - SparseMatrix ida_spdfdy, ida_spdfdyp, *spdfdy, *spdfdyp; - octave_function *ida_jac; - Cell jaccell; - - if (havejac) - { - if (havejacfun) - { - ida_jac = options.getfield ("Jacobian").function_value (); - if (havejacsparse) - dae.set_jacobian (ida_jac, ida_sparse_jac); - else - dae.set_jacobian (ida_jac, ida_dense_jac); - } - else - { - jaccell = options.getfield ("Jacobian").cell_value (); - - if (havejacsparse) - { - ida_spdfdy = jaccell(0).sparse_matrix_value (); - ida_spdfdyp = jaccell(1).sparse_matrix_value (); - spdfdy = &ida_spdfdy; - spdfdyp = &ida_spdfdyp; - dae.set_jacobian (spdfdy, spdfdyp, ida_sparse_cell_jac); - } - else - { - ida_dfdy = jaccell(0).matrix_value (); - ida_dfdyp = jaccell(1).matrix_value (); - dfdy = &ida_dfdy; - dfdyp = &ida_dfdyp; - dae.set_jacobian (dfdy, dfdyp, ida_dense_cell_jac); - } - } - } - - // Initialize IDA - dae.initialize (); - - // Set tolerances - realtype rel_tol = - options.getfield("RelTol").double_value (); - - bool haveabstolvec = - options.getfield ("haveabstolvec").bool_value (); - - if (haveabstolvec) - { - ColumnVector abs_tol = - options.getfield("AbsTol").vector_value (); - - dae.set_tolerance (abs_tol, rel_tol); - } - else - { - realtype abs_tol = - options.getfield("AbsTol").double_value (); - - dae.set_tolerance (abs_tol, rel_tol); - } - - //Set max step - realtype maxstep = - options.getfield("MaxStep").double_value (); - - dae.set_maxstep (maxstep); - - //Set initial step - if (!(options.getfield("InitialStep").is_empty ())) - { - realtype initialstep = - options.getfield("InitialStep").double_value (); - - dae.set_initialstep (initialstep); - } - - //Set max order FIXME: it doesn't work - int maxorder = - options.getfield("MaxOrder").int_value (); - - dae.set_maxorder (maxorder); - - //Set Refine - const int refine = - options.getfield("Refine").int_value (); - - bool haverefine = - (refine > 1); - - octave_function *output_fcn = NULL; - ColumnVector outputsel; - - // OutputFcn - bool haveoutputfunction = - options.getfield("haveoutputfunction").bool_value (); - - if (haveoutputfunction) - output_fcn = options.getfield("OutputFcn").function_value (); - - // OutputSel - bool haveoutputsel = - options.getfield("haveoutputselection").bool_value (); - - if (haveoutputsel) - outputsel = options.getfield("OutputSel").vector_value (); - - octave_function *event_fcn = NULL; - - // Events - bool haveeventfunction = - options.getfield("haveeventfunction").bool_value (); - - if (haveeventfunction) - event_fcn = options.getfield("Events").function_value (); - - // Set up linear solver - dae.set_up (); - - // Integrate - retval = dae.integrate (numt, tspan, y0, yp0, refine, - haverefine, haveoutputfunction, - output_fcn, haveoutputsel, outputsel, - haveeventfunction, event_fcn); - - // Statistics - bool havestats = - options.getfield("havestats").bool_value (); - - if (havestats) - dae.print_stat (); - - return retval; - } -} -#endif - - -DEFUN_DLD (__ode15__, args, nargout, doc: /* -*- texinfo -*- -@deftypefn {} {@var{t}, @var{y} =} __ode15__ (@var{fun}, @ -@var{tspan}, @var{y0}, @var{yp0}, @var{options}) -Undocumented internal function. -@end deftypefn */) -{ - -#if defined (HAVE_SUNDIALS) - - // Check number of parameters - int nargin = args.length (); - - if (nargin != 5 || nargout != 5) - print_usage (); - - // Check odefun - octave_function *ida_fcn = NULL; - - octave_value f_arg = args(0); - - if (f_arg.is_function_handle ()) - ida_fcn = f_arg.function_value (); - else - error ("__ode15__: odefun must be a function handle"); - - // Check input tspan - ColumnVector tspan = - args(1).xvector_value ("__ode15__: TRANGE must be a vector of numbers"); - - int numt = - tspan.numel (); - - realtype t0 = - tspan (0); - - if (numt < 2) - error ("__ode15__: TRANGE must contain at least 2 elements"); - else if (!(tspan.is_sorted ()) || (tspan(0) == tspan(numt - 1))) - error ("__ode15__: TRANGE must be strictly monotonic"); - - // input y0 and yp0 - ColumnVector y0 = - args(2).xvector_value ("__ode15__: initial state y0 must be a vector"); - - ColumnVector yp0 = - args(3).xvector_value ("__ode15__: initial state yp0 must be a vector"); - - - if (y0.numel () != yp0.numel ()) - error ("__ode15__: initial state y0 and yp0 must have the same length"); - else if (y0.numel () < 1) - error ("__ode15__: initial state yp0 must be a vector or a scalar"); - - - if (! args(4).is_map ()) - error ("__ode15__: OPTS argument must be a structure"); - - octave_scalar_map options = - args(4).xscalar_map_value ("__ode15__:", - "OPTS argument must be a scalar structure"); - - - return octave::do_ode15 (ida_fcn, tspan, numt, t0, - y0, yp0, options); - - -#else - - octave_unused_parameter (args); - err_disabled_feature ("__ode15__", "sundials_ida, sundials_nvecserial"); - -#endif -} - - - diff --git a/libinterp/dldfcn/module-files b/libinterp/dldfcn/module-files --- a/libinterp/dldfcn/module-files +++ b/libinterp/dldfcn/module-files @@ -5,7 +5,6 @@ __glpk__.cc|$(GLPK_CPPFLAGS)|$(GLPK_LDFLAGS)|$(GLPK_LIBS) __init_fltk__.cc|$(FLTK_CPPFLAGS) $(FT2_CPPFLAGS) $(FONTCONFIG_CPPFLAGS)|$(FLTK_LDFLAGS) $(FT2_LDFLAGS)|$(FLTK_LIBS) $(FT2_LIBS) $(OPENGL_LIBS) __init_gnuplot__.cc|$(FT2_CPPFLAGS) $(FONTCONFIG_CPPFLAGS)|| -__ode15__.cc|$(SUNDIALS_XCPPFLAGS)|$(SUNDIALS_XLDFLAGS)|$(SUNDIALS_XLIBS) __osmesa_print__.cc|$(OSMESA_CPPFLAGS) $(FT2_CPPFLAGS)|$(OSMESA_LDFLAGS)|$(OSMESA_LIBS) __voronoi__.cc|$(QHULL_CPPFLAGS)|$(QHULL_LDFLAGS)|$(QHULL_LIBS) amd.cc|$(SPARSE_XCPPFLAGS)|$(SPARSE_XLDFLAGS)|$(SPARSE_XLIBS) diff --git a/m4/acinclude.m4 b/m4/acinclude.m4 --- a/m4/acinclude.m4 +++ b/m4/acinclude.m4 @@ -1697,32 +1697,6 @@ ]) ]) dnl -dnl Check whether sundials_ida library is configured with double precision realtype -dnl -AC_DEFUN([OCTAVE_CHECK_SIZEOF_SUNDIALS_IDA_REALTYPE], [ - AC_CACHE_CHECK([whether sundials_ida is configured with double precision realtype], - [octave_cv_sizeof_ida_double], - [AC_RUN_IFELSE([AC_LANG_PROGRAM([[ - #if defined (HAVE_IDA_IDA_H) - #include - #else - #include - #endif - #include - ]], [[ - assert (sizeof (double) == sizeof (realtype)); - ]])], - octave_cv_sizeof_ida_double=yes, - octave_cv_sizeof_ida_double=no) - ]) - if test $octave_cv_sizeof_ida_double = yes; then - AC_DEFINE(HAVE_SUNDIALS_IDA_DOUBLE, 1, - [Define to 1 if sundials_ida is configured with realtype as double.]) - else - warn_sundials_ida="sundials_ida library found, but it's not configured with double precision realtype; function ode15i and ode15s will be disabled" - fi -]) -dnl dnl Add warning to final summary. dnl AC_DEFUN([OCTAVE_CONFIGURE_WARNING], [ diff --git a/scripts/ode/decic.m b/scripts/ode/decic.m deleted file mode 100644 --- a/scripts/ode/decic.m +++ /dev/null @@ -1,223 +0,0 @@ -## Copyright (C) 2016, Francesco Faccio -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -## -*- texinfo -*- -## @deftypefn {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fun}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0}) -## @deftypefnx {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fun}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0}, @var{options}) -## @deftypefnx {} {[@var{y0_new}, @var{yp0_new}, @var{resnorm}] =} decic (@dots{}) -## -## Compute consistent initial conditions @var{y0_new} and @var{yp0_new}, given -## initial guesses @var{y0} and @var{yp0}. Choose a maximum of -## @code{length(@var{y0})} components between @var{fixed_y0} and -## @var{fixed_yp0} as fixed values. -## -## @var{fun} is a function handle. The function must accept three inputs where -## the first is time @var{t}, the second is a column vector of unknowns @var{y} -## and the third is a column vector of unknowns @var{yp}. -## -## @var{t0} is the initial time such that @code{@var{fun}(@var{t0}, -## @var{y0_new}, @var{yp0_new}) = 0}, specified as a scalar. -## -## @var{y0} is a vector used as initial guess for @var{y}. -## -## @var{fixed_y0} is a vector which specifies the components of @var{y0} to be -## hold fixed. Choose a maximum of @code{length(@var{y0})} components between -## @var{fixed_y0} and @var{fixed_yp0} as fixed values. -## Set @var{fixed_y0}(i) component to 1 if you want to fix the value of -## @var{y0}(i). -## Set @var{fixed_y0}(i) component to 0 if you want to allow the value of -## @var{y0}(i) to change. -## -## @var{yp0} is a vector used as initial guess for @var{yp}. -## -## @var{fixed_yp0} is a vector which specifies the components of @var{yp0} to -## be hold fixed. Choose a maximum of @code{length(@var{yp0})} components -## between @var{fixed_y0} and @var{fixed_yp0} as fixed values. -## Set @var{fixed_yp0}(i) component to 1 if you want to fix the value of -## @var{yp0}(i). -## Set @var{fixed_yp0}(i) component to 0 if you want to allow the value of -## @var{yp0}(i) to change. -## -## The optional seventh argument @var{options} is a structure array. -## Use @code{odeset} to generate this structure. The relevant options are -## @code{RelTol} and @code{AbsTol} which specify the error thresholds used to -## compute the initial conditions. -## -## The function typically returns two outputs. Variable @var{y0_new} is a -## column vector and contains the consistent initial value of y. The -## output @var{yp0_new} is a column vector and contains the consistent initial -## value of yp. -## -## The optional third output @var{resnorm} is the vector of norm of the -## residuals. If @var{resnorm} is small, @code{decic} has successfully computed -## the initial conditions. If the value of @var{resnorm} is large, use -## @code{RelTol} and @code{AbsTol} to adjust it. -## -## Example: Compute initial conditions of @nospell{Robetson}'s equations: -## -## @example -## @group -## function res = robertsidae(@var{t}, @var{y}, @var{yp}) -## res = [-(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3)); -## -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) + -## 3e7*@var{y}(2)^2); -## @var{y}(1) + @var{y}(2) + @var{y}(3) - 1]; -## endfunction -## @end group -## [@var{y0_new},@var{yp0_new}] = decic (@@robertsidae, 0, [1; 0; 0], [1; 1; 0], -## [-1e-4; 1; 0], [0; 0; 0]); -## @end example -## @seealso{odeset} -## @end deftypefn - -function [y0_new, yp0_new, resnrm] = decic (odefun, t0, y0, fixed_y0, yp0, - fixed_yp0, options) - - if (nargin < 6 || nargin > 7 || nargout > 3) - print_usage (); - endif - - #Check input - if (! isa (odefun, "function_handle")) - error ("Octave:invalid-input-arg", - "decic: FUN must be a valid function handle"); - endif - - if (! isnumeric (t0) || numel (t0) != 1) - error ("Octave:invalid-input-arg", - "decic: INIT must be a numeric scalar value"); - endif - - if (! isnumeric (y0) || ! isvector (y0) || ! isnumeric (fixed_y0) || - ! isvector (fixed_y0) || ! isnumeric (yp0) || ! isvector (yp0)|| - ! isnumeric (fixed_yp0) || ! isvector (fixed_yp0)) - error ("Octave:invalid-input-arg", - "decic: y0, fixed_y0, yp0 and fixed_yp0 must be numeric vectors"); - - elseif (! isequal (numel (y0), numel (fixed_y0), numel (yp0), - numel (fixed_yp0))) - error ("Octave:invalid-input-arg", - "decic: length of y0, fixed_y0, yp0 and fixed_yp0 must be equal"); - endif - - for i = 1:numel (y0) - if (! (fixed_y0(i) == 0 || fixed_y0(i) == 1) || ! (fixed_yp0(i) == 0 - || fixed_yp0(i) == 1)) - error ("Octave:invalid-input-arg", - "decic: fixed_y0 and fixed_yp0 must be boolean vectors"); - endif - endfor - - n = numel (y0); - nl = sum (~fixed_y0); - nu = sum (~fixed_yp0); - - if (n - nl - nu > 0) - error ("Octave:invalid-input-arg", - "decic: you cannot fix more than length(y0) components"); - endif - - #Set default value - TolFun = 0; - TolX = eps; - - #Check AbsTol and RelTol - if (nargin == 7) - if (! isempty (options.AbsTol)) - if (! isscalar (options.AbsTol)) - error ("Octave:invalid-input-arg", - "decic: AbsTol must be a scalar value"); - else - TolFun = options.AbsTol; - endif - endif - - if (! isempty (options.RelTol)) - if (! isscalar (options.RelTol)) - error ("Octave:invalid-input-arg", - "decic: RelTol must be a scalar value"); - else - TolX = options.RelTol; - endif - endif - endif - - x0 = [y0(~fixed_y0); yp0(~fixed_yp0)]; - opt = optimset ('tolfun', TolFun, 'tolx', TolX, 'display', 'iter-detailed'); - x = fminunc (@(x) objective (x, t0, y0, fixed_y0, yp0, fixed_yp0, nl, nu, - odefun), - x0, opt); - - y0_new = y0; - yp0_new = yp0; - - y0_new(~fixed_y0) = x(1:nl); - yp0_new(~fixed_yp0) = x(nl+1:nl+nu); - resnrm = odefun (t0, y0_new, yp0_new); - -endfunction - -function res = objective (x, t0, y0, fixed_y0, yp0, - fixed_yp0, nl, nu, odefun) - y = y0; - y(~fixed_y0) = x(1:nl); - yp = yp0; - yp(~fixed_yp0) = x(nl+1:nl+nu); - res = sqrt (sum (odefun (t0, y, yp) .^ 2)); -endfunction - -%!function res = rob (t, y, yp) -%! res =[-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); -(yp(2) - 0.04*y(1) + -%! + 1e4*y(2)*y(3) + 3e7*y(2)^2); y(1) + y(2) + y(3) - 1]; -%!endfunction - -%!test # Without options -%! ref1 = [1;0;0]; -%! ref2 = [-4e-2; 4e-2; 0]; -%! [ynew,ypnew] = decic (@rob,0,[1;0;0],[1;1;0],[23;110;0],[0;0;1]); -%! assert ([ynew(1:end), ypnew(1:end)], [ref1(1:end), ref2(1:end)], 1e-10); -%!test # With options -%! ref1 = [1;0;0]; -%! ref2 = [-4e-2; 4e-2; 0]; -%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-4); -%! [ynew,ypnew] = decic (@rob,0,[1;0;0],[1;1;0],[23;110;0],[0;0;1],opt); -%! assert ([ynew(1:end), ypnew(1:end)], [ref1(1:end), ref2(1:end)], 1e-5); - -## Test input validation -%!error decic () -%!error decic (1) -%!error decic (1,2) -%!error decic (1,2,3,4,5,6,7,8) -%!error -%! decic (1, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]); -%!error -%! decic (@rob, [1, 1], [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]); -%!error -%! decic (@rob, 0, [0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]); -%!error -%! decic (@rob, 0, [1; 0; 0], [1; 0],"", [0; 0; 0]); -%!error -%! decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; "1"]); -%!error -%! decic (@rob, 0, [1; 0; 0], [5; 5; 0], [-1e-4; 1; 0], [0; 0; 0]); -%!error -%! decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 4; 0]); -%!error #too many components fixed -%! decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 1; 1]); - - diff --git a/scripts/ode/module.mk b/scripts/ode/module.mk --- a/scripts/ode/module.mk +++ b/scripts/ode/module.mk @@ -4,7 +4,6 @@ scripts_ode_PRIVATE_FCN_FILES = \ scripts/ode/private/AbsRel_norm.m \ - scripts/ode/private/check_default_input.m \ scripts/ode/private/integrate_adaptive.m \ scripts/ode/private/kahan.m \ scripts/ode/private/odedefaults.m \ @@ -16,14 +15,11 @@ scripts/ode/private/starting_stepsize.m scripts_ode_FCN_FILES = \ - scripts/ode/ode15i.m \ - scripts/ode/ode15s.m \ scripts/ode/ode23.m \ scripts/ode/ode45.m \ scripts/ode/odeset.m \ scripts/ode/odeget.m \ - scripts/ode/odeplot.m \ - scripts/ode/decic.m + scripts/ode/odeplot.m scripts_odedir = $(fcnfiledir)/ode diff --git a/scripts/ode/ode15i.m b/scripts/ode/ode15i.m deleted file mode 100644 --- a/scripts/ode/ode15i.m +++ /dev/null @@ -1,540 +0,0 @@ -## Copyright (C) 2016, Francesco Faccio -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -## -*- texinfo -*- -## @deftypefn {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}) -## @deftypefnx {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}, @var{ode_opt}) -## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15i (@dots{}) -## @deftypefnx {} {@var{solution} =} ode15i (@dots{}) -## -## Solve a set of full-implicit Ordinary Differential Equations and -## Differential Algebraic Equations (DAEs) of index 1, with the variable-step, -## variable order BDF (Backward Differentiation Formula) method, which -## ranges from order 1 to 5. -## -## @var{fun} is a function handle, inline function, or string containing the -## name of the function that defines the ODE: @code{f(@var{t},@var{y},@var{yp})}. -## The function must accept three inputs where the first is time @var{t}, the -## second is a column vector of unknowns @var{y} and the third is a column -## vector of unknowns @var{yp}. -## -## @var{trange} specifies the time interval over which the ODE will be -## evaluated. Typically, it is a two-element vector specifying the initial and -## final times (@code{[tinit, tfinal]}). If there are more than two elements -## then the solution will also be evaluated at these intermediate time. -## -## @var{y0} and @var{yp0} contain the initial values for the unknowns @var{y} -## and @var{yp}. If they are row vectors then the solution @var{y} will be a -## matrix in which each column is the solution for the corresponding initial -## value in @var{y0} and @var{yp0}. -## -## @var{y0} and @var{yp0} must be consistent initial conditions, meaning that -## @code{f(@var{t},@var{y0},@var{yp0})=0} is satisfied. You can use function -## decic to compute consistent initial conditions, given initial guesses. -## -## The optional fifth argument @var{ode_opt} specifies non-default options to -## the ODE solver. It is a structure generated by @code{odeset}. -## -## The function typically returns two outputs. Variable @var{t} is a -## column vector and contains the times where the solution was found. The -## output @var{y} is a matrix in which each column refers to a different -## unknown of the problem and each row corresponds to a time in @var{t}. -## -## The output can also be returned as a structure @var{solution} which -## has field @var{x} containing the time where the solution was evaluated and -## field @var{y} containing the solution matrix for the times in @var{x}. -## Use @code{fieldnames (@var{solution})} to see the other fields and -## additional information returned. -## -## If using the @qcode{"Events"} option then three additional outputs may -## be returned. @var{te} holds the time when an Event function returned a -## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} -## contains an index indicating which Event function was triggered in the case -## of multiple Event functions. -## -## This function can be called with two output arguments: @var{t} and @var{y}. -## Variable @var{t} is a column vector and contains the time stamps, instead -## @var{y} is a matrix in which each column refers to a different unknown of -## the problem and the rows number is the same of @var{t} rows number so -## that each row of @var{y} contains the values of all unknowns at the time -## value contained in the corresponding row in @var{t}. -## -## Example: Solve the @nospell{Robetson}'s equations: -## -## @example -## @group -## function res = robertsidae(@var{t}, @var{y}, @var{yp}) -## res = [-(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3)); -## -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) + -## 3e7*@var{y}(2)^2); -## @var{y}(1) + @var{y}(2) + @var{y}(3) - 1]; -## endfunction -## [@var{t},@var{y}] = ode15i (@@robertsidae, [0 1e3], [1; 0; 0],[-1e-4; 1e-4; 0]); -## @end group -## @end example -## @seealso{decic, odeset, odeget} -## @end deftypefn - -function varargout = ode15i (fun, trange, y0, yp0, varargin) - - solver = 'ode15i'; - - if (nargin < 4) - print_usage (); - endif - - n = numel (y0); - - if (nargin > 4) - options = varargin{1}; - else - options = odeset (); - endif - - ## Check fun, trange, y0, yp0 - fun = check_default_input (fun, trange, solver, y0, yp0); - - if (! isempty (options.Jacobian)) - if (ischar (options.Jacobian)) - try - options.Jacobian = str2func (options.Jacobian); - catch - warning (lasterr); - end_try_catch - if (! isa (options.Jacobian, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - endif - endif - - if (! isempty (options.OutputFcn)) - if (ischar (options.OutputFcn)) - try - options.OutputFcn = str2func (options.OutputFcn); - catch - warning (lasterr); - end_try_catch - if (! isa (options.OutputFcn, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "OutputFcn"); - endif - endif - endif - - if (! isempty (options.Events)) - if (ischar (options.Events)) - try - options.Events = str2func (options.Events); - catch - warning (lasterr); - end_try_catch - if (! isa (options.Events, "function_handle") && ! ismatrix (options.Events)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Events"); - endif - endif - endif - - [defaults, classes, attributes] = ... - odedefaults (n, trange(1), trange(end)); - - persistent ignorefields = {"NonNegative", "Mass", ... - "MStateDependence", "MvPattern", ... - "MassSingular", "InitialSlope", "BDF"}; - - defaults = rmfield (defaults, ignorefields); - classes = rmfield (classes, ignorefields); - attributes = rmfield (attributes, ignorefields); - - classes = odeset (classes, 'Vectorized', {}); - attributes = ... - odeset (attributes, 'Jacobian', {}, 'Vectorized', {}); - - options = ... - odemergeopts ("ode15i", options, defaults, - classes, attributes, solver); - - ## Jacobian - options.havejac = false; - options.havejacsparse = false; - options.havejacfun = false; - - if (! isempty (options.Jacobian)) - options.havejac = true; - if (iscell (options.Jacobian)) - if (numel (options.Jacobian) == 2) - if (issparse (options.Jacobian{1}) && issparse (options.Jacobian{2})) ## Jac is sparse cell - options.havejacsparse = true; - endif - - if (any (size (options.Jacobian{1}) != [n n]) - || any (size (options.Jacobian{2}) != [n n]) - || ! isnumeric (options.Jacobian{1}) - || ! isnumeric (options.Jacobian{2}) - || ! isreal (options.Jacobian{1}) - || ! isreal (options.Jacobian{2})) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - - elseif (isa (options.Jacobian, "function_handle")) - options.havejacfun = true; - if (nargin (options.Jacobian) == 3) - [A, B] = options.Jacobian (trange(1), y0, yp0); - if (issparse (A) && issparse (B)) - options.havejacsparse = true; ## Jac is sparse fun - endif - - if (any (size (A) != [n n]) || any (size (B) != [n n]) - || ! isnumeric (A) || ! isnumeric (B) || ! isreal (A) - || ! isreal (B)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - endif - - ## Abstol and Reltol - - options.haveabstolvec = false; - - if (numel (options.AbsTol) != 1 && numel (options.AbsTol) != n) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "AbsTol"); - - elseif (numel (options.AbsTol) == n) - options.haveabstolvec = true; - endif - - ## Stats - options.havestats = false; - if (strcmp (options.Stats, "on")) - options.havestats = true; - endif - - ## Don't use Refine when the output is a structure - if (nargout == 1) - options.Refine = 1; - endif - - ## OutputFcn and OutputSel - if (isempty (options.OutputFcn) && nargout == 0) - options.OutputFcn = @odeplot; - options.haveoutputfunction = true; - else - options.haveoutputfunction = ! isempty (options.OutputFcn); - endif - - options.haveoutputselection = ! isempty (options.OutputSel); - if (options.haveoutputselection) - options.OutputSel = options.OutputSel - 1; - endif - - ## Events - options.haveeventfunction = ! isempty (options.Events); - - - [t, y, te, ye, ie] = __ode15__ (fun, trange, y0, yp0, options); - - - if (nargout == 2) - varargout{1} = t; - varargout{2} = y; - elseif (nargout == 1) - varargout{1}.x = t; # Time stamps are saved in field x - varargout{1}.y = y; # Results are saved in field y - varargout{1}.solver = solver; - if (options.haveeventfunction) - varargout{1}.xe = te; # Time info when an event occurred - varargout{1}.ye = ye; # Results when an event occurred - varargout{1}.ie = ie; # Index info which event occurred - endif - elseif (nargout == 5) - varargout = cell (1,5); - varargout{1} = t; - varargout{2} = y; - if (options.haveeventfunction) - varargout{3} = te; # Time info when an event occurred - varargout{4} = ye; # Results when an event occurred - varargout{5} = ie; # Index info which event occurred - endif - endif - -endfunction - -%!demo -%! -%! ##Solve Robertson's equations with ode15i -%! fun = @ (t, y, yp) [-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); -%! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2); -%! y(1) + y(2) + y(3) - 1]; -%! -%! opt = odeset ('RelTol',1e-4, 'AbsTol', [1e-8, 1e-14, 1e-6]); -%! y0 = [1; 0; 0]; -%! yp0 = [-1e-4; 1e-4; 0]; -%! tspan = [0 4*logspace(-6, 6)]; -%! -%! [t, y] = ode15i (fun, tspan, y0, yp0, opt); -%! -%! y (:,2) = 1e4 * y (:, 2); -%! figure (2); -%! semilogx (t, y, 'o') -%! xlabel ('time'); -%! ylabel ('species concentration'); -%! title ('Robertson DAE problem with a Conservation Law'); -%! legend ('y1', 'y2', 'y3'); - -%!function res = rob (t, y, yp) -%! res =[-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); -%! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2); -%! y(1) + y(2) + y(3) - 1]; -%!endfunction -%! -%!function ref = fref() -%! ref = [100, 0.617234887614937, 0.000006153591397, 0.382758958793666]; -%!endfunction -%! -%!function ref2 = fref2() -%! ref2 = [4e6 0 0 1]; -%!endfunction -%! -%!function [DFDY, DFDYP] = jacfundense(t, y, yp) -%! DFDY = [-0.04, 1e4*y(3), 1e4*y(2); -%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); -%! 1, 1, 1]; -%! DFDYP = [-1, 0, 0; -%! 0, -1, 0; -%! 0, 0, 0]; -%!endfunction -%! -%!function [DFDY, DFDYP] = jacfunsparse(t, y, yp) -%! DFDY = sparse ([-0.04, 1e4*y(3), 1e4*y(2); -%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); -%! 1, 1, 1]); -%! DFDYP = sparse ([-1, 0, 0; -%! 0, -1, 0; -%! 0, 0, 0]); -%!endfunction -%!function [DFDY, DFDYP] = jacwrong(t, y, yp) -%! DFDY = [-0.04, 1e4*y(3); -%! 0.04, -1e4*y(3)-6e7*y(2)]; -%! DFDYP = [-1, 0; -%! 0, -1]; -%!endfunction -%!function [DFDY, DFDYP, A] = jacwrong2(t, y, yp) -%! DFDY = [-0.04, 1e4*y(3), 1e4*y(2); -%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); -%! 1, 1, 1]; -%! DFDYP = [-1, 0, 0; -%! 0, -1, 0; -%! 0, 0, 0]; -%! A = DFDY; -%!endfunction -%!function [val, isterminal, direction] = ff (t, y, yp) -%! isterminal = [0 1]; -%! if (t < 1e1) -%! val = [-1, -2]; -%! else -%! val = [1 3]; -%! endif -%! -%! direction = [1 0]; -%!endfunction - -%!test # anonymous function instead of real function -%! ref = [0.049787079136413]; -%! ff = @(t, u, udot) udot + 3 * u; -%! [t, y] = ode15i (ff, 0:1, 1, -3); -%! assert ([t(end), y(end)], [1, ref], 1e-3); -%!test # function passed as string -%! [t, y] = ode15i ('rob',[0 100 200], [1;0;0], [-1e-4;1e-4;0]); -%! assert ([t(2), y(2,:)], fref, 1e-3); -%!test # solve in intermidiate step -%! [t, y] = ode15i (@rob,[0 100 200], [1;0;0], [-1e-4;1e-4;0]); -%! assert ([t(2), y(2,:)], fref, 1e-3); -%!test # numel(trange) = 2 final value -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); -%! assert ([t(end), y(end,:)], fref, 1e-5); -%!test # With empty options -%! opt = odeset(); -%! [t, y] = ode15i (@rob,[0 1e6 2e6 3e6 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref2, 1e-3); -%! opt = odeset(); -%!test # Without options -%! [t, y] = ode15i (@rob,[0 1e6 2e6 3e6 4e6], [1;0;0], [-1e-4;1e-4;0]); -%! assert ([t(end), y(end,:)], fref2, 1e-3); -%!test # InitialStep option -%! opt = odeset ("InitialStep", 1e-8); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(2)-t(1)], [1e-8], 1e-9); -%!test # MaxStep option -%! opt = odeset ("MaxStep", 1e-3); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); -%! assert ([t(5)-t(4)], [1e-3], 1e-3); -%!test # AbsTol scalar option -%! opt = odeset ("AbsTol", 1e-8); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # AbsTol scalar and RelTol option -%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-6); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # AbsTol vector option -%! opt = odeset ("AbsTol", [1e-8, 1e-14,1e-6]); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # AbsTol vector and RelTol option -%! opt = odeset ("AbsTol", [1e-8, 1e-14,1e-6], "RelTol", 1e-6); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # RelTol option -%! opt = odeset ("RelTol", 1e-6); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # Jacobian fun dense -%! opt = odeset ("Jacobian", @jacfundense); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # Jacobian fun dense as string -%! opt = odeset ("Jacobian", 'jacfundense'); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # Jacobian fun sparse -%! opt = odeset ("Jacobian", @jacfunsparse, "AbsTol", 1e-7, "RelTol", 1e-7); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), y(end,:)], fref, 1e-3); -%!test # Solve in backward direction starting at t=100 -%! YPref = [-0.001135972751027; -0.000000027483627; 0.001136000234654]; -%! Yref = [0.617234887614937, 0.000006153591397, 0.382758958793666]; -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); -%! [t2, y2] = ode15i (@rob,[100 0], Yref', YPref); -%! assert ([t2(end), y2(end,:)], [0 1 0 0], 2e-2); -#%!test # Solve in backward direction with MaxStep option -%! YPref = [-0.001135972751027; -0.000000027483627; 0.001136000234654]; -%! Yref = [0.617234887614937, 0.000006153591397, 0.382758958793666]; -%! opt = odeset ('MaxStep', 1e-2); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); -%! [t2, y2] = ode15i (@rob,[100 0], Yref', YPref, opt); -%! assert ([t2(end), y2(end,:)], [0 1 0 0], 2e-2); -%! assert ([t2(9)-t2(10)], [1e-2], 1e-2); -#%!test # Solve in backward direction starting with intermidiate step -%! YPref = [-0.001135972751027; -0.000000027483627; 0.001136000234654]; -%! Yref = [0.617234887614937, 0.000006153591397, 0.382758958793666]; -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); -%! [t2, y2] = ode15i (@rob,[100 5 0], Yref', YPref); -%! assert ([t2(end), y2(end,:)], [0 1 0 0], 2e-2); -%!test # Refine -%! opt = odeset ("Refine", 3); -%! [t, y] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0]); -%! [t2, y2] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([numel(t2)], numel(t)*3, 3); -%!test # Refine ignored if numel (trange) > 2 -%! opt = odeset ("Refine", 3); -%! [t, y] = ode15i (@rob,[0 10 100], [1;0;0], [-1e-4;1e-4;0]); -%! [t2, y2] = ode15i (@rob,[0 10 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([numel(t2)], numel(t)); -%!test # Events option add further elements in sol -%! opt = odeset ("Events", @ff); -%! sol = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert (isfield (sol, "ie")); -%! assert (sol.ie, [0;1]); -%! assert (isfield (sol, "xe")); -%! assert (isfield (sol, "ye")); -%! assert (sol.x(end), 10, 1); -%!test # Events option, five output arguments -%! opt = odeset ("Events", @ff); -%! [t, y, te, ye, ie] = ode15i (@rob,[0 100], [1;0;0], [-1e-4;1e-4;0], opt); -%! assert ([t(end), te', ie'], [10, 10, 10, 0, 1], [1, 0.2, 0.2, 0, 0]); - -%!error # Jacobian fun wrong dimension -%! opt = odeset ("Jacobian", @jacwrong); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!error # Jacobian cell dense wrong dimension -%! DFDY = [-0.04, 1; -%! 0.04, 1]; -%! DFDYP = [-1, 0, 0; -%! 0, -1, 0; -%! 0, 0, 0]; -%! opt = odeset ("Jacobian", {DFDY, DFDYP}); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!error # Jacobian cell sparse wrong dimension -%! DFDY = sparse ([-0.04, 1; -%! 0.04, 1]); -%! DFDYP = sparse ([-1, 0, 0; -%! 0, -1, 0; -%! 0, 0, 0]); -%! opt = odeset ("Jacobian", {DFDY, DFDYP}); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!error # Jacobian cell wrong number of matrices -%! A = [1 2 3; 4 5 6; 7 8 9]; -%! opt = odeset ("Jacobian", {A,A,A}); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!error # Jacobian single matrix -%! opt = odeset ("Jacobian", [1 2 3; 4 5 6; 7 8 9]); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!error # Jacobian single matrix wrong dimension -%! opt = odeset ("Jacobian", [1 2 3; 4 5 6]); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!error # Jacobian strange field -%! opt = odeset ("Jacobian", "foo"); -%! [t, y] = ode15i (@rob,[0 4e6], [1;0;0], [-1e-4;1e-4;0], opt); -%!function ydot = fun (t, y, yp) -%! ydot = [y - yp]; -%!endfunction -%!error ode15i (); -%!error ode15i (1); -%!error ode15i (1, 1, 1); -%!error ode15i (1, 1, 1); -%!error ode15i (1, 1, 1, 1); -%!error ode15i (1, 1, 1, 1, 1); -%!error ode15i (1, 1, 1, 1, 1, 1); -%!error ode15i (@fun, 1, 1, 1); -%!error ode15i (@fun, [1, 1], 1, 1); -%!error ode15i (@fun, [1, 2], [1], [1, 2]); -%!error -%! opt = odeset ('RelTol', "foo"); -%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); -%!error -%! opt = odeset ('RelTol', [1, 2]); -%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); -%!error -%! opt = odeset ('RelTol', -2); -%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); -%!error -%! opt = odeset ('AbsTol', "foo"); -%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); -%!error -%! opt = odeset ('AbsTol', -1); -%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); -%!error -%! opt = odeset ('AbsTol', [1, 1, 1]); -%! [t, y] = ode15i (@fun, [0 2], [2], [2], opt); - - diff --git a/scripts/ode/ode15s.m b/scripts/ode/ode15s.m deleted file mode 100644 --- a/scripts/ode/ode15s.m +++ /dev/null @@ -1,685 +0,0 @@ -## Copyright (C) 2016, Francesco Faccio -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -## -*- texinfo -*- -## @deftypefn {} {[@var{t}, @var{y}] =} ode15s (@var{fun}, @var{trange}, @var{y0}) -## @deftypefnx {} {[@var{t}, @var{y}] =} ode15s (@var{fun}, @var{trange}, @var{y0}, @var{ode_opt}) -## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15s (@dots{}) -## @deftypefnx {} {@var{solution} =} ode15s (@dots{}) -## -## Solve a set of stiff Ordinary Differential Equations and stiff semi-explicit -## Differential Algebraic Equations (DAEs) of index 1, with the variable-step, -## variable order BDF (Backward Differentiation Formula) method, which -## ranges from order 1 to 5. -## -## @var{fun} is a function handle, inline function, or string containing the -## name of the function that defines the ODE: @code{f(@var{t},@var{y})}. -## The function must accept two inputs where the first is time @var{t} and the -## second is a column vector of unknowns @var{y}. -## -## @var{trange} specifies the time interval over which the ODE will be -## evaluated. Typically, it is a two-element vector specifying the initial and -## final times (@code{[tinit, tfinal]}). If there are more than two elements -## then the solution will also be evaluated at these intermediate time. -## -## @var{init} contains the initial value for the unknowns. If it is a row -## vector then the solution @var{y} will be a matrix in which each column is -## the solution for the corresponding initial value in @var{init}. -## -## The optional fourth argument @var{ode_opt} specifies non-default options to -## the ODE solver. It is a structure generated by @code{odeset}. -## -## The function typically returns two outputs. Variable @var{t} is a -## column vector and contains the times where the solution was found. The -## output @var{y} is a matrix in which each column refers to a different -## unknown of the problem and each row corresponds to a time in @var{t}. -## -## The output can also be returned as a structure @var{solution} which -## has field @var{x} containing the time where the solution was evaluated and -## field @var{y} containing the solution matrix for the times in @var{x}. -## Use @code{fieldnames (@var{solution})} to see the other fields and -## additional information returned. -## -## If using the @qcode{"Events"} option then three additional outputs may -## be returned. @var{te} holds the time when an Event function returned a -## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} -## contains an index indicating which Event function was triggered in the case -## of multiple Event functions. -## -## This function can be called with two output arguments: @var{t} and @var{y}. -## Variable @var{t} is a column vector and contains the time stamps, instead -## @var{y} is a matrix in which each column refers to a different unknown of -## the problem and the rows number is the same of @var{t} rows number so -## that each row of @var{y} contains the values of all unknowns at the time -## value contained in the corresponding row in @var{t}. -## -## Example: Solve the @nospell{Robetson}'s equations: -## @example -## @group -## function res = robertsidae(@var{t}, @var{y}) -## res = [-0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3); -## 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3) - 3e7*@var{y}(2)^2; -## @var{y}(1) + @var{y}(2) + @var{y}(3) - 1]; -## endfunction -## opt = odeset ("Mass", [1 0 0; 0 1 0; 0 0 0], "MStateDependence", "none"); -## [@var{t},@var{y}] = ode15s (@@robertsidae, [0 1e3], [1; 0; 0], opt); -## @end group -## @end example -## @seealso{decic, odeset, odeget} -## @end deftypefn - -function varargout = ode15s (fun, trange, y0, varargin) - - solver = 'ode15s'; - - if (nargin < 3) - print_usage (); - endif - - ## Check fun, trange, y0, yp0 - fun = check_default_input (fun, trange, solver, y0); - - n = numel (y0); - - if (nargin > 3) - options = varargin{1}; - else - options = odeset (); - endif - - if (! isempty (options.Mass)) - if (ischar (options.Mass)) - try - options.Mass = str2func (options.Mass); - catch - warning (lasterr); - end_try_catch - if (! isa (options.Mass, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Mass"); - endif - endif - endif - - if (! isempty (options.Jacobian)) - if (ischar (options.Jacobian)) - try - options.Jacobian = str2func (options.Jacobian); - catch - warning (lasterr); - end_try_catch - if (! isa (options.Jacobian, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - endif - endif - - if (! isempty (options.OutputFcn)) - if (ischar (options.OutputFcn)) - try - options.OutputFcn = str2func (options.OutputFcn); - catch - warning (lasterr); - end_try_catch - if (! isa (options.OutputFcn, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "OutputFcn"); - endif - endif - endif - - if (! isempty (options.Events)) - if (ischar (options.Events)) - try - options.Events = str2func (options.Events); - catch - warning (lasterr); - end_try_catch - if (! isa (options.Events, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Events"); - endif - endif - endif - - - [defaults, classes, attributes] = ... - odedefaults (n, trange(1), trange(end)); - - classes = odeset (classes, 'Vectorized', {}); - attributes = odeset (attributes, 'Jacobian', {}, 'Vectorized', {}); - - options = ... - odemergeopts ("ode15s", options, defaults, - classes, attributes, solver); - - ## Mass - - options.havemassfun = false; - options.havestatedep = false; - options.havetimedep = false; - options.havemasssparse = false; - - if (! isempty (options.Mass)) - if (isa (options.Mass, "function_handle")) - options.havemassfun = true; - if (nargin (options.Mass) == 2) - options.havestatedep = true; - M = options.Mass (trange(1), y0); - options.havemasssparse = issparse (M); - if (any (size (M) != [n n]) || ! isnumeric (M) || ! isreal (M)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Mass"); - endif - elseif (nargin (options.Mass) == 1) - options.havetimedep = true; - M = options.Mass (trange(1)); - options.havemasssparse = issparse (M); - if (any (size (M) != [n n]) || ! isnumeric (M) || ! isreal (M)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Mass"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Mass"); - endif - elseif (ismatrix (options.Mass)) - options.havemasssparse = issparse (options.Mass); - if (any (size (options.Mass) != [n n]) || - ! isnumeric (options.Mass) || ! isreal (options.Mass)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Mass"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Mass"); - endif - endif - - - ## Jacobian - options.havejac = false; - options.havejacsparse = false; - options.havejacfun = false; - - if (! isempty (options.Jacobian)) - options.havejac = true; - if (isa (options.Jacobian, "function_handle")) - options.havejacfun = true; - if (nargin (options.Jacobian) == 2) - [A] = options.Jacobian (trange(1), y0); - if (issparse (A)) - options.havejacsparse = true; ## Jac is sparse fun - endif - if (any (size (A) != [n n]) || ! isnumeric (A) || ! isreal (A)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - elseif (ismatrix (options.Jacobian)) - if (issparse (options.Jacobian)) - options.havejacsparse = true; ## Jac is sparse matrix - endif - if (! issquare (options.Jacobian)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - else - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "Jacobian"); - endif - endif - - - ## Derivative of M(t,y) for implicit problem not implemented yet - if (! isempty (options.Mass) && ! isempty (options.Jacobian)) - if (options.MStateDependence != "none" || options.havestatedep == true) - options.havejac = false; - options.Jacobian = []; - warning ("ode15s:mass_state_dependent_provided", - ["with MStateDependence != 'none' an internal", ... - " approximation of Jacobian Matrix will be used.", ... - " Set MStateDependence equal to 'none' if you want ", ... - " to provide a constant or time-dependent Jacobian"]); - endif - endif - - ## Use sparse methods only if all matrices are sparse - if (! options.havemasssparse) - options.havejacsparse = false; - endif - - ## If Mass or Jacobian is fun, then new Jacobian is fun - if (options.havejac) - if (options.havejacfun || options.havetimedep) - options.Jacobian = @ (t, y, yp) wrapjacfun (t, y, yp, - options.Jacobian, - options.Mass, - options.havetimedep, - options.havejacfun); - options.havejacfun = true; - else ## All matrices are constant - options.Jacobian = {[- options.Jacobian], [options.Mass]}; - - endif - endif - - - - ## Abstol and Reltol - - options.haveabstolvec = false; - - if (numel (options.AbsTol) != 1 && numel (options.AbsTol) != n) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "AbsTol"); - - elseif (numel (options.AbsTol) == n) - options.haveabstolvec = true; - endif - - ## Stats - options.havestats = false; - if (strcmp (options.Stats, "on")) - options.havestats = true; - endif - - ## Don't use Refine when the output is a structure - if (nargout == 1) - options.Refine = 1; - endif - - ## OutputFcn and OutputSel - if (isempty (options.OutputFcn) && nargout == 0) - options.OutputFcn = @odeplot; - options.haveoutputfunction = true; - else - options.haveoutputfunction = ! isempty (options.OutputFcn); - endif - - options.haveoutputselection = ! isempty (options.OutputSel); - if (options.haveoutputselection) - options.OutputSel = options.OutputSel - 1; - endif - - ## Events - options.haveeventfunction = ! isempty (options.Events); - - yp0 = options.InitialSlope; - - - [t, y, te, ye, ie] = __ode15__ (@ (t, y, yp) wrap (t, y, yp, options.Mass, - options.havetimedep, - options.havestatedep, - fun), - trange, y0, yp0, options); - - if (nargout == 2) - varargout{1} = t; - varargout{2} = y; - elseif (nargout == 1) - varargout{1}.x = t; # Time stamps are saved in field x - varargout{1}.y = y; # Results are saved in field y - varargout{1}.solver = solver; - if (options.haveeventfunction) - varargout{1}.xe = te; # Time info when an event occurred - varargout{1}.ye = ye; # Results when an event occurred - varargout{1}.ie = ie; # Index info which event occurred - endif - elseif (nargout == 5) - varargout = cell (1,5); - varargout{1} = t; - varargout{2} = y; - if (options.haveeventfunction) - varargout{3} = te; # Time info when an event occurred - varargout{4} = ye; # Results when an event occurred - varargout{5} = ie; # Index info which event occurred - endif - endif - -endfunction - -function res = wrap (t, y, yp, Mass, havetimedep, havestatedep, fun) - - if (! isempty (Mass) && havestatedep) - res = Mass (t, y) * yp - fun (t, y); - elseif (! isempty (Mass) && havetimedep) - res = Mass (t) * yp - fun (t, y); - elseif (! isempty (Mass)) - res = Mass * yp - fun (t, y); - else - res = yp - fun (t, y); - endif - -endfunction - -function [jac, jact] = wrapjacfun (t, y, yp, Jac, Mass, - havetimedep, havejacfun) - if (havejacfun) - jac = - Jac (t, y); - else - jac = - Jac; - endif - - if (! isempty (Mass) && havetimedep) - jact = Mass (t); - elseif (! isempty (Mass)) - jact = Mass; - else - jact = speye (numel (y)); - endif - -endfunction - - -%!demo -%! -%! ##Solve Robertson's equations with ode15s -%! fun = @ (t, y) [-0.04*y(1) + 1e4*y(2).*y(3); -%! 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2; -%! y(1) + y(2) + y(3) - 1 ]; -%! -%! y0 = [1; 0; 0]; -%! tspan = [0 4*logspace(-6, 6)]; -%! M = [1 0 0; 0 1 0; 0 0 0]; -%! -%! options = odeset ('RelTol', 1e-4, 'AbsTol', [1e-6 1e-10 1e-6], -%! 'MStateDependence', 'none', 'Mass', M); -%! -%! [t, y] = ode15s (fun, tspan, y0, options); -%! -%! y (:,2) = 1e4 * y (:,2); -%! figure (2); -%! semilogx (t, y, 'o') -%! xlabel ('time'); -%! ylabel ('species concentration'); -%! title ('Robertson DAE problem with a Conservation Law'); -%! legend ('y1', 'y2', 'y3'); - -%!function ydot = fpol (t, y) # The Van der Pol -%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; -%!endfunction -%!function ref = fref () # The computed reference sol -%! ref = [0.32331666704577, -1.83297456798624]; -%!endfunction -%!function jac = fjac (t, y) # its Jacobian -%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; -%!endfunction -%!function jac = fjcc (t, y) # sparse type -%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]); -%!endfunction -%!function mas = fmas (t, y) -%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests -%!endfunction -%!function mas = fmsa (t, y) -%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix -%!endfunction -%!function res = rob (t, y) -%! res = [-0.04*y(1) + 1e4*y(2).*y(3); -%! 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2; -%! y(1) + y(2) + y(3) - 1 ]; -%!endfunction -%! -%!function refrob = frefrob() -%! refrob = [100, 0.617234887614937, 0.000006153591397, 0.382758958793666]; -%!endfunction -%!function [val, isterminal, direction] = feve (t, y) -%! isterminal = [0 1]; -%! if (t < 1e1) -%! val = [-1, -2]; -%! else -%! val = [1 3]; -%! endif -%! -%! direction = [1 0]; -%!endfunction -%!function masrob = massdensefunstate (t, y) -%! masrob = [1 0 0; 0 1 0; 0 0 0]; -%!endfunction -%!function masrob = masssparsefunstate (t, y) -%! masrob = sparse([1 0 0; 0 1 0; 0 0 0]); -%!endfunction -%!function masrob = massdensefuntime (t) -%! masrob = [1 0 0; 0 1 0; 0 0 0]; -%!endfunction -%!function masrob = masssparsefuntime (t) -%! masrob = sparse([1 0 0; 0 1 0; 0 0 0]); -%!endfunction -%!function jac = jacfundense (t, y) -%! jac = [-0.04, 1e4*y(3), 1e4*y(2); -%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); -%! 1, 1, 1]; -%!endfunction -%!function jac = jacfunsparse (t, y) -%! jac = sparse([-0.04, 1e4*y(3), 1e4*y(2); -%! 0.04, -1e4*y(3)-6e7*y(2), -1e4*y(2); -%! 1, 1, 1]); -%!endfunction - - - -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", [1 0 0; 0 1 0; 0 0 0]); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", sparse([1 0 0; 0 1 0; 0 0 0])); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @massdensefunstate); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @masssparsefunstate); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", 'massdensefuntime'); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", [1 0 0; 0 1 0; 0 0 0], -%! "Jacobian", 'jacfundense'); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", sparse([1 0 0; 0 1 0; 0 0 0]), -%! "Jacobian", @jacfundense); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @massdensefunstate, -%! "Jacobian", @jacfundense); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @masssparsefunstate, -%! "Jacobian", @jacfundense); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @massdensefuntime, -%! "Jacobian", @jacfundense); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%! opt = odeset ("MStateDependence", "none", -%! "Mass", 'masssparsefuntime', -%! "Jacobian", 'jacfundense'); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", [1 0 0; 0 1 0; 0 0 0], -%! "Jacobian", @jacfunsparse); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", sparse([1 0 0; 0 1 0; 0 0 0]), -%! "Jacobian", @jacfunsparse); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @massdensefunstate, -%! "Jacobian", @jacfunsparse); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! warning ("off", "ode15s:mass_state_dependent_provided", "local"); -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @masssparsefunstate, -%! "Jacobian", @jacfunsparse); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @massdensefuntime, -%! "Jacobian", @jacfunsparse); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test -%! opt = odeset ("MStateDependence", "none", -%! "Mass", @masssparsefuntime, -%! "Jacobian", @jacfunsparse); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!test # two output arguments -%! [t, y] = ode15s (@fpol, [0 2], [2 0]); -%! assert ([t(end), y(end,:)], [2, fref], 1e-2); -%!test # anonymous function instead of real function -%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; -%! [t, y] = ode15s (fvdb, [0 2], [2 0]); -%! assert ([t(end), y(end,:)], [2, fref], 1e-2); -%!test # Solve another anonymous function below zero -%! ref = [0, 14.77810590694212]; -%! [t, y] = ode15s (@(t,y) y, [-2 0], 2); -%! assert ([t(end), y(end,:)], ref, 5e-2); -%!test # InitialStep option -%! opt = odeset ("InitialStep", 1e-8); -%! [t, y] = ode15s (@fpol, [0 0.2], [2 0], opt); -%! assert ([t(2)-t(1)], [1e-8], 1e-9); -%!test # MaxStep option -%! opt = odeset ("MaxStep", 1e-3); -%! sol = ode15s (@fpol, [0 0.2], [2 0], opt); -%! assert ([sol.x(5)-sol.x(4)], [1e-3], 1e-3); -%!test # Solve in backward direction starting at t=0 -%! ref = [-1.205364552835178, 0.951542399860817]; -%! sol = ode15s (@fpol, [0 -2], [2 0]); -%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 5e-3); -%!test # Solve in backward direction starting at t=2 -%! ref = [-1.205364552835178, 0.951542399860817]; -%! sol = ode15s (@fpol, [2 0 -2], fref); -%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 3e-2); -%!test # Solve another anonymous function in backward direction -%! ref = [-1, 0.367879437558975]; -%! sol = ode15s (@(t,y) y, [0 -1], 1); -%! assert ([sol.x(end), sol.y(end,:)], ref, 1e-2); -%!test # Solve another anonymous function below zero -%! ref = [0, 14.77810590694212]; -%! sol = ode15s (@(t,y) y, [-2 0], 2); -%! assert ([sol.x(end), sol.y(end,:)], ref, 5e-2); -%!test # Solve in backward direction starting at t=0 with MaxStep option -%! ref = [-1.205364552835178, 0.951542399860817]; -%! opt = odeset ("MaxStep", 1e-3); -%! sol = ode15s (@fpol, [0 -2], [2 0], opt); -%! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3); -%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 1e-3); -%!test # AbsTol option -%! opt = odeset ("AbsTol", 1e-5); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 4e-3); -%!test # AbsTol and RelTol option -%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); -%!test # RelTol option -- higher accuracy -%! opt = odeset ("RelTol", 1e-8); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-4); -%!test # Mass option as function -%! opt = odeset ("Mass", @fmas, "MStateDependence", "none"); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); -%!test # Mass option as matrix -%! opt = odeset ("Mass", eye (2,2), "MStateDependence", "none"); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); -%!test # Mass option as sparse matrix -%! opt = odeset ("Mass", sparse (eye (2,2)), "MStateDependence", "none"); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); -%!test # Mass option as function and sparse matrix -%! opt = odeset ("Mass", 'fmsa', "MStateDependence", "none"); -%! sol = ode15s (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 3e-3); -%!test # Refine -%! opt2 = odeset ("Refine", 3, "Mass", @massdensefunstate, -%! "MStateDependence", "none"); -%! opt1 = odeset ("Mass", @massdensefunstate, "MStateDependence", "none"); -%! [t, y] = ode15s (@rob,[0 100], [1;0;0], opt1); -%! [t2, y2] = ode15s (@rob,[0 100], [1;0;0], opt2); -%! assert ([numel(t2)], numel(t)*3, 3); -%!test # Refine ignored if numel (trange) > 2 -%! opt2 = odeset ("Refine", 3, "Mass", 'massdensefunstate', -%! "MStateDependence", "none"); -%! opt1 = odeset ("Mass", @massdensefunstate, "MStateDependence", "none"); -%! [t, y] = ode15s ('rob',[0 10 100], [1;0;0], opt1); -%! [t2, y2] = ode15s ('rob',[0 10 100], [1;0;0], opt2); -%! assert ([numel(t2)], numel(t)); -%!test # Events option add further elements in sol -%! opt = odeset ("Events", @feve, "Mass", @massdensefunstate, -%! "MStateDependence", "none"); -%! sol = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert (isfield (sol, "ie")); -%! assert (sol.ie, [0;1]); -%! assert (isfield (sol, "xe")); -%! assert (isfield (sol, "ye")); -%! assert (sol.x(end), 10, 1); -%!test # Events option, five output arguments -%! opt = odeset ("Events", @feve, "Mass", @massdensefunstate, -%! "MStateDependence", "none"); -%! [t, y, te, ye, ie] = ode15s (@rob,[0 100], [1;0;0], opt); -%! assert ([t(end), te', ie'], [10, 10, 10, 0, 1], [1, 0.5, 0.5, 0, 0]); - - - - - - - - - diff --git a/scripts/ode/private/check_default_input.m b/scripts/ode/private/check_default_input.m deleted file mode 100644 --- a/scripts/ode/private/check_default_input.m +++ /dev/null @@ -1,65 +0,0 @@ -## Copyright (C) 2016, Francesco Faccio -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -function [fun] = check_default_input (fun, trange, solver, varargin); - - ## Check fun - validateattributes (fun, {"function_handle", "char"}, {}, solver, "fun"); - - if (! (nargin (fun) == nargin - 2)) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "fun"); - endif - - if (ischar (fun)) - try - fun = str2func (fun); - catch - warning (lasterr); - end_try_catch - endif - if (! isa (fun, "function_handle")) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "fun"); - endif - - ## Check trange - validateattributes (trange, {"float"}, {"vector", "real"}, solver, "trange"); - - if (numel (trange) < 2) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "trange"); - elseif (! ((all (diff (trange) > 0)) || all (diff (-trange) > 0))) - error ("Octave:invalid-input-arg", - [solver ": invalid value assigned to field '%s'"], "trange"); - endif - - ## Check y0 and yp0 - y0 = varargin{1}; - - if (nargin == 5) - yp0 = varargin{2}; - n = numel (feval (fun, trange(1), y0, yp0)); - validateattributes (yp0, {"float"}, {"numel", n}, solver, "yp0"); - else - n = numel (feval (fun, trange (1), y0)); - endif - - validateattributes (y0, {"float"}, {"numel", n}, solver, "y0"); - -endfunction diff --git a/src/mkoctfile.in.cc b/src/mkoctfile.in.cc --- a/src/mkoctfile.in.cc +++ b/src/mkoctfile.in.cc @@ -198,9 +198,6 @@ vars["FFTW3F_LDFLAGS"] = get_variable ("FFTW3F_LDFLAGS", %OCTAVE_CONF_FFTW3F_LDFLAGS%); vars["FFTW3F_LIBS"] = get_variable ("FFTW3F_LIBS", %OCTAVE_CONF_FFTW3F_LIBS%); - vars["KLU_LDFLAGS"] = get_variable ("KLU_LDFLAGS", - %OCTAVE_CONF_KLU_LDFLAGS%); - vars["KLU_LIBS"] = get_variable ("KLU_LIBS", %OCTAVE_CONF_KLU_LIBS%); vars["LIBS"] = get_variable ("LIBS", %OCTAVE_CONF_LIBS%); vars["FLIBS"] = get_variable ("FLIBS", %OCTAVE_CONF_FLIBS%); vars["OCTAVE_LINK_DEPS"] = get_variable ("OCTAVE_LINK_DEPS", @@ -218,13 +215,7 @@ vars["LFLAGS"] = get_variable ("LFLAGS", DEFAULT_LFLAGS); vars["F77_INTEGER8_FLAG"] = get_variable ("F77_INTEGER8_FLAG", %OCTAVE_CONF_F77_INTEGER_8_FLAG%); - vars["SUNDIALS_IDA_LDFLAGS"] = get_variable ("SUNDIALS_IDA_LDFLAGS", - %OCTAVE_CONF_SUNDIALS_IDA_LDFLAGS%); - vars["SUNDIALS_IDA_LIBS"] = get_variable ("SUNDIALS_IDA_LIBS", %OCTAVE_CONF_SUNDIALS_IDA_LIBS%); - vars["SUNDIALS_NVECSERIAL_LDFLAGS"] = get_variable ("SUNDIALS_NVECSERIAL_LDFLAGS", - %OCTAVE_CONF_SUNDIALS_NVECSERIAL_LDFLAGS%); - vars["SUNDIALS_NVECSERIAL_LIBS"] = get_variable ("SUNDIALS_NVECSERIAL_LIBS", %OCTAVE_CONF_SUNDIALS_NVECSERIAL_LIBS%); - + vars["ALL_FFLAGS"] = vars["FFLAGS"] + " " + vars["F77_INTEGER8_FLAG"]; vars["ALL_CFLAGS"] = vars["INCFLAGS"] + " " + vars["XTRA_CFLAGS"] + " " @@ -241,10 +232,6 @@ vars["FFTW_LIBS"] = vars["FFTW3_LDFLAGS"] + " " + vars["FFTW3_LIBS"] + " " + vars["FFTW3F_LDFLAGS"] + " " + vars["FFTW3F_LIBS"]; - - vars["SUNDIALS_LIBS"] = vars["SUNDIALS_IDA_LDFLAGS"] + " " + vars["SUNDIALS_IDA_LIBS"] + " " - + vars["KLU_LDFLAGS"] + " " + vars["KLU_LIBS"] + " " - + vars["SUNDIALS_NVECSERIAL_LDFLAGS"] + " " + vars["SUNDIALS_NVECSERIAL_LIBS"]; } static std::string usage_msg = "usage: mkoctfile [options] file ...";