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 <francesco.faccio@mail.polimi.it>
+Copyright (C) 2018 William Greene <w.h.greene@gmail.com>
This file is part of Octave.
@@ -44,15 +45,34 @@ along with Octave; see the file COPYING.
# include <ida/ida.h>
# endif
-# if defined (HAVE_IDA_IDA_DENSE_H)
-# include <ida/ida_dense.h>
+# if defined (HAVE_SUNDIALS_SUNDIALS_MATRIX_H)
+# include <sundials/sundials_matrix.h>
# endif
-# if defined (HAVE_IDA_IDA_KLU_H)
-# include <ida/ida_klu.h>
+# if defined (HAVE_SUNDIALS_SUNDIALS_LINEARSOLVER_H)
+# include <sundials/sundials_linearsolver.h>
+# endif
+
+# if defined (HAVE_SUNLINSOL_SUNLINSOL_DENSE_H)
+# include <sunlinsol/sunlinsol_dense.h>
+# endif
+
+# if defined (HAVE_IDA_IDA_DIRECT_H)
+# include <ida/ida_direct.h>
+# endif
+
+# if defined (HAVE_SUNDIALS_SUNDIALS_SPARSE_H)
# include <sundials/sundials_sparse.h>
# endif
+# if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
+# include <sunlinsol/sunlinsol_klu.h>
+# endif
+
+# if defined (HAVE_SUNMATRIX_SUNMATRIX_SPARSE_H)
+# include <sunmatrix/sunmatrix_sparse.h>
+# endif
+
# if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H)
# include <nvector/nvector_serial.h>
# 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 <IDA *> (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 <IDA *> (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 <ida/ida.h>
- #else
- #include <ida.h>
#endif
#include <assert.h>
]], [[
@@ -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 <ida/ida_klu.h>
- #else
- #include <ida_klu.h>
+ #if defined (HAVE_IDA_IDA_H)
+ #include <ida/ida.h>
+ #endif
+ #if defined (HAVE_SUNDIALS_SUNDIALS_SPARSE_H)
+ #include <sundials/sundials_sparse.h>
+ #endif
+ #if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
+ #include <sunlinsol/sunlinsol_klu.h>
#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 <ida/ida_dense.h>
- #else
- #include <ida_dense.h>
+ #if defined (HAVE_IDA_IDA_H)
+ #include <ida/ida.h>
+ #endif
+ #if defined (HAVE_SUNDIALS_SUNDIALS_MATRIX_H)
+ #include <sundials/sundials_matrix.h>
+ #endif
+ #if defined (HAVE_SUNDIALS_SUNDIALS_LINEARSOLVER_H)
+ #include <sundials/sundials_linearsolver.h>
#endif
+ #if defined (HAVE_IDA_IDA_DIRECT_H)
+ #include <ida/ida_direct.h>
+ #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);