Mercurial > octave

changeset 27182:0beeb6817376

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
__ode15__.cc style fixes * __ode15__.cc (class IDA): Use m_prefix for member variables. Minor style fixes.
author John W. Eaton <jwe@octave.org>
date Thu, 13 Jun 2019 13:49:53 -0500
parents 5754e8d8199a
children 2d9decd77e58
files libinterp/dldfcn/__ode15__.cc
diffstat 1 files changed, 151 insertions(+), 150 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/dldfcn/__ode15__.cc	Thu Jun 13 09:23:11 2019 -0500
+++ b/libinterp/dldfcn/__ode15__.cc	Thu Jun 13 13:49:53 2019 -0500
@@ -161,64 +161,67 @@
 
     //Default
     IDA (void)
-      : t0 (0.0), y0 (), yp0 (), havejac (false), havejacfun (false),
-        havejacsparse (false), mem (nullptr), num (), ida_fun (),
-        ida_jac (), dfdy (nullptr), dfdyp (nullptr), spdfdy (nullptr),
-        spdfdyp (nullptr), fun (nullptr), jacfun (nullptr), jacspfun (nullptr),
-        jacdcell (nullptr), jacspcell (nullptr),
-        sunJacMatrix (nullptr), sunLinearSolver (nullptr)
+      : 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)
-      : t0 (t), y0 (y), yp0 (yp), havejac (false), havejacfun (false),
-        havejacsparse (false), mem (nullptr), num (), ida_fun (ida_fcn),
-        ida_jac (), dfdy (nullptr), dfdyp (nullptr), spdfdy (nullptr),
-        spdfdyp (nullptr), fun (daefun), jacfun (nullptr), jacspfun (nullptr),
-        jacdcell (nullptr), jacspcell (nullptr),
-        sunJacMatrix (nullptr), sunLinearSolver (nullptr)
+      : 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 (&mem);
-      SUNLinSolFree(sunLinearSolver);
-      SUNMatDestroy(sunJacMatrix);
+      IDAFree (&m_mem);
+      SUNLinSolFree (m_sunLinearSolver);
+      SUNMatDestroy (m_sunJacMatrix);
     }
 
     IDA&
     set_jacobian (const octave_value& jac, DAEJacFuncDense j)
     {
-      jacfun = j;
-      ida_jac = jac;
-      havejac = true;
-      havejacfun = true;
-      havejacsparse = false;
+      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)
     {
-      jacspfun = j;
-      ida_jac = jac;
-      havejac = true;
-      havejacfun = true;
-      havejacsparse = true;
+      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)
     {
-      jacdcell = j;
-      dfdy = dy;
-      dfdyp = dyp;
-      havejac = true;
-      havejacfun = false;
-      havejacsparse = false;
+      m_jacdcell = j;
+      m_dfdy = dy;
+      m_dfdyp = dyp;
+      m_havejac = true;
+      m_havejacfun = false;
+      m_havejacsparse = false;
+
       return *this;
     }
 
@@ -226,12 +229,13 @@
     set_jacobian (SparseMatrix *dy, SparseMatrix *dyp,
                   DAEJacCellSparse j)
     {
-      jacspcell = j;
-      spdfdy = dy;
-      spdfdyp = dyp;
-      havejac = true;
-      havejacfun = false;
-      havejacsparse = true;
+      m_jacspcell = j;
+      m_spdfdy = dy;
+      m_spdfdyp = dyp;
+      m_havejac = true;
+      m_havejacfun = false;
+      m_havejacsparse = true;
+
       return *this;
     }
 
@@ -331,27 +335,27 @@
 
   private:
 
-    realtype t0;
-    ColumnVector y0;
-    ColumnVector yp0;
-    bool havejac;
-    bool havejacfun;
-    bool havejacsparse;
-    void *mem;
-    int num;
-    octave_value ida_fun;
-    octave_value ida_jac;
-    Matrix *dfdy;
-    Matrix *dfdyp;
-    SparseMatrix *spdfdy;
-    SparseMatrix *spdfdyp;
-    DAERHSFuncIDA fun;
-    DAEJacFuncDense jacfun;
-    DAEJacFuncSparse jacspfun;
-    DAEJacCellDense jacdcell;
-    DAEJacCellSparse jacspcell;
-    SUNMatrix sunJacMatrix;
-    SUNLinearSolver sunLinearSolver;
+    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
@@ -367,41 +371,41 @@
   IDA::resfun_impl (realtype t, N_Vector& yy,
                     N_Vector& yyp, N_Vector& rr)
   {
-    ColumnVector y = IDA::NVecToCol (yy, num);
+    ColumnVector y = IDA::NVecToCol (yy, m_num);
 
-    ColumnVector yp = IDA::NVecToCol (yyp, num);
+    ColumnVector yp = IDA::NVecToCol (yyp, m_num);
 
-    ColumnVector res = (*fun) (y, yp, t, ida_fun);
+    ColumnVector res = (*m_fun) (y, yp, t, m_ida_fun);
 
     realtype *puntrr = nv_data_s (rr);
 
-    for (octave_idx_type i = 0; i < num; i++)
+    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, num);
+    N_Vector yy = ColToNVec(y, m_num);
 
-    if (havejacsparse)
+    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.
-        sunJacMatrix = SUNSparseMatrix (num, num, num*num, CSC_MAT);
-        if (! sunJacMatrix)
+        m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, m_num*m_num, CSC_MAT);
+        if (! m_sunJacMatrix)
           error ("Unable to create sparse Jacobian for Sundials");
 
-        sunLinearSolver = SUNLinSol_KLU (yy, sunJacMatrix);
-        if (! sunLinearSolver)
+        m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix);
+        if (! m_sunLinearSolver)
           error ("Unable to create KLU sparse solver");
 
-        if (IDASetLinearSolver (mem, sunLinearSolver, sunJacMatrix))
+        if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
           error ("Unable to set sparse linear solver");
 
-        IDASetJacFn (mem, IDA::jacsparse);
+        IDASetJacFn (m_mem, IDA::jacsparse);
 
 #  else
         error ("SUNDIALS SUNLINSOL KLU is not available in this version of Octave");
@@ -411,18 +415,18 @@
     else
       {
 
-        sunJacMatrix = SUNDenseMatrix (num, num);
-        if (! sunJacMatrix)
+        m_sunJacMatrix = SUNDenseMatrix (m_num, m_num);
+        if (! m_sunJacMatrix)
           error ("Unable to create dense Jacobian for Sundials");
 
-        sunLinearSolver = SUNLinSol_Dense (yy, sunJacMatrix);
-        if (! sunLinearSolver)
+        m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix);
+        if (! m_sunLinearSolver)
           error ("Unable to create dense linear solver");
 
-        if (IDASetLinearSolver (mem, sunLinearSolver, sunJacMatrix))
+        if (IDASetLinearSolver (m_mem, m_sunLinearSolver, m_sunJacMatrix))
           error ("Unable to set dense linear solver");
 
-        if (havejac && IDASetJacFn (mem, IDA::jacdense) != 0)
+        if (m_havejac && IDASetJacFn (m_mem, IDA::jacdense) != 0)
           error ("Unable to set dense Jacobian function");
 
       }
@@ -441,14 +445,14 @@
 
     Matrix jac;
 
-    if (havejacfun)
-      jac = (*jacfun) (y, yp, t, cj, ida_jac);
+    if (m_havejacfun)
+      jac = (*m_jacfun) (y, yp, t, cj, m_ida_jac);
     else
-      jac = (*jacdcell) (dfdy, dfdyp, cj);
+      jac = (*m_jacdcell) (m_dfdy, m_dfdyp, cj);
 
     std::copy (jac.fortran_vec (),
                jac.fortran_vec () + jac.numel (),
-      SUNDenseMatrix_Data(JJ));
+               SUNDenseMatrix_Data (JJ));
   }
 
 #  if defined (HAVE_SUNDIALS_SUNLINSOL_KLU)
@@ -457,22 +461,22 @@
                        SUNMatrix& Jac)
 
   {
-    ColumnVector y = NVecToCol (yy, num);
+    ColumnVector y = NVecToCol (yy, m_num);
 
-    ColumnVector yp = NVecToCol (yyp, num);
+    ColumnVector yp = NVecToCol (yyp, m_num);
 
     SparseMatrix jac;
 
-    if (havejacfun)
-      jac = (*jacspfun) (y, yp, t, cj, ida_jac);
+    if (m_havejacfun)
+      jac = (*m_jacspfun) (y, yp, t, cj, m_ida_jac);
     else
-      jac = (*jacspcell) (spdfdy, spdfdyp, cj);
+      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 < num + 1; i++)
+    for (int i = 0; i < m_num + 1; i++)
       colptrs[i] = jac.cidx(i);
 
     double *d = SUNSparseMatrix_Data (Jac);
@@ -514,32 +518,32 @@
   {
     void *userdata = this;
 
-    if (IDASetUserData (mem, userdata) != 0)
+    if (IDASetUserData (m_mem, userdata) != 0)
       error ("User data not set");
   }
 
   void
   IDA::initialize (void)
   {
-    num = y0.numel ();
-    mem = IDACreate ();
+    m_num = m_y0.numel ();
+    m_mem = IDACreate ();
 
-    N_Vector yy = ColToNVec (y0, num);
+    N_Vector yy = ColToNVec (m_y0, m_num);
 
-    N_Vector yyp = ColToNVec (yp0, num);
+    N_Vector yyp = ColToNVec (m_yp0, m_num);
 
     IDA::set_userdata ();
 
-    if (IDAInit (mem, IDA::resfun, t0, yy, yyp) != 0)
+    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, num);
+    N_Vector abs_tol = ColToNVec (abstol, m_num);
 
-    if (IDASVtolerances (mem, reltol, abs_tol) != 0)
+    if (IDASVtolerances (m_mem, reltol, abs_tol) != 0)
       error ("IDA: Tolerance not set");
 
     N_VDestroy_Serial (abs_tol);
@@ -548,7 +552,7 @@
   void
   IDA::set_tolerance (realtype abstol, realtype reltol)
   {
-    if (IDASStolerances (mem, reltol, abstol) != 0)
+    if (IDASStolerances (m_mem, reltol, abstol) != 0)
       error ("IDA: Tolerance not set");
   }
 
@@ -561,7 +565,7 @@
                   const octave_value& event_fcn)
   {
     // Set up output
-    ColumnVector tout, yout (num), ypout (num), ysel (outputsel.numel ());
+    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;
@@ -572,9 +576,9 @@
     realtype tsol = tspan(0);
     realtype tend = tspan(numt-1);
 
-    N_Vector yyp = ColToNVec (yp, num);
+    N_Vector yyp = ColToNVec (yp, m_num);
 
-    N_Vector yy = ColToNVec (y, num);
+    N_Vector yy = ColToNVec (y, m_num);
 
     // Initialize OutputFcn
     if (haveoutputfcn)
@@ -592,9 +596,9 @@
         // First output value
         tout.resize (numt);
         tout(0) = tsol;
-        output.resize (numt, num);
+        output.resize (numt, m_num);
 
-        for (octave_idx_type i = 0; i < num; i++)
+        for (octave_idx_type i = 0; i < m_num; i++)
           output.elem (0, i) = y.elem (i);
 
         //Main loop
@@ -602,14 +606,14 @@
           {
             // IDANORMAL already interpolates tspan(j)
 
-            if (IDASolve (mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0)
+            if (IDASolve (m_mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0)
               error ("IDASolve failed");
 
-            yout = NVecToCol (yy, num);
-            ypout = NVecToCol (yyp, num);
+            yout = NVecToCol (yy, m_num);
+            ypout = NVecToCol (yyp, m_num);
             tout(j) = tsol;
 
-            for (octave_idx_type i = 0; i < num; i++)
+            for (octave_idx_type i = 0; i < m_num; i++)
               output.elem (j, i) = yout.elem (i);
 
             if (haveoutputfcn)
@@ -624,7 +628,7 @@
             // If integration is stopped, return only the reached steps
             if (status == 1)
               {
-                output.resize (j + 1, num);
+                output.resize (j + 1, m_num);
                 tout.resize (j + 1);
               }
 
@@ -635,9 +639,9 @@
         // First output value
         tout.resize (1);
         tout(0) = tsol;
-        output.resize (1, num);
+        output.resize (1, m_num);
 
-        for (octave_idx_type i = 0; i < num; i++)
+        for (octave_idx_type i = 0; i < m_num; i++)
           output.elem (0, i) = y.elem (i);
 
         bool posdirection = (tend > tsol);
@@ -647,7 +651,7 @@
                || (posdirection == 0 && tsol > tend))
                && status == 0)
           {
-            if (IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0)
+            if (IDASolve (m_mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0)
               error ("IDASolve failed");
 
             if (haverefine)
@@ -658,14 +662,14 @@
                                          ye, ie, oldval, oldisterminal,
                                          olddir, temp, yold);
 
-            ypout = NVecToCol (yyp, num);
+            ypout = NVecToCol (yyp, m_num);
             cont += 1;
-            output.resize (cont + 1, num); // This may be not efficient
+            output.resize (cont + 1, m_num); // This may be not efficient
             tout.resize (cont + 1);
             tout(cont) = tsol;
-            yout = NVecToCol (yy, num);
+            yout = NVecToCol (yy, m_num);
 
-            for (octave_idx_type i = 0; i < num; i++)
+            for (octave_idx_type i = 0; i < m_num; i++)
               output.elem (cont, i) = yout.elem (i);
 
             if (haveoutputfcn && ! haverefine && tout(cont) < tend)
@@ -681,15 +685,15 @@
         if (status == 0)
           {
             // Interpolate in tend
-            N_Vector dky = N_VNew_Serial (num);
+            N_Vector dky = N_VNew_Serial (m_num);
 
-            if (IDAGetDky (mem, tend, 0, dky) != 0)
+            if (IDAGetDky (m_mem, tend, 0, dky) != 0)
               error ("IDA failed to interpolate y");
 
             tout(cont) = tend;
-            yout = NVecToCol (dky, num);
+            yout = NVecToCol (dky, m_num);
 
-            for (octave_idx_type i = 0; i < num; i++)
+            for (octave_idx_type i = 0; i < m_num; i++)
               output.elem (cont, i) = yout.elem (i);
 
             // Plot final value
@@ -763,7 +767,7 @@
           {
             temp = 1; // register only the first event
             te.resize (1);
-            ye.resize (1, num);
+            ye.resize (1, m_num);
             ie.resize (1);
 
             // Linear interpolation
@@ -774,7 +778,7 @@
             ColumnVector ytemp
               = y - ((tsol - te(0)) * (y - yold) / (tsol - told));
 
-            for (octave_idx_type i = 0; i < num; i++)
+            for (octave_idx_type i = 0; i < m_num; i++)
               ye.elem (0, i) = ytemp.elem (i);
 
           }
@@ -782,7 +786,7 @@
           // Not first step: register all events and test if stop integration or not
           {
             te.resize (temp + index.numel ());
-            ye.resize (temp + index.numel (), num);
+            ye.resize (temp + index.numel (), m_num);
             ie.resize (temp + index.numel ());
 
             for (octave_idx_type i = 0; i < index.numel (); i++)
@@ -799,7 +803,7 @@
                 ColumnVector ytemp
                   = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told);
 
-                for (octave_idx_type j = 0; j < num; j++)
+                for (octave_idx_type j = 0; j < m_num; j++)
                   ye.elem (temp + i, j) = ytemp.elem (j);
 
               }
@@ -831,18 +835,18 @@
     realtype h = 0, tcur = 0;
     bool status = 0;
 
-    N_Vector dky = N_VNew_Serial (num);
+    N_Vector dky = N_VNew_Serial (m_num);
 
-    N_Vector dkyp = N_VNew_Serial (num);
+    N_Vector dkyp = N_VNew_Serial (m_num);
 
-    ColumnVector yout (num);
-    ColumnVector ypout (num);
+    ColumnVector yout (m_num);
+    ColumnVector ypout (m_num);
     std::string string = "";
 
-    if (IDAGetLastStep (mem, &h) != 0)
+    if (IDAGetLastStep (m_mem, &h) != 0)
       error ("IDA failed to return last step");
 
-    if (IDAGetCurrentTime (mem, &tcur) != 0)
+    if (IDAGetCurrentTime (m_mem, &tcur) != 0)
       error ("IDA failed to return the current time");
 
     realtype tin = tcur - h;
@@ -853,21 +857,21 @@
          i < refine && tin + step * i < tend && status == 0;
          i++)
       {
-        if (IDAGetDky (mem, tin + step*i, 0, dky) != 0)
+        if (IDAGetDky (m_mem, tin + step*i, 0, dky) != 0)
           error ("IDA failed to interpolate y");
 
-        if (IDAGetDky (mem, tin + step*i, 1, dkyp) != 0)
+        if (IDAGetDky (m_mem, tin + step*i, 1, dkyp) != 0)
           error ("IDA failed to interpolate yp");
 
         cont += 1;
-        output.resize (cont + 1, num);
+        output.resize (cont + 1, m_num);
         tout.resize (cont + 1);
 
         tout(cont) = tin + step * i;
-        yout = NVecToCol (dky, num);
-        ypout = NVecToCol (dkyp, num);
+        yout = NVecToCol (dky, m_num);
+        ypout = NVecToCol (dkyp, m_num);
 
-        for (octave_idx_type j = 0; j < num; j++)
+        for (octave_idx_type j = 0; j < m_num; j++)
           output.elem (cont, j) = yout.elem (j);
 
         if (haveoutputfcn)
@@ -936,21 +940,21 @@
   void
   IDA::set_maxstep (realtype maxstep)
   {
-    if (IDASetMaxStep (mem, maxstep) != 0)
+    if (IDASetMaxStep (m_mem, maxstep) != 0)
       error ("IDA: Max Step not set");
   }
 
   void
   IDA::set_initialstep (realtype initialstep)
   {
-    if (IDASetInitStep (mem, initialstep) != 0)
+    if (IDASetInitStep (m_mem, initialstep) != 0)
       error ("IDA: Initial Step not set");
   }
 
   void
   IDA::set_maxorder (int maxorder)
   {
-    if (IDASetMaxOrd (mem, maxorder) != 0)
+    if (IDASetMaxOrd (m_mem, maxorder) != 0)
       error ("IDA: Max Order not set");
   }
 
@@ -959,13 +963,13 @@
   {
     long int nsteps = 0, netfails = 0, nrevals = 0;
 
-    if (IDAGetNumSteps (mem, &nsteps) != 0)
+    if (IDAGetNumSteps (m_mem, &nsteps) != 0)
       error ("IDA failed to return the number of internal steps");
 
-    if (IDAGetNumErrTestFails (mem, &netfails) != 0)
+    if (IDAGetNumErrTestFails (m_mem, &netfails) != 0)
       error ("IDA failed to return the number of internal errors");
 
-    if (IDAGetNumResEvals (mem, &nrevals) != 0)
+    if (IDAGetNumResEvals (m_mem, &nrevals) != 0)
       error ("IDA failed to return the number of residual evaluations");
 
     octave_stdout << nsteps << " successful steps\n";
@@ -1064,16 +1068,14 @@
 
     bool havejacfun = options.getfield ("havejacfun").bool_value ();
 
-    Matrix ida_dfdy, ida_dfdyp, *dfdy, *dfdyp;
-    SparseMatrix ida_spdfdy, ida_spdfdyp, *spdfdy, *spdfdyp;
-    octave_value ida_jac;
-    Cell jaccell;
+    Matrix ida_dfdy, ida_dfdyp;
+    SparseMatrix ida_spdfdy, ida_spdfdyp;
 
     if (havejac)
       {
         if (havejacfun)
           {
-            ida_jac = options.getfield ("Jacobian");
+            octave_value ida_jac = options.getfield ("Jacobian");
 
             if (havejacsparse)
               dae.set_jacobian (ida_jac, ida_sparse_jac);
@@ -1082,23 +1084,22 @@
           }
         else
           {
-            jaccell = options.getfield ("Jacobian").cell_value ();
+            Cell 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);
+
+                dae.set_jacobian (&ida_spdfdy, &ida_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);
+
+                dae.set_jacobian (&ida_dfdy, &ida_dfdyp, ida_dense_cell_jac);
               }
           }
       }