////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2016-2020 The Octave Project Developers
//
// See the file COPYRIGHT.md in the top-level directory of this
// distribution or <https://octave.org/copyright/>.
//
// 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
// <https://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include "dColVector.h"
#include "dMatrix.h"
#include "dSparse.h"

#include "Cell.h"
#include "defun-dld.h"
#include "error.h"
#include "errwarn.h"
#include "oct-map.h"
#include "ov.h"
#include "ovl.h"
#include "pager.h"
#include "parse.h"

#if defined (HAVE_SUNDIALS)

#  if defined (HAVE_NVECTOR_NVECTOR_SERIAL_H)
#    include <nvector/nvector_serial.h>
#  endif

#  if defined (HAVE_IDA_IDA_H)
#    include <ida/ida.h>
#  elif defined (HAVE_IDA_H)
#    include <ida.h>
#  endif
#  if defined (HAVE_IDA_IDA_DIRECT_H)
#    include <ida/ida_direct.h>
#  elif defined (HAVE_IDA_DIRECT_H)
#    include <ida_direct.h>
#  endif

#  if defined (HAVE_SUNLINSOL_SUNLINSOL_DENSE_H)
#    include <sunlinsol/sunlinsol_dense.h>
#  endif

#  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

namespace octave
{
#  if ! defined (HAVE_IDASETJACFN) && defined (HAVE_IDADLSSETJACFN)
  static inline int
  IDASetJacFn (void *ida_mem, IDADlsJacFn jac)
  {
    return IDADlsSetJacFn (ida_mem, jac);
  }
#  endif

#  if ! defined (HAVE_IDASETLINEARSOLVER) && defined (HAVE_IDADLSSETLINEARSOLVER)
  static inline int
  IDASetLinearSolver (void *ida_mem, SUNLinearSolver LS, SUNMatrix A)
  {
    return IDADlsSetLinearSolver (ida_mem, LS, A);
  }
#  endif

#  if ! defined (HAVE_SUNLINSOL_DENSE) && defined (HAVE_SUNDENSELINEARSOLVER)
  static inline SUNLinearSolver
  SUNLinSol_Dense (N_Vector y, SUNMatrix A)
  {
    return SUNDenseLinearSolver (y, A);
  }
#  endif

#  if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
#    if ! defined (HAVE_SUNLINSOL_KLU) && defined (HAVE_SUNKLU)
  static inline SUNLinearSolver
  SUNLinSol_KLU (N_Vector y, SUNMatrix A)
  {
    return SUNKLU (y, A);
  }
#    endif
#  endif

  static inline realtype *
  nv_data_s (N_Vector& v)
  {
#if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC)
    // Disable warning from GCC about old-style casts in Sundials
    // macro expansions.  Do this in a function so that this
    // diagnostic may still be enabled for the rest of the file.
#  pragma GCC diagnostic push
#  pragma GCC diagnostic ignored "-Wold-style-cast"
#endif

    return NV_DATA_S (v);

#if defined (HAVE_PRAGMA_GCC_DIAGNOSTIC)
    // Restore prevailing warning state for remainder of the file.
#  pragma GCC diagnostic pop
#endif
  }

  class IDA
  {
  public:

    typedef
    ColumnVector (*DAERHSFuncIDA) (const ColumnVector& x,
                                   const ColumnVector& xdot,
                                   realtype t, const octave_value& idaf);

    typedef
    Matrix (*DAEJacFuncDense) (const ColumnVector& x,
                               const ColumnVector& xdot, realtype t,
                               realtype cj, const octave_value& idaj);

    typedef
    SparseMatrix (*DAEJacFuncSparse) (const ColumnVector& x,
                                      const ColumnVector& xdot,
                                      realtype t, realtype cj,
                                      const octave_value& idaj);

    typedef
    Matrix (*DAEJacCellDense) (Matrix *dfdy, Matrix *dfdyp,
                               realtype cj);

    typedef
    SparseMatrix (*DAEJacCellSparse) (SparseMatrix *dfdy,
                                      SparseMatrix *dfdyp, realtype cj);

    //Default
    IDA (void)
      : m_t0 (0.0), m_y0 (), m_yp0 (), m_havejac (false), m_havejacfun (false),
        m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fun (),
        m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr),
        m_spdfdyp (nullptr), m_fun (nullptr), m_jacfun (nullptr), m_jacspfun (nullptr),
        m_jacdcell (nullptr), m_jacspcell (nullptr),
        m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr)
    { }


    IDA (realtype t, ColumnVector y, ColumnVector yp,
         const octave_value& ida_fcn, DAERHSFuncIDA daefun)
      : m_t0 (t), m_y0 (y), m_yp0 (yp), m_havejac (false), m_havejacfun (false),
        m_havejacsparse (false), m_mem (nullptr), m_num (), m_ida_fun (ida_fcn),
        m_ida_jac (), m_dfdy (nullptr), m_dfdyp (nullptr), m_spdfdy (nullptr),
        m_spdfdyp (nullptr), m_fun (daefun), m_jacfun (nullptr), m_jacspfun (nullptr),
        m_jacdcell (nullptr), m_jacspcell (nullptr),
        m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr)
    { }


    ~IDA (void)
    {
      IDAFree (&m_mem);
      SUNLinSolFree (m_sunLinearSolver);
      SUNMatDestroy (m_sunJacMatrix);
    }

    IDA&
    set_jacobian (const octave_value& jac, DAEJacFuncDense j)
    {
      m_jacfun = j;
      m_ida_jac = jac;
      m_havejac = true;
      m_havejacfun = true;
      m_havejacsparse = false;

      return *this;
    }

    IDA&
    set_jacobian (const octave_value& jac, DAEJacFuncSparse j)
    {
      m_jacspfun = j;
      m_ida_jac = jac;
      m_havejac = true;
      m_havejacfun = true;
      m_havejacsparse = true;

      return *this;
    }

    IDA&
    set_jacobian (Matrix *dy, Matrix *dyp, DAEJacCellDense j)
    {
      m_jacdcell = j;
      m_dfdy = dy;
      m_dfdyp = dyp;
      m_havejac = true;
      m_havejacfun = false;
      m_havejacsparse = false;

      return *this;
    }

    IDA&
    set_jacobian (SparseMatrix *dy, SparseMatrix *dyp,
                  DAEJacCellSparse j)
    {
      m_jacspcell = j;
      m_spdfdy = dy;
      m_spdfdyp = dyp;
      m_havejac = true;
      m_havejacfun = false;
      m_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 (const ColumnVector& y);

    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 (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 (t, cj, yy, yyp, JJ);
      return 0;
    }

    void
    jacdense_impl (realtype t, realtype cj,
                   N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ);

#  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,
               N_Vector, N_Vector)
    {
      IDA *self = static_cast <IDA *> (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, SUNMatrix& Jac);
#  endif

    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, const octave_value& output_fcn,
                 ColumnVector& outputsel, bool haveeventfunction,
                 const octave_value& event_fcn, ColumnVector& te,
                 Matrix& ye, ColumnVector& ie, ColumnVector& oldval,
                 ColumnVector& oldisterminal, ColumnVector& olddir,
                 int& temp, ColumnVector& yold);

    bool
    outputfun (const octave_value& output_fcn, bool haveoutputsel,
               const ColumnVector& output, realtype tout, realtype tend,
               ColumnVector& outputsel, const std::string& flag);


    bool
    event (const octave_value& event_fcn,
           ColumnVector& te, Matrix& ye, ColumnVector& ie,
           realtype tsol, const ColumnVector& y, const 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,
               const octave_value& output_fcn, bool haveoutputsel,
               ColumnVector& outputsel, bool haveeventfunction,
               const octave_value& event_fcn);

    void print_stat (void);

  private:

    realtype m_t0;
    ColumnVector m_y0;
    ColumnVector m_yp0;
    bool m_havejac;
    bool m_havejacfun;
    bool m_havejacsparse;
    void *m_mem;
    int m_num;
    octave_value m_ida_fun;
    octave_value m_ida_jac;
    Matrix *m_dfdy;
    Matrix *m_dfdyp;
    SparseMatrix *m_spdfdy;
    SparseMatrix *m_spdfdyp;
    DAERHSFuncIDA m_fun;
    DAEJacFuncDense m_jacfun;
    DAEJacFuncSparse m_jacspfun;
    DAEJacCellDense m_jacdcell;
    DAEJacCellSparse m_jacspcell;
    SUNMatrix m_sunJacMatrix;
    SUNLinearSolver m_sunLinearSolver;
  };

  int
  IDA::resfun (realtype t, N_Vector yy, N_Vector yyp, N_Vector rr,
               void *user_data)
  {
    IDA *self = static_cast <IDA *> (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)
  {
    ColumnVector y = IDA::NVecToCol (yy, m_num);

    ColumnVector yp = IDA::NVecToCol (yyp, m_num);

    ColumnVector res = (*m_fun) (y, yp, t, m_ida_fun);

    realtype *puntrr = nv_data_s (rr);

    for (octave_idx_type i = 0; i < m_num; i++)
      puntrr[i] = res(i);
  }

  void
  IDA::set_up (const ColumnVector& y)
  {
    N_Vector yy = ColToNVec(y, m_num);

    if (m_havejacsparse)
      {
#  if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
        // FIXME : one should not allocate space for a full Jacobian
        // when using a sparse format.  Consider allocating less space
        // then possibly using SUNSparseMatrixReallocate to increase it.
        m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, m_num*m_num, CSC_MAT);
        if (! m_sunJacMatrix)
          error ("Unable to create sparse Jacobian for Sundials");

        m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix);
        if (! m_sunLinearSolver)
          error ("Unable to create KLU sparse solver");

        if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
          error ("Unable to set sparse linear solver");

        IDASetJacFn (m_mem, IDA::jacsparse);

#  else
        error ("SUNDIALS SUNLINSOL KLU is not available in this version of Octave");
#  endif

      }
    else
      {

        m_sunJacMatrix = SUNDenseMatrix (m_num, m_num);
        if (! m_sunJacMatrix)
          error ("Unable to create dense Jacobian for Sundials");

        m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix);
        if (! m_sunLinearSolver)
          error ("Unable to create dense linear solver");

        if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
          error ("Unable to set dense linear solver");

        if (m_havejac && IDASetJacFn (m_mem, IDA::jacdense) != 0)
          error ("Unable to set dense Jacobian function");

      }
  }

  void
  IDA::jacdense_impl (realtype t, realtype cj,
                      N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ)

  {
    long int Neq = NV_LENGTH_S(yy);

    ColumnVector y = NVecToCol (yy, Neq);

    ColumnVector yp = NVecToCol (yyp, Neq);

    Matrix jac;

    if (m_havejacfun)
      jac = (*m_jacfun) (y, yp, t, cj, m_ida_jac);
    else
      jac = (*m_jacdcell) (m_dfdy, m_dfdyp, cj);

    std::copy (jac.fortran_vec (),
               jac.fortran_vec () + jac.numel (),
               SUNDenseMatrix_Data (JJ));
  }

#  if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
  void
  IDA::jacsparse_impl (realtype t, realtype cj, N_Vector& yy, N_Vector& yyp,
                       SUNMatrix& Jac)

  {
    ColumnVector y = NVecToCol (yy, m_num);

    ColumnVector yp = NVecToCol (yyp, m_num);

    SparseMatrix jac;

    if (m_havejacfun)
      jac = (*m_jacspfun) (y, yp, t, cj, m_ida_jac);
    else
      jac = (*m_jacspcell) (m_spdfdy, m_spdfdyp, cj);

    SUNMatZero_Sparse (Jac);
    sunindextype *colptrs = SUNSparseMatrix_IndexPointers (Jac);
    sunindextype *rowvals = SUNSparseMatrix_IndexValues (Jac);

    for (int i = 0; i < m_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);
        d[i] = jac.data(i);
      }
  }
#  endif

  ColumnVector
  IDA::NVecToCol (N_Vector& v, long int n)
  {
    ColumnVector data (n);
    realtype *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 = 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;

    if (IDASetUserData (m_mem, userdata) != 0)
      error ("User data not set");
  }

  void
  IDA::initialize (void)
  {
    m_num = m_y0.numel ();
    m_mem = IDACreate ();

    N_Vector yy = ColToNVec (m_y0, m_num);

    N_Vector yyp = ColToNVec (m_yp0, m_num);

    IDA::set_userdata ();

    if (IDAInit (m_mem, IDA::resfun, m_t0, yy, yyp) != 0)
      error ("IDA not initialized");
  }

  void
  IDA::set_tolerance (ColumnVector& abstol, realtype reltol)
  {
    N_Vector abs_tol = ColToNVec (abstol, m_num);

    if (IDASVtolerances (m_mem, reltol, abs_tol) != 0)
      error ("IDA: Tolerance not set");

    N_VDestroy_Serial (abs_tol);
  }

  void
  IDA::set_tolerance (realtype abstol, realtype reltol)
  {
    if (IDASStolerances (m_mem, reltol, abstol) != 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,
                  const octave_value& output_fcn, bool haveoutputsel,
                  ColumnVector& outputsel, bool haveeventfunction,
                  const octave_value& event_fcn)
  {
    // Set up output
    ColumnVector tout, yout (m_num), ypout (m_num), ysel (outputsel.numel ());
    ColumnVector ie, te, oldval, oldisterminal, olddir;
    Matrix output, ye;
    int cont = 0, temp = 0;
    bool status = 0;
    std::string string = "";
    ColumnVector yold = y;

    realtype tsol = tspan(0);
    realtype tend = tspan(numt-1);

    N_Vector yyp = ColToNVec (yp, m_num);

    N_Vector yy = ColToNVec (y, m_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, m_num);

        for (octave_idx_type i = 0; i < m_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)

            if (IDASolve (m_mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0)
              error ("IDASolve failed");

            yout = NVecToCol (yy, m_num);
            ypout = NVecToCol (yyp, m_num);
            tout(j) = tsol;

            for (octave_idx_type i = 0; i < m_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, m_num);
                tout.resize (j + 1);
              }

          }
      }
    else // numel (tspan) == 2
      {
        // First output value
        tout.resize (1);
        tout(0) = tsol;
        output.resize (1, m_num);

        for (octave_idx_type i = 0; i < m_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)
          {
            if (IDASolve (m_mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 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, m_num);
            cont += 1;
            output.resize (cont + 1, m_num); // This may be not efficient
            tout.resize (cont + 1);
            tout(cont) = tsol;
            yout = NVecToCol (yy, m_num);

            for (octave_idx_type i = 0; i < m_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 (m_num);

            if (IDAGetDky (m_mem, tend, 0, dky) != 0)
              error ("IDA failed to interpolate y");

            tout(cont) = tend;
            yout = NVecToCol (dky, m_num);

            for (octave_idx_type i = 0; i < m_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");

      }

    return ovl (tout, output, te, ye, ie);
  }

  bool
  IDA::event (const octave_value& event_fcn,
              ColumnVector& te, Matrix& ye, ColumnVector& ie,
              realtype tsol, const ColumnVector& y, const 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 = ovl (tsol, y, yp);

    // cont is the number of steps reached by the solver
    // temp is the number of events registered

    if (flag == "init")
      {
        octave_value_list 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);
        octave_value_list output = feval (event_fcn, args, 3);
        ColumnVector val = output(0).vector_value ();
        ColumnVector isterminal = output(1).vector_value ();
        ColumnVector 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, m_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 < m_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 (), m_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 < m_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, const octave_value& output_fcn,
                    ColumnVector& outputsel, bool haveeventfunction,
                    const octave_value& 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 (m_num);

    N_Vector dkyp = N_VNew_Serial (m_num);

    ColumnVector yout (m_num);
    ColumnVector ypout (m_num);
    std::string string = "";

    if (IDAGetLastStep (m_mem, &h) != 0)
      error ("IDA failed to return last step");

    if (IDAGetCurrentTime (m_mem, &tcur) != 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++)
      {
        if (IDAGetDky (m_mem, tin + step*i, 0, dky) != 0)
          error ("IDA failed to interpolate y");

        if (IDAGetDky (m_mem, tin + step*i, 1, dkyp) != 0)
          error ("IDA failed to interpolate yp");

        cont += 1;
        output.resize (cont + 1, m_num);
        tout.resize (cont + 1);

        tout(cont) = tin + step * i;
        yout = NVecToCol (dky, m_num);
        ypout = NVecToCol (dkyp, m_num);

        for (octave_idx_type j = 0; j < m_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 (const octave_value& output_fcn, bool haveoutputsel,
                  const ColumnVector& yout, realtype tsol,
                  realtype tend, ColumnVector& outputsel,
                  const std::string& flag)
  {
    bool status = 0;

    octave_value_list output;
    output(2) = flag;

    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;

    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;
        octave_value_list 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)
  {
    if (IDASetMaxStep (m_mem, maxstep) != 0)
      error ("IDA: Max Step not set");
  }

  void
  IDA::set_initialstep (realtype initialstep)
  {
    if (IDASetInitStep (m_mem, initialstep) != 0)
      error ("IDA: Initial Step not set");
  }

  void
  IDA::set_maxorder (int maxorder)
  {
    if (IDASetMaxOrd (m_mem, maxorder) != 0)
      error ("IDA: Max Order not set");
  }

  void
  IDA::print_stat (void)
  {
    long int nsteps = 0, netfails = 0, nrevals = 0;

    if (IDAGetNumSteps (m_mem, &nsteps) != 0)
      error ("IDA failed to return the number of internal steps");

    if (IDAGetNumErrTestFails (m_mem, &netfails) != 0)
      error ("IDA failed to return the number of internal errors");

    if (IDAGetNumResEvals (m_mem, &nrevals) != 0)
      error ("IDA failed to return the number of residual evaluations");

    octave_stdout << nsteps << " successful steps\n";
    octave_stdout << netfails << " failed attempts\n";
    octave_stdout << nrevals << " function evaluations\n";
    // octave_stdout << " partial derivatives\n";
    // octave_stdout << " LU decompositions\n";
    // octave_stdout << " solutions of linear systems\n";
  }

  ColumnVector
  ida_user_function (const ColumnVector& x, const ColumnVector& xdot,
                     double t, const octave_value& ida_fc)
  {
    octave_value_list tmp;

    try
      {
        tmp = feval (ida_fc, ovl (t, x, xdot), 1);
      }
    catch (execution_exception& e)
      {
        err_user_supplied_eval (e, "__ode15__");
      }

    return tmp(0).vector_value ();
  }

  Matrix
  ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot,
                 double t, double cj, const octave_value& ida_jc)
  {
    octave_value_list tmp;

    try
      {
        tmp = feval (ida_jc, ovl (t, x, xdot), 2);
      }
    catch (execution_exception& e)
      {
        err_user_supplied_eval (e, "__ode15__");
      }

    return tmp(0).matrix_value () + cj * tmp(1).matrix_value ();
  }

  SparseMatrix
  ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot,
                  double t, double cj, const octave_value& ida_jc)
  {
    octave_value_list tmp;

    try
      {
        tmp = feval (ida_jc, ovl (t, x, xdot), 2);
      }
    catch (execution_exception& e)
      {
        err_user_supplied_eval (e, "__ode15__");
      }

    return tmp(0).sparse_matrix_value () + cj * tmp(1).sparse_matrix_value ();
  }

  Matrix
  ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj)
  {
    return (*dfdy) + cj * (*dfdyp);
  }

  SparseMatrix
  ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp,
                       double cj)
  {
    return (*spdfdy) + cj * (*spdfdyp);
  }

  octave_value_list
  do_ode15 (const octave_value& 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;
    SparseMatrix ida_spdfdy, ida_spdfdyp;

    if (havejac)
      {
        if (havejacfun)
          {
            octave_value ida_jac = options.getfield ("Jacobian");

            if (havejacsparse)
              dae.set_jacobian (ida_jac, ida_sparse_jac);
            else
              dae.set_jacobian (ida_jac, ida_dense_jac);
          }
        else
          {
            Cell jaccell = options.getfield ("Jacobian").cell_value ();

            if (havejacsparse)
              {
                ida_spdfdy = jaccell(0).sparse_matrix_value ();
                ida_spdfdyp = jaccell(1).sparse_matrix_value ();

                dae.set_jacobian (&ida_spdfdy, &ida_spdfdyp,
                                  ida_sparse_cell_jac);
              }
            else
              {
                ida_dfdy = jaccell(0).matrix_value ();
                ida_dfdyp = jaccell(1).matrix_value ();

                dae.set_jacobian (&ida_dfdy, &ida_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").isempty ())
      {
        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_value output_fcn;
    ColumnVector outputsel;

    // OutputFcn
    bool haveoutputfunction
      = options.getfield("haveoutputfunction").bool_value ();

    if (haveoutputfunction)
      output_fcn = options.getfield("OutputFcn");

    // OutputSel
    bool haveoutputsel = options.getfield("haveoutputselection").bool_value ();

    if (haveoutputsel)
      outputsel = options.getfield("OutputSel").vector_value ();

    octave_value event_fcn;

    // Events
    bool haveeventfunction
      = options.getfield("haveeventfunction").bool_value ();

    if (haveeventfunction)
      event_fcn = options.getfield("Events");

    // Set up linear solver
    dae.set_up (y0);

    // 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, ,
           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)
    print_usage ();

  // Check odefun
  octave_value ida_fcn = args(0);

  if (! ida_fcn.is_function_handle ())
    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.issorted () == UNSORTED || 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).isstruct ())
    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
}

/*
## No test needed for internal helper function.
%!assert (1)
*/