main/science/mathematics/octave/files/octave-sundials4.patch

# HG changeset patch
# User Bill Greene <w.h.greene@gmail.com>
# Date 1550229868 -3600
#      Fri Feb 15 12:24:28 2019 +0100
# Node ID 4bf27c090f5695bcf545fc4af15e2a61a3941d46
# Parent  f034b29320ad5034ad5c66480f64411e9d773440
Update DAE/IDE solvers to work with SUNDIALS 3 (bug #52475).

* libinterp/dldfcn/__ode15__.cc : use SUNDIALS API version 3.x

diff --git a/libinterp/dldfcn/__ode15__.cc b/libinterp/dldfcn/__ode15__.cc
--- a/libinterp/dldfcn/__ode15__.cc
+++ b/libinterp/dldfcn/__ode15__.cc
@@ -1,6 +1,7 @@
 /*
 
 Copyright (C) 2016-2019 Francesco Faccio <francesco.faccio@mail.polimi.it>
+Copyright (C) 2019 William Greene <w.h.greene@gmail.com>
 
 This file is part of Octave.
 
@@ -112,7 +113,8 @@
         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 +124,17 @@
         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 +192,7 @@
     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 +207,24 @@
     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)
     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 +234,7 @@
 
     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 +298,8 @@
     DAEJacFuncSparse jacspfun;
     DAEJacCellDense jacdcell;
     DAEJacCellSparse jacspcell;
+    SUNMatrix sunJacMatrix;
+    SUNLinearSolver sunLinearSolver;
   };
 
   int
@@ -323,36 +332,61 @@
   }
 
   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");
+
+        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");
 
-        IDASlsSetSparseJacFn (mem, IDA::jacsparse);
+        if (IDADlsSetLinearSolver (mem, sunLinearSolver, sunJacMatrix))
+          error ("Unable to set sparse linear solver");
+
+        IDADlsSetJacFn (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");
+
+        sunJacMatrix = SUNDenseMatrix (num, num);
+        if (! sunJacMatrix)
+          error ("Unable to create dense Jacobian for Sundials");
 
-        if (havejac && IDADlsSetDenseJacFn (mem, IDA::jacdense) != 0)
-          error ("Dense Jacobian not set");
+        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,7 +400,7 @@
 
     std::copy (jac.fortran_vec (),
                jac.fortran_vec () + jac.numel (),
-               JJ->data);
+      SUNDenseMatrix_Data(JJ));
 
     END_INTERRUPT_WITH_EXCEPTIONS;
   }
@@ -374,7 +408,7 @@
 #  if defined (HAVE_SUNDIALS_IDAKLU)
   void
   IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp,
-                       SlsMat& Jac)
+                       SUNMatrix& Jac)
 
   {
     BEGIN_INTERRUPT_WITH_EXCEPTIONS;
@@ -390,17 +424,18 @@
     else
       jac = (*jacspcell) (spdfdy, spdfdyp, cj);
 
-    SparseSetMatToZero (Jac);
-    int *colptrs = *(Jac->colptrs);
-    int *rowvals = *(Jac->rowvals);
+    SUNMatZero_Sparse (Jac);
+    octave_idx_type *colptrs = SUNSparseMatrix_IndexPointers (Jac);
+    octave_idx_type *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 +602,7 @@
 
         //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 +727,7 @@
             // 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 +752,7 @@
                 // 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 +1131,7 @@
       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,
# HG changeset patch
# User Carlo de Falco <carlo.defalco@polimi.it>
# Date 1550230515 -3600
#      Fri Feb 15 12:35:15 2019 +0100
# Node ID 9b27b77d2fec1a334baa7f4b2ec70751f4bf2b99
# Parent  4bf27c090f5695bcf545fc4af15e2a61a3941d46
Update detection of sundials in the build system (bug #52475).

* m4/acnclude.m4 : check for sunlinsol_klu.h and sunlinsol_dense.h
* configure.ac : update check for sundials features
* script/ode/ode15i.m : update conditionals in tests
* script/ode/ode15s.m : update conditionals in tests
* libinterp/dldfcn/__ode15__.cc : udate conditionals in preprocessor directives

diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -2206,15 +2206,15 @@
   [], [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 --git a/libinterp/dldfcn/__ode15__.cc b/libinterp/dldfcn/__ode15__.cc
--- a/libinterp/dldfcn/__ode15__.cc
+++ b/libinterp/dldfcn/__ode15__.cc
@@ -45,15 +45,31 @@
 #    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_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_KLU_H)
-#    include <ida/ida_klu.h>
+#  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_NVECTOR_NVECTOR_SERIAL_H)
 #    include <nvector/nvector_serial.h>
 #  endif
@@ -221,7 +237,7 @@
     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, SUNMatrix Jac, void *user_data, N_Vector,
@@ -338,7 +354,7 @@
 
     if (havejacsparse)
       {
-#  if defined (HAVE_SUNDIALS_IDAKLU)
+#if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
 
         sunJacMatrix = SUNSparseMatrix (num, num, num*num, CSC_MAT);
         if (! sunJacMatrix)
@@ -405,7 +421,7 @@
     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,
                        SUNMatrix& Jac)
diff --git a/m4/acinclude.m4 b/m4/acinclude.m4
--- a/m4/acinclude.m4
+++ b/m4/acinclude.m4
@@ -2210,14 +2210,11 @@
 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 @@
   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])
+  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 --git a/scripts/ode/ode15i.m b/scripts/ode/ode15i.m
--- a/scripts/ode/ode15i.m
+++ b/scripts/ode/ode15i.m
@@ -452,7 +452,7 @@
 %! 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 @@
 %!       "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 --git a/scripts/ode/ode15s.m b/scripts/ode/ode15s.m
--- a/scripts/ode/ode15s.m
+++ b/scripts/ode/ode15s.m
@@ -545,21 +545,21 @@
 %! [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 @@
 %! [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);
# HG changeset patch
# Parent  6ef7a85c3b1908f57665a08d1d7205c15018c2af

diff --git a/libinterp/dldfcn/__ode15__.cc b/libinterp/dldfcn/__ode15__.cc
--- a/libinterp/dldfcn/__ode15__.cc
+++ b/libinterp/dldfcn/__ode15__.cc
@@ -67,6 +67,18 @@
 
 
 #  if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
+#    if defined (HAVE_KLU_H)
+#      include <klu.h>
+#    endif
+#    if defined (HAVE_KLU_KLU_H)
+#      include <klu/klu.h>
+#    endif
+#    if defined (HAVE_SUITESPARSE_KLU_H)
+#      include <suitesparse/klu.h>
+#    endif
+#    if defined (HAVE_UFPARSE_KLU_H)
+#      include <ufsparse/klu.h>
+#    endif
 #    include <sunlinsol/sunlinsol_klu.h>
 #  endif
 
diff --git a/m4/acinclude.m4 b/m4/acinclude.m4
--- a/m4/acinclude.m4
+++ b/m4/acinclude.m4
@@ -2244,6 +2244,18 @@
          #if defined (HAVE_SUNDIALS_SUNDIALS_SPARSE_H)
          #include <sundials/sundials_sparse.h>
          #endif
+         #if defined (HAVE_KLU_H)
+         #include <klu.h>
+         #endif
+         #if defined (HAVE_KLU_KLU_H)
+         #include <klu/klu.h>
+         #endif
+         #if defined (HAVE_SUITESPARSE_KLU_H)
+         #include <suitesparse/klu.h>
+         #endif
+         #if defined (HAVE_UFPARSE_KLU_H)
+         #include <ufsparse/klu.h>
+         #endif
          #if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
          #include <sunlinsol/sunlinsol_klu.h>
          #endif