diff -up octave-5.0.91/configure.ac.sundials3 octave-5.0.91/configure.ac --- octave-5.0.91/configure.ac.sundials3 2019-02-04 10:50:20.000000000 -0700 +++ octave-5.0.91/configure.ac 2019-02-05 22:05:44.096260529 -0700 @@ -2220,15 +2220,15 @@ OCTAVE_CHECK_LIB(sundials_ida, [SUNDIALS [], [don't use SUNDIALS IDA library, solvers ode15i and ode15s will be disabled], [warn_sundials_ida= OCTAVE_CHECK_SUNDIALS_SIZEOF_REALTYPE - OCTAVE_CHECK_SUNDIALS_IDA_DENSE - OCTAVE_CHECK_SUNDIALS_IDAKLU]) + OCTAVE_CHECK_SUNDIALS_SUNLINSOL_DENSE + OCTAVE_CHECK_SUNDIALS_SUNLINSOL_KLU]) LIBS="$save_LIBS" dnl Define this way instead of with an #if in oct-conf-post.h so that dnl the build features script will get the correct value. if test -n "$SUNDIALS_IDA_LIBS" \ && test -n "$SUNDIALS_NVECSERIAL_LIBS" \ - && test $octave_cv_sundials_ida_dense = yes \ + && test $octave_cv_sundials_sunlinsol_dense = yes \ && test $octave_cv_sundials_realtype_is_double = yes; then AC_DEFINE(HAVE_SUNDIALS, 1, [Define to 1 if SUNDIALS is available.]) fi diff -up octave-5.0.91/libinterp/dldfcn/__ode15__.cc.sundials3 octave-5.0.91/libinterp/dldfcn/__ode15__.cc --- octave-5.0.91/libinterp/dldfcn/__ode15__.cc.sundials3 2019-02-04 10:50:20.000000000 -0700 +++ octave-5.0.91/libinterp/dldfcn/__ode15__.cc 2019-02-05 22:06:48.074012827 -0700 @@ -1,6 +1,7 @@ /* Copyright (C) 2016-2019 Francesco Faccio +Copyright (C) 2018 William Greene This file is part of Octave. @@ -44,15 +45,34 @@ along with Octave; see the file COPYING. # include # endif -# if defined (HAVE_IDA_IDA_DENSE_H) -# include +# if defined (HAVE_SUNDIALS_SUNDIALS_MATRIX_H) +# include # endif -# if defined (HAVE_IDA_IDA_KLU_H) -# include +# if defined (HAVE_SUNDIALS_SUNDIALS_LINEARSOLVER_H) +# include +# endif + +# if defined (HAVE_SUNLINSOL_SUNLINSOL_DENSE_H) +# include +# endif + +# if defined (HAVE_IDA_IDA_DIRECT_H) +# include +# endif + +# if defined (HAVE_SUNDIALS_SUNDIALS_SPARSE_H) # include # endif +# if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H) +# include +# endif + +# if defined (HAVE_SUNMATRIX_SUNMATRIX_SPARSE_H) +# include +# endif + # if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H) # include # endif @@ -112,7 +132,8 @@ namespace octave havejacsparse (false), mem (nullptr), num (), ida_fun (nullptr), ida_jac (nullptr), dfdy (nullptr), dfdyp (nullptr), spdfdy (nullptr), spdfdyp (nullptr), fun (nullptr), jacfun (nullptr), jacspfun (nullptr), - jacdcell (nullptr), jacspcell (nullptr) + jacdcell (nullptr), jacspcell (nullptr), + sunJacMatrix (nullptr), sunLinearSolver (nullptr) { } @@ -122,11 +143,17 @@ namespace octave havejacsparse (false), mem (nullptr), num (), ida_fun (ida_fcn), ida_jac (nullptr), dfdy (nullptr), dfdyp (nullptr), spdfdy (nullptr), spdfdyp (nullptr), fun (daefun), jacfun (nullptr), jacspfun (nullptr), - jacdcell (nullptr), jacspcell (nullptr) + jacdcell (nullptr), jacspcell (nullptr), + sunJacMatrix (nullptr), sunLinearSolver (nullptr) { } - ~IDA (void) { IDAFree (&mem); } + ~IDA (void) + { + IDAFree (&mem); + SUNLinSolFree(sunLinearSolver); + SUNMatDestroy(sunJacMatrix); + } IDA& set_jacobian (octave_function *jac, DAEJacFuncDense j) @@ -184,7 +211,7 @@ namespace octave static N_Vector ColToNVec (const ColumnVector& data, long int n); void - set_up (void); + set_up (const ColumnVector& y); void set_tolerance (ColumnVector& abstol, realtype reltol); @@ -199,25 +226,24 @@ namespace octave 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, DlsMat JJ, void *user_data, + jacdense (realtype t, realtype cj, N_Vector yy, + N_Vector yyp, N_Vector, SUNMatrix JJ, void *user_data, N_Vector, N_Vector, N_Vector) { IDA *self = static_cast (user_data); - self->jacdense_impl (Neq, t, cj, yy, yyp, JJ); + self->jacdense_impl (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); + jacdense_impl (realtype t, realtype cj, + N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ); -# if defined (HAVE_SUNDIALS_IDAKLU) +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) static int jacsparse (realtype t, realtype cj, N_Vector yy, N_Vector yyp, - N_Vector, SlsMat Jac, void *user_data, N_Vector, + N_Vector, SUNMatrix Jac, void *user_data, N_Vector, N_Vector, N_Vector) { IDA *self = static_cast (user_data); @@ -227,7 +253,7 @@ namespace octave void jacsparse_impl (realtype t, realtype cj, N_Vector& yy, - N_Vector& yyp, SlsMat& Jac); + N_Vector& yyp, SUNMatrix& Jac); #endif void set_maxstep (realtype maxstep); @@ -291,6 +317,8 @@ namespace octave DAEJacFuncSparse jacspfun; DAEJacCellDense jacdcell; DAEJacCellSparse jacspcell; + SUNMatrix sunJacMatrix; + SUNLinearSolver sunLinearSolver; }; int @@ -323,36 +351,61 @@ namespace octave } void - IDA::set_up (void) + IDA::set_up (const ColumnVector& y) { + N_Vector yy = ColToNVec(y, num); + if (havejacsparse) { -# if defined (HAVE_SUNDIALS_IDAKLU) - if (IDAKLU (mem, num, num*num, CSC_MAT) != 0) - error ("IDAKLU solver not initialized"); +#if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) + + sunJacMatrix = SUNSparseMatrix (num, num, num*num, CSC_MAT); + if (! sunJacMatrix) + error ("Unable to create sparse Jacobian for Sundials"); + + sunLinearSolver = SUNKLU (yy, sunJacMatrix); + if (! sunLinearSolver) + error ("Unable to create KLU sparse solver"); + + if (IDADlsSetLinearSolver (mem, sunLinearSolver, sunJacMatrix)) + error ("Unable to set sparse linear solver"); + + IDADlsSetJacFn(mem, IDA::jacsparse); - IDASlsSetSparseJacFn (mem, IDA::jacsparse); # else - error ("IDAKLU is not available in this version of Octave"); + error ("SUNDIALS SUNLINSOL KLU is not available in this version of Octave"); # endif + } else { - if (IDADense (mem, num) != 0) - error ("IDADense solver not initialized"); - if (havejac && IDADlsSetDenseJacFn (mem, IDA::jacdense) != 0) - error ("Dense Jacobian not set"); + sunJacMatrix = SUNDenseMatrix (num, num); + if (! sunJacMatrix) + error ("Unable to create dense Jacobian for Sundials"); + + sunLinearSolver = SUNDenseLinearSolver(yy, sunJacMatrix); + if (! sunLinearSolver) + error ("Unable to create dense linear solver"); + + if (IDADlsSetLinearSolver (mem, sunLinearSolver, sunJacMatrix)) + error ("Unable to set dense linear solver"); + + if (havejac && IDADlsSetJacFn (mem, IDA::jacdense) != 0) + error("Unable to set dense Jacobian function"); + } } void - IDA::jacdense_impl (long int Neq, realtype t, realtype cj, - N_Vector& yy, N_Vector& yyp, DlsMat& JJ) + IDA::jacdense_impl (realtype t, realtype cj, + N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; + long int Neq = NV_LENGTH_S(yy); + ColumnVector y = NVecToCol (yy, Neq); ColumnVector yp = NVecToCol (yyp, Neq); @@ -366,15 +419,15 @@ namespace octave std::copy (jac.fortran_vec (), jac.fortran_vec () + jac.numel (), - JJ->data); + SUNDenseMatrix_Data(JJ)); END_INTERRUPT_WITH_EXCEPTIONS; } -# if defined (HAVE_SUNDIALS_IDAKLU) +# if defined (HAVE_SUNDIALS_SUNLINSOL_KLU) void IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp, - SlsMat& Jac) + SUNMatrix& Jac) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; @@ -390,17 +443,18 @@ namespace octave else jac = (*jacspcell) (spdfdy, spdfdyp, cj); - SparseSetMatToZero (Jac); - int *colptrs = *(Jac->colptrs); - int *rowvals = *(Jac->rowvals); + SUNMatZero_Sparse (Jac); + sunindextype *colptrs = SUNSparseMatrix_IndexPointers (Jac); + sunindextype *rowvals = SUNSparseMatrix_IndexValues (Jac); for (int i = 0; i < num + 1; i++) colptrs[i] = jac.cidx(i); + double *d = SUNSparseMatrix_Data (Jac); for (int i = 0; i < jac.nnz (); i++) { rowvals[i] = jac.ridx(i); - Jac->data[i] = jac.data(i); + d[i] = jac.data(i); } END_INTERRUPT_WITH_EXCEPTIONS; @@ -567,7 +621,7 @@ namespace octave //main loop while (((posdirection == 1 && tsol < tend) - || (posdirection == 0 && tsol > tend)) + || (posdirection == 0 && tsol > tend)) && status == 0) { if (IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0) @@ -692,7 +746,7 @@ namespace octave // Linear interpolation ie(0) = index(0); te(0) = tsol - val (index(0)) * (tsol - told) - / (val (index(0)) - oldval (index(0))); + / (val (index(0)) - oldval (index(0))); ColumnVector ytemp = y - ((tsol - te(0)) * (y - yold) / (tsol - told)); @@ -717,7 +771,7 @@ namespace octave // Linear interpolation ie(temp+i) = index(i); te(temp+i) = tsol - val(index(i)) * (tsol - told) - / (val(index(i)) - oldval(index(i))); + / (val(index(i)) - oldval(index(i))); ColumnVector ytemp = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told); @@ -1096,7 +1150,7 @@ namespace octave event_fcn = options.getfield("Events").function_value (); // Set up linear solver - dae.set_up (); + dae.set_up (y0); // Integrate retval = dae.integrate (numt, tspan, y0, yp0, refine, diff -up octave-5.0.91/m4/acinclude.m4.sundials3 octave-5.0.91/m4/acinclude.m4 --- octave-5.0.91/m4/acinclude.m4.sundials3 2019-02-04 10:50:20.000000000 -0700 +++ octave-5.0.91/m4/acinclude.m4 2019-02-05 22:05:44.100260576 -0700 @@ -2210,14 +2210,11 @@ dnl Check whether SUNDIALS IDA library i dnl precision realtype. dnl AC_DEFUN([OCTAVE_CHECK_SUNDIALS_SIZEOF_REALTYPE], [ - AC_CHECK_HEADERS([ida/ida.h ida.h]) AC_CACHE_CHECK([whether SUNDIALS IDA is configured with double precision realtype], [octave_cv_sundials_realtype_is_double], [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ #if defined (HAVE_IDA_IDA_H) #include - #else - #include #endif #include ]], [[ @@ -2233,61 +2230,72 @@ AC_DEFUN([OCTAVE_CHECK_SUNDIALS_SIZEOF_R fi ]) dnl -dnl Check whether SUNDIALS IDA library is configured with IDAKLU +dnl Check whether SUNDIALS IDA library is configured with SUNLINSOL_KLU dnl enabled. dnl -AC_DEFUN([OCTAVE_CHECK_SUNDIALS_IDAKLU], [ - AC_CHECK_HEADERS([ida/ida_klu.h ida_klu.h]) - AC_CACHE_CHECK([whether SUNDIALS IDA is configured with IDAKLU enabled], - [octave_cv_sundials_idaklu], +AC_DEFUN([OCTAVE_CHECK_SUNDIALS_SUNLINSOL_KLU], [ + AC_CHECK_HEADERS([sundials/sundials_sparse.h sunlinsol/sunlinsol_klu.h sunmatrix/sunmatrix_sparse.h]) + AC_CACHE_CHECK([whether SUNDIALS IDA is configured with SUNLINSOL_KLU enabled], + [octave_cv_sundials_sunlinsol_klu], [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ - #if defined (HAVE_IDA_IDA_KLU_H) - #include - #else - #include + #if defined (HAVE_IDA_IDA_H) + #include + #endif + #if defined (HAVE_SUNDIALS_SUNDIALS_SPARSE_H) + #include + #endif + #if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H) + #include #endif ]], [[ - IDAKLU (0, 0, 0, 0); + SUNKLU (0, 0); ]])], - octave_cv_sundials_idaklu=yes, - octave_cv_sundials_idaklu=no) + octave_cv_sundials_sunlinsol_klu=yes, + octave_cv_sundials_sunlinsol_klu=no) ]) - if test $octave_cv_sundials_idaklu = yes; then - AC_DEFINE(HAVE_SUNDIALS_IDAKLU, 1, - [Define to 1 if SUNDIALS IDA is configured with IDAKLU enabled.]) + if test $octave_cv_sundials_sunlinsol_klu = yes; then + AC_DEFINE(HAVE_SUNDIALS_SUNLINSOL_KLU, 1, + [Define to 1 if SUNDIALS IDA is configured with SUNLINSOL_KLU enabled.]) else - warn_sundials_idaklu="SUNDIALS IDA library not configured with IDAKLU, ode15i and ode15s will not support the sparse Jacobian feature" - OCTAVE_CONFIGURE_WARNING([warn_sundials_idaklu]) + warn_sundials_idaklu="SUNDIALS IDA library not configured with SUNLINSOL_KLU, ode15i and ode15s will not support the sparse Jacobian feature" + OCTAVE_CONFIGURE_WARNING([warn_sundials_sunlinsol_klu]) fi ]) dnl -dnl Check whether SUNDIALS IDA library has the IDADENSE linear solver. +dnl Check whether SUNDIALS IDA library has the SUNLINSOL_DENSE linear solver. dnl The IDADENSE API was removed in SUNDIALS version 3.0.0. dnl -AC_DEFUN([OCTAVE_CHECK_SUNDIALS_IDA_DENSE], [ - AC_CHECK_HEADERS([ida/ida_dense.h ida_dense.h]) - AC_CACHE_CHECK([whether SUNDIALS IDA includes the IDADENSE linear solver], - [octave_cv_sundials_ida_dense], +AC_DEFUN([OCTAVE_CHECK_SUNDIALS_SUNLINSOL_DENSE], [ + AC_CHECK_HEADERS([sunlinsol/sunlinsol_dense.h sundials/sundials_matrix.h sundials/sundials_linearsolver.h ida/ida_direct.h]) + AC_CACHE_CHECK([whether SUNDIALS IDA includes the SUNLINSOL_DENSE linear solver], + [octave_cv_sundials_sunlinsol_dense], [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ - #if defined (HAVE_IDA_IDA_DENSE_H) - #include - #else - #include + #if defined (HAVE_IDA_IDA_H) + #include + #endif + #if defined (HAVE_SUNDIALS_SUNDIALS_MATRIX_H) + #include + #endif + #if defined (HAVE_SUNDIALS_SUNDIALS_LINEARSOLVER_H) + #include #endif + #if defined (HAVE_IDA_IDA_DIRECT_H) + #include + #endif ]], [[ void *mem = 0; long int num = 0; IDADense (mem, num); ]])], - octave_cv_sundials_ida_dense=yes, - octave_cv_sundials_ida_dense=no) + octave_cv_sundials_sunlinsol_dense=yes, + octave_cv_sundials_sunlinsol_dense=no) ]) - if test $octave_cv_sundials_ida_dense = yes; then - AC_DEFINE(HAVE_SUNDIALS_IDADENSE, 1, - [Define to 1 if SUNDIALS IDA includes the IDADENSE linear solver.]) + if test $octave_cv_sundials_sunlinsol_dense = yes; then + AC_DEFINE(HAVE_SUNDIALS_SUNLINSOL_DENSE, 1, + [Define to 1 if SUNDIALS IDA includes the SUNLINSOL_DENSE linear solver.]) else - warn_sundials_ida_dense="SUNDIALS IDA library does not include the IDADENSE linear solver, ode15i and ode15s will be disabled" - OCTAVE_CONFIGURE_WARNING([warn_sundials_ida_dense]) + warn_sundials_ida_dense="SUNDIALS IDA library does not include the SUNLINSOL_DENSE linear solver, ode15i and ode15s will be disabled" + OCTAVE_CONFIGURE_WARNING([warn_sundials_sunlinsol_dense]) fi ]) dnl diff -up octave-5.0.91/scripts/ode/ode15i.m.sundials3 octave-5.0.91/scripts/ode/ode15i.m --- octave-5.0.91/scripts/ode/ode15i.m.sundials3 2019-02-04 10:50:20.000000000 -0700 +++ octave-5.0.91/scripts/ode/ode15i.m 2019-02-05 22:05:44.101260588 -0700 @@ -452,7 +452,7 @@ endfunction %! assert ([t(end), y(end,:)], fref, 1e-3); ## Jacobian fun sparse -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! 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); @@ -545,7 +545,7 @@ endfunction %! "invalid value assigned to field 'Jacobian'"); ## Jacobian cell sparse wrong dimension -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! DFDY = sparse ([-0.04, 1; %! 0.04, 1]); %! DFDYP = sparse ([-1, 0, 0; diff -up octave-5.0.91/scripts/ode/ode15s.m.sundials3 octave-5.0.91/scripts/ode/ode15s.m --- octave-5.0.91/scripts/ode/ode15s.m.sundials3 2019-02-04 10:50:20.000000000 -0700 +++ octave-5.0.91/scripts/ode/ode15s.m 2019-02-05 22:05:44.102260599 -0700 @@ -545,21 +545,21 @@ endfunction %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! 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); -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! 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); -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! warning ("off", "ode15s:mass_state_dependent_provided", "local"); %! opt = odeset ("MStateDependence", "none", %! "Mass", @massdensefunstate, @@ -575,14 +575,14 @@ endfunction %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! 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); -%!testif HAVE_SUNDIALS_IDAKLU +%!testif HAVE_SUNDIALS_SUNLINSOL_KLU %! opt = odeset ("MStateDependence", "none", %! "Mass", @masssparsefuntime, %! "Jacobian", @jacfunsparse);