Mercurial > octave

changeset 22917:c1ac4ac3ebbc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
* __ode15__.cc: Style fixes.
author John W. Eaton <jwe@octave.org>
date Fri, 16 Dec 2016 12:44:57 -0500
parents 1632fb20b60b
children 0b5d9978d7b1
files libinterp/dldfcn/__ode15__.cc
diffstat 1 files changed, 207 insertions(+), 371 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/dldfcn/__ode15__.cc	Fri Dec 16 12:04:27 2016 -0500
+++ b/libinterp/dldfcn/__ode15__.cc	Fri Dec 16 12:44:57 2016 -0500
@@ -104,7 +104,7 @@
         ida_jac (0), dfdy (0), dfdyp (0), spdfdy (0),
         spdfdyp (0), fun (0), jacfun (0), jacspfun (0),
         jacdcell (0), jacspcell (0)
-    { };
+    { }
 
 
     IDA (realtype t, ColumnVector y, ColumnVector yp,
@@ -114,7 +114,7 @@
         ida_jac (0), dfdy (0), dfdyp (0), spdfdy (0),
         spdfdyp (0), fun (daefun), jacfun (0), jacspfun (0),
         jacdcell (0), jacspcell (0)
-    { };
+    { }
 
 
     ~IDA (void) { IDAFree (&mem); }
@@ -166,17 +166,13 @@
       return *this;
     }
 
-    void
-    set_userdata (void);
+    void set_userdata (void);
 
-    void
-    initialize (void);
+    void initialize (void);
 
-    static ColumnVector
-    NVecToCol (N_Vector& v, long int n);
+    static ColumnVector NVecToCol (N_Vector& v, long int n);
 
-    static N_Vector
-    ColToNVec (const ColumnVector& data, long int n);
+    static N_Vector ColToNVec (const ColumnVector& data, long int n);
 
     void
     set_up (void);
@@ -200,10 +196,8 @@
               N_Vector yyp, N_Vector, DlsMat 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);
+      IDA *self = static_cast <IDA *> (user_data);
+      self->jacdense_impl (Neq, t, cj, yy, yyp, JJ);
       return 0;
     }
 
@@ -216,10 +210,8 @@
                N_Vector, SlsMat 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);
+      IDA *self = static_cast <IDA *> (user_data);
+      self->jacsparse_impl (t, cj, yy, yyp, Jac);
       return 0;
     }
 
@@ -227,11 +219,9 @@
     jacsparse_impl (realtype t, realtype cj, N_Vector& yy,
                     N_Vector& yyp, SlsMat& Jac);
 
-    void
-    set_maxstep (realtype maxstep);
+    void set_maxstep (realtype maxstep);
 
-    void
-    set_initialstep (realtype initialstep);
+    void set_initialstep (realtype initialstep);
 
     bool
     interpolate (int& cont, Matrix& output, ColumnVector& tout,
@@ -257,8 +247,7 @@
            ColumnVector& oldisterminal, ColumnVector& olddir,
            int cont, int& temp, realtype told, ColumnVector& yold);
 
-    void
-    set_maxorder (int maxorder);
+    void set_maxorder (int maxorder);
 
     octave_value_list
     integrate (const int numt, const ColumnVector& tt,
@@ -268,8 +257,7 @@
                ColumnVector& outputsel, bool haveeventfunction,
                octave_function *event_fcn);
 
-    void
-    print_stat (void);
+    void print_stat (void);
 
   private:
 
@@ -292,17 +280,14 @@
     DAEJacFuncSparse jacspfun;
     DAEJacCellDense jacdcell;
     DAEJacCellSparse jacspcell;
-
   };
 
   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);
+    IDA *self = static_cast <IDA *> (user_data);
+    self->resfun_impl (t, yy, yyp, rr);
     return 0;
   }
 
@@ -312,20 +297,16 @@
   {
     BEGIN_INTERRUPT_WITH_EXCEPTIONS;
 
-    ColumnVector y =
-      IDA::NVecToCol (yy, num);
+    ColumnVector y = IDA::NVecToCol (yy, num);
 
-    ColumnVector yp =
-      IDA::NVecToCol (yyp, num);
+    ColumnVector yp = IDA::NVecToCol (yyp, num);
 
-    ColumnVector res =
-      (*fun) (y, yp, t, ida_fun);
+    ColumnVector res = (*fun) (y, yp, t, ida_fun);
 
-    realtype *puntrr =
-      nv_data_s (rr);
+    realtype *puntrr = nv_data_s (rr);
 
     for (octave_idx_type i = 0; i < num; i++)
-      puntrr [i] = res (i);
+      puntrr[i] = res(i);
 
     END_INTERRUPT_WITH_EXCEPTIONS;
   }
@@ -335,30 +316,18 @@
   {
     if (havejacsparse)
       {
-        int flag =
-          IDAKLU (mem, num, num*num, CSC_MAT);
-        if (flag != 0)
-          {
-            error ("IDAKLU solver not initialized");
-          }
-        flag = IDASlsSetSparseJacFn (mem, IDA::jacsparse);
+        if (IDAKLU (mem, num, num*num, CSC_MAT) != 0)
+          error ("IDAKLU solver not initialized");
+
+        IDASlsSetSparseJacFn (mem, IDA::jacsparse);
       }
     else
       {
-        int flag =
-          IDADense (mem, num);
-        if (flag != 0)
-          {
-            error ("IDADense solver not initialized");
-          }
-        if (havejac)
-          {
-            flag = IDADlsSetDenseJacFn (mem, IDA::jacdense);
-            if (flag != 0)
-              {
-                error ("Dense Jacobian not set");
-              }
-          }
+        if (IDADense (mem, num) != 0)
+          error ("IDADense solver not initialized");
+
+        if (havejac && IDADlsSetDenseJacFn (mem, IDA::jacdense) != 0)
+          error ("Dense Jacobian not set");
       }
   }
 
@@ -369,11 +338,9 @@
   {
     BEGIN_INTERRUPT_WITH_EXCEPTIONS;
 
-    ColumnVector y =
-      NVecToCol (yy, Neq);
+    ColumnVector y = NVecToCol (yy, Neq);
 
-    ColumnVector yp =
-      NVecToCol (yyp, Neq);
+    ColumnVector yp = NVecToCol (yyp, Neq);
 
     Matrix jac;
 
@@ -384,7 +351,7 @@
 
     std::copy (jac.fortran_vec (),
                jac.fortran_vec () + jac.numel (),
-               JJ -> data);
+               JJ->data);
 
     END_INTERRUPT_WITH_EXCEPTIONS;
   }
@@ -396,11 +363,9 @@
   {
     BEGIN_INTERRUPT_WITH_EXCEPTIONS;
 
-    ColumnVector y =
-      NVecToCol (yy, num);
+    ColumnVector y = NVecToCol (yy, num);
 
-    ColumnVector yp =
-      NVecToCol (yyp, num);
+    ColumnVector yp = NVecToCol (yyp, num);
 
     SparseMatrix jac;
 
@@ -410,16 +375,16 @@
       jac = (*jacspcell) (spdfdy, spdfdyp, cj);
 
     SparseSetMatToZero (Jac);
-    int *colptrs = *(Jac -> colptrs);
-    int *rowvals = *(Jac -> rowvals);
+    int *colptrs = *(Jac->colptrs);
+    int *rowvals = *(Jac->rowvals);
 
     for (int i = 0; i < num + 1; i++)
       colptrs[i] = jac.cidx(i);
 
     for (int i = 0; i < jac.nnz (); i++)
       {
-        rowvals[i] = jac.ridx (i);
-        Jac -> data[i] = jac.data (i);
+        rowvals[i] = jac.ridx(i);
+        Jac->data[i] = jac.data(i);
       }
 
     END_INTERRUPT_WITH_EXCEPTIONS;
@@ -429,11 +394,10 @@
   IDA::NVecToCol (N_Vector& v, long int n)
   {
     ColumnVector data (n);
-    realtype *punt;
-    punt = nv_data_s (v);
+    realtype *punt = nv_data_s (v);
 
     for (octave_idx_type i = 0; i < n; i++)
-      data (i) = punt [i];
+      data(i) = punt[i];
 
     return data;
   }
@@ -441,14 +405,12 @@
   N_Vector
   IDA::ColToNVec (const ColumnVector &data, long int n)
   {
-    N_Vector v =
-      N_VNew_Serial (n);
+    N_Vector v = N_VNew_Serial (n);
 
-    realtype * punt;
-    punt = nv_data_s (v);
+    realtype *punt = nv_data_s (v);
 
     for (octave_idx_type i = 0; i < n; i++)
-      punt [i] = data (i);
+      punt[i] = data(i);
 
     return v;
   }
@@ -456,62 +418,44 @@
   void
   IDA::set_userdata (void)
   {
-    void * userdata = this;
+    void *userdata = this;
 
-    int flag =
-      IDASetUserData(mem, userdata);
-    if (flag != 0)
-      {
-        error ("User data not set");
-      }
+    if (IDASetUserData (mem, userdata) != 0)
+      error ("User data not set");
   }
 
   void
   IDA::initialize (void)
   {
-    num = y0.numel();
+    num = y0.numel ();
     mem = IDACreate ();
 
-    N_Vector yy =
-      ColToNVec(y0, num);
+    N_Vector yy = ColToNVec (y0, num);
 
-    N_Vector yyp =
-      ColToNVec(yp0, num);
+    N_Vector yyp = ColToNVec (yp0, num);
 
     IDA::set_userdata ();
 
-    int flag =
-      IDAInit (mem, IDA::resfun, t0, yy, yyp);
-    if (flag != 0)
-      {
-        error ("IDA not initialized");
-      }
+    if (IDAInit (mem, IDA::resfun, 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, num);
 
-    int flag =
-      IDASVtolerances (mem, reltol, abs_tol);
-    if (flag != 0)
-      {
-        error ("IDA: Tolerance not set");
-      }
+    if (IDASVtolerances (mem, reltol, abs_tol) != 0)
+      error ("IDA: Tolerance not set");
+
     N_VDestroy_Serial (abs_tol);
   }
 
   void
   IDA::set_tolerance (realtype abstol, realtype reltol)
   {
-    int flag =
-      IDASStolerances (mem, reltol, abstol);
-    if (flag != 0)
-      {
-        error ("IDA: Tolerance not set");
-      }
+    if (IDASStolerances (mem, reltol, abstol) != 0)
+      error ("IDA: Tolerance not set");
   }
 
   octave_value_list
@@ -526,22 +470,17 @@
     ColumnVector tout, yout (num), ypout (num), ysel (outputsel.numel ());
     ColumnVector ie, te, oldval, oldisterminal, olddir;
     Matrix output, ye;
-    int cont = 0, flag = 0, temp = 0;
+    int cont = 0, temp = 0;
     bool status = 0;
     std::string string = "";
     ColumnVector yold = y;
-    octave_value_list retval;
+
+    realtype tsol = tspan(0);
+    realtype tend = tspan(numt-1);
 
-    realtype tsol =
-      tspan (0);
-    realtype tend =
-      tspan (numt - 1);
+    N_Vector yyp = ColToNVec (yp, num);
 
-    N_Vector yyp =
-      ColToNVec (yp, num);
-
-    N_Vector yy =
-      ColToNVec (y, num);
+    N_Vector yy = ColToNVec (y, num);
 
     // Initialize OutputFcn
     if (haveoutputfcn)
@@ -558,7 +497,7 @@
       {
         // First output value
         tout.resize (numt);
-        tout (0) = tsol;
+        tout(0) = tsol;
         output.resize (numt, num);
 
         for (octave_idx_type i = 0; i < num; i++)
@@ -568,15 +507,13 @@
         for (octave_idx_type j = 1; j < numt && status == 0; j++)
           {
             // IDANORMAL already interpolates tspan(j)
-            flag = IDASolve (mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL);
-            if (flag != 0)
-              {
-                error ("IDASolve failed");
-              }
+
+            if (IDASolve (mem, tspan (j), &tsol, yy, yyp, IDA_NORMAL) != 0)
+              error ("IDASolve failed");
 
             yout = NVecToCol (yy, num);
             ypout = NVecToCol (yyp, num);
-            tout (j) = tsol;
+            tout(j) = tsol;
 
             for (octave_idx_type i = 0; i < num; i++)
               output.elem (j, i) = yout.elem (i);
@@ -586,9 +523,9 @@
                                        tend, outputsel, string);
 
             if (haveeventfunction)
-              status = IDA::event (event_fcn, te, ye, ie, tout (j), yout,
+              status = IDA::event (event_fcn, te, ye, ie, tout(j), yout,
                                    string, ypout, oldval, oldisterminal,
-                                   olddir, j, temp, tout (j - 1), yold);
+                                   olddir, j, temp, tout(j-1), yold);
 
             // If integration is stopped, return only the reached steps
             if (status == 1)
@@ -603,25 +540,21 @@
       {
         // First output value
         tout.resize (1);
-        tout (0) = tsol;
+        tout(0) = tsol;
         output.resize (1, num);
 
         for (octave_idx_type i = 0; i < num; i++)
           output.elem (0, i) = y.elem (i);
 
-        bool posdirection =
-          (tend > tsol);
+        bool posdirection = (tend > tsol);
 
         //main loop
         while (((posdirection == 1 && tsol < tend)
                || (posdirection == 0 && tsol > tend))
                && status == 0)
           {
-            flag = IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP);
-            if (flag != 0)
-              {
-                error ("IDASolve failed");
-              }
+            if (IDASolve (mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0)
+              error ("IDASolve failed");
 
             if (haverefine)
               status = IDA::interpolate (cont, output, tout, refine, tend,
@@ -635,34 +568,31 @@
             cont += 1;
             output.resize (cont + 1, num); // This may be not efficient
             tout.resize (cont + 1);
-            tout (cont) = tsol;
+            tout(cont) = tsol;
             yout = NVecToCol (yy, num);
 
             for (octave_idx_type i = 0; i < num; i++)
               output.elem (cont, i) = yout.elem (i);
 
-            if (haveoutputfcn && ! haverefine && tout (cont) < tend)
+            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,
+            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);
+                                   olddir, cont, temp, tout(cont-1), yold);
           }
+
         if (status == 0)
           {
             // Interpolate in tend
-            N_Vector dky =
-              N_VNew_Serial (num);
+            N_Vector dky = N_VNew_Serial (num);
 
-            flag = IDAGetDky (mem, tend, 0, dky);
-            if (flag != 0)
-              {
-                error ("IDA failed to interpolate y");
-              }
+            if (IDAGetDky (mem, tend, 0, dky) != 0)
+              error ("IDA failed to interpolate y");
 
-            tout (cont) = tend;
+            tout(cont) = tend;
             yout = NVecToCol (dky, num);
 
             for (octave_idx_type i = 0; i < num; i++)
@@ -678,25 +608,20 @@
                 if (haveeventfunction)
                   status = IDA::event (event_fcn, te, ye, ie, tend, yout,
                                        string, ypout, oldval, oldisterminal,
-                                       olddir, cont, temp, tout (cont - 1),
+                                       olddir, cont, temp, tout(cont-1),
                                        yold);
               }
 
             N_VDestroy_Serial (dky);
           }
+
         // Cleanup plotter
         status = IDA::outputfun (output_fcn, haveoutputsel, yout, tend, tend,
                                  outputsel, "done");
 
       }
 
-    retval (0) = tout;
-    retval (1) = output;
-    retval (2) = te;
-    retval (3) = ye;
-    retval (4) = ie;
-
-    return retval;
+    return ovl (tout, output, te, ye, ie);
   }
 
   bool
@@ -708,19 +633,15 @@
               int& temp, realtype told, ColumnVector& yold)
   {
     bool status = 0;
-    octave_value_list args (3);
-    octave_value_list output (3);
-    ColumnVector val, isterminal, dir;
-    args (0) = tsol;
-    args (1) = y;
-    args (2) = yp;
+
+    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")
       {
-        output = feval (event_fcn, args, 3);
+        octave_value_list output = feval (event_fcn, args, 3);
         oldval = output(0).vector_value ();
         oldisterminal = output(1).vector_value ();
         olddir = output(2).vector_value ();
@@ -728,10 +649,10 @@
     else if (flag == "")
       {
         ColumnVector index (0);
-        output = feval (event_fcn, args, 3);
-        val = output(0).vector_value ();
-        isterminal = output(1).vector_value ();
-        dir = output(2).vector_value ();
+        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++)
@@ -752,12 +673,12 @@
             ie.resize (1);
 
             // Linear interpolation
-            ie (0) = index (0);
-            te (0) = tsol - val (index (0)) * (tsol - told)
-              / (val (index (0)) - oldval (index (0)));
+            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));
+            ColumnVector ytemp
+              = y - ((tsol - te(0)) * (y - yold) / (tsol - told));
 
             for (octave_idx_type i = 0; i < num; i++)
               ye.elem (0, i) = ytemp.elem (i);
@@ -773,17 +694,16 @@
             for (octave_idx_type i = 0; i < index.numel (); i++)
               {
 
-                if (isterminal (index (i)) == 1)
+                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)));
+                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);
+                ColumnVector ytemp
+                  = y - (tsol - te (temp + i)) * (y - yold) / (tsol - told);
 
                 for (octave_idx_type j = 0; j < num; j++)
                   ye.elem (temp + i, j) = ytemp.elem (j);
@@ -817,53 +737,39 @@
     realtype h = 0, tcur = 0;
     bool status = 0;
 
-    N_Vector dky =
-      N_VNew_Serial (num);
+    N_Vector dky = N_VNew_Serial (num);
 
-    N_Vector dkyp =
-      N_VNew_Serial (num);
+    N_Vector dkyp = N_VNew_Serial (num);
 
     ColumnVector yout (num);
     ColumnVector ypout (num);
     std::string string = "";
 
-    int flag =
-      IDAGetLastStep (mem, &h);
+    if (IDAGetLastStep (mem, &h) != 0)
+      error ("IDA failed to return last step");
 
-    if (flag != 0)
-      {
-        error ("IDA failed to return last step");
-      }
-    flag = IDAGetCurrentTime (mem, &tcur);
-    if (flag != 0)
-      {
-        error ("IDA failed to return the current time");
-      }
+    if (IDAGetCurrentTime (mem, &tcur) != 0)
+      error ("IDA failed to return the current time");
 
-    realtype tin =
-      tcur - h;
+    realtype tin = tcur - h;
+
+    realtype step = h / refine;
 
-    realtype step =
-      h / refine;
-
-    for (octave_idx_type i = 1; i < refine && tin + step * i < tend
-           && status == 0; i++)
+    for (octave_idx_type i = 1;
+         i < refine && tin + step * i < tend && status == 0;
+         i++)
       {
-        flag = IDAGetDky (mem, tin + step*i, 0, dky);
-        if (flag != 0)
-          {
-            error ("IDA failed to interpolate y");
-          }
-        flag = IDAGetDky (mem, tin + step*i, 1, dkyp);
-        if (flag != 0)
-          {
-            error ("IDA failed to interpolate yp");
-          }
+        if (IDAGetDky (mem, tin + step*i, 0, dky) != 0)
+          error ("IDA failed to interpolate y");
+
+        if (IDAGetDky (mem, tin + step*i, 1, dkyp) != 0)
+          error ("IDA failed to interpolate yp");
+
         cont += 1;
         output.resize (cont + 1, num);
         tout.resize (cont + 1);
 
-        tout (cont) = tin + step * i;
+        tout(cont) = tin + step * i;
         yout = NVecToCol (dky, num);
         ypout = NVecToCol (dkyp, num);
 
@@ -872,15 +778,17 @@
 
         if (haveoutputfcn)
           status = IDA::outputfun (output_fcn, haveoutputsel, yout,
-                                   tout (cont), tend, outputsel, "");
+                                   tout(cont), tend, outputsel, "");
 
         if (haveeventfunction)
-          status = IDA::event (event_fcn, te, ye, ie, tout (cont),
+          status = IDA::event (event_fcn, te, ye, ie, tout(cont),
                                yout, string, ypout, oldval,
                                oldisterminal, olddir, cont, temp,
-                               tout (cont - 1), yold);
+                               tout(cont-1), yold);
       }
+
     N_VDestroy_Serial (dky);
+
     return status;
   }
 
@@ -891,39 +799,39 @@
                   std::string flag)
   {
     bool status = 0;
-    octave_value_list output, val;
+
+    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));
+          ysel(i) = yout(outputsel(i));
 
-        output (1) = ysel;
+        output(1) = ysel;
       }
     else
-      output (1) = yout;
-
-    output (2) = flag;
+      output(1) = yout;
 
     if (flag == "init")
       {
-        ColumnVector toutput (2);
-        toutput (0) = tsol;
-        toutput (1) = tend;
-        output (0) = toutput;
+        ColumnVector toutput(2);
+        toutput(0) = tsol;
+        toutput(1) = tend;
+        output(0) = toutput;
 
         feval (output_fcn, output, 0);
       }
     else if (flag == "")
       {
-        output (0) = tsol;
-        val = feval (output_fcn, output, 1);
+        output(0) = tsol;
+        octave_value_list val = feval (output_fcn, output, 1);
         status = val(0).bool_value ();
       }
     else
       {  // Cleanup plotter
-        output (0) = tend;
+        output(0) = tend;
         feval (output_fcn, output, 0);
       }
 
@@ -933,166 +841,111 @@
   void
   IDA::set_maxstep (realtype maxstep)
   {
-    int flag =
-      IDASetMaxStep (mem, maxstep);
-
-    if (flag != 0)
-      {
-        error ("IDA: Max Step not set");
-      }
+    if (IDASetMaxStep (mem, maxstep) != 0)
+      error ("IDA: Max Step not set");
   }
 
   void
   IDA::set_initialstep (realtype initialstep)
   {
-    int flag =
-      IDASetInitStep (mem, initialstep);
-
-    if (flag != 0)
-      {
-        error ("IDA: Initial Step not set");
-      }
+    if (IDASetInitStep (mem, initialstep) != 0)
+      error ("IDA: Initial Step not set");
   }
 
   void
   IDA::set_maxorder (int maxorder)
   {
-    int flag =
-      IDASetMaxOrd (mem, maxorder);
-
-    if (flag != 0)
-      {
-        error ("IDA: Max Order not set");
-      }
+    if (IDASetMaxOrd (mem, maxorder) != 0)
+      error ("IDA: Max Order not set");
   }
 
   void
   IDA::print_stat (void)
   {
     long int nsteps = 0, netfails = 0, nrevals = 0;
-    int flag =
-      IDAGetNumSteps(mem, &nsteps);
+
+    if (IDAGetNumSteps (mem, &nsteps) != 0)
+      error ("IDA failed to return the number of internal steps");
+
+    if (IDAGetNumErrTestFails (mem, &netfails) != 0)
+      error ("IDA failed to return the number of internal errors");
 
-    if (flag != 0)
-      {
-        error ("IDA failed to return the number of internal steps");
-      }
-    flag = IDAGetNumErrTestFails(mem, &netfails);
-    if (flag != 0)
-      {
-        error ("IDA failed to return the number of internal errors");
-      }
-    flag = IDAGetNumResEvals(mem, &nrevals);
-    if (flag != 0)
-      {
-        error ("IDA failed to return the number of residual evaluations");
-      }
+    if (IDAGetNumResEvals (mem, &nrevals) != 0)
+      error ("IDA failed to return the number of residual evaluations");
 
-    std::cout<<nsteps<<" successful steps\n";
-    std::cout<<netfails<<" failed attempts\n";
-    std::cout<<nrevals<<" function evaluations\n";
-    //std::cout<<<<" partial derivatives\n";
-    //std::cout<<<<" LU decompositions\n";
-    //std::cout<<<<" solutions of linear systems\n";
+    std::cout << nsteps << " successful steps\n";
+    std::cout << netfails << " failed attempts\n";
+    std::cout << nrevals << " function evaluations\n";
+    // std::cout << " partial derivatives\n";
+    // std::cout << " LU decompositions\n";
+    // std::cout << " solutions of linear systems\n";
   }
 
   ColumnVector
   ida_user_function (const ColumnVector& x, const ColumnVector& xdot,
                      double t, octave_function *ida_fc)
   {
-    ColumnVector retval;
-    octave_value_list args;
-
-    args(2) = xdot;
-    args(1) = x;
-    args(0) = t;
-
     octave_value_list tmp;
 
     try
       {
-        tmp = ida_fc -> do_multi_index_op (1, args);
+        tmp = ida_fc->do_multi_index_op (1, ovl (xdot, x, t));
       }
     catch (octave::execution_exception& e)
       {
         err_user_supplied_eval (e, "__ode15__");
       }
 
-    retval = tmp(0).vector_value ();
-
-    return retval;
+    return tmp(0).vector_value ();
   }
 
   Matrix
   ida_dense_jac (const ColumnVector& x, const ColumnVector& xdot,
                  double t, double cj, octave_function *ida_jc)
   {
-    Matrix retval;
-    octave_value_list newargs (3);
-
-    newargs (0) = t;
-    newargs (1) = x;
-    newargs (2) = xdot;
-
     octave_value_list tmp;
 
     try
       {
-        tmp = ida_jc -> do_multi_index_op (2, newargs);
+        tmp = ida_jc->do_multi_index_op (2, ovl (t, x, xdot));
       }
     catch (octave::execution_exception& e)
       {
         err_user_supplied_eval (e, "__ode15__");
       }
 
-    retval = tmp(0).matrix_value () + cj * tmp(1).matrix_value ();
-
-    return retval;
+    return tmp(0).matrix_value () + cj * tmp(1).matrix_value ();
   }
 
   SparseMatrix
   ida_sparse_jac (const ColumnVector& x, const ColumnVector& xdot,
                   double t, double cj, octave_function *ida_jc)
   {
-    SparseMatrix retval;
-    octave_value_list newargs (3);
-
-    newargs (0) = t;
-    newargs (1) = x;
-    newargs (2) = xdot;
-
     octave_value_list tmp;
 
     try
       {
-        tmp = ida_jc -> do_multi_index_op (2, newargs);
+        tmp = ida_jc->do_multi_index_op (2, ovl (t, x, xdot));
       }
     catch (octave::execution_exception& e)
       {
         err_user_supplied_eval (e, "__ode15__");
       }
 
-    retval = tmp(0).sparse_matrix_value () +
-      cj * tmp(1).sparse_matrix_value ();
-
-    return retval;
+    return tmp(0).sparse_matrix_value () + cj * tmp(1).sparse_matrix_value ();
   }
 
   Matrix
   ida_dense_cell_jac (Matrix *dfdy, Matrix *dfdyp, double cj)
   {
-    Matrix retval;
-    retval = (*dfdy) + cj * (*dfdyp);
-    return retval;
+    return (*dfdy) + cj * (*dfdyp);
   }
 
   SparseMatrix
   ida_sparse_cell_jac (SparseMatrix *spdfdy, SparseMatrix *spdfdyp,
                        double cj)
   {
-    SparseMatrix retval;
-    retval = (*spdfdy) + cj * (*spdfdyp);
-    return retval;
+    return (*spdfdy) + cj * (*spdfdyp);
   }
 
   octave_value_list
@@ -1110,14 +963,11 @@
     IDA dae (t0, y0, yp0, ida_fcn, ida_user_function);
 
     // Set Jacobian
-    bool havejac =
-      options.getfield ("havejac").bool_value ();
+    bool havejac = options.getfield ("havejac").bool_value ();
 
-    bool havejacsparse =
-      options.getfield ("havejacsparse").bool_value ();
+    bool havejacsparse = options.getfield ("havejacsparse").bool_value ();
 
-    bool havejacfun =
-      options.getfield ("havejacfun").bool_value ();
+    bool havejacfun = options.getfield ("havejacfun").bool_value ();
 
     Matrix ida_dfdy, ida_dfdyp, *dfdy, *dfdyp;
     SparseMatrix ida_spdfdy, ida_spdfdyp, *spdfdy, *spdfdyp;
@@ -1129,6 +979,7 @@
         if (havejacfun)
           {
             ida_jac = options.getfield ("Jacobian").function_value ();
+
             if (havejacsparse)
               dae.set_jacobian (ida_jac, ida_sparse_jac);
             else
@@ -1161,68 +1012,58 @@
     dae.initialize ();
 
     // Set tolerances
-    realtype rel_tol =
-      options.getfield("RelTol").double_value ();
+    realtype rel_tol = options.getfield("RelTol").double_value ();
 
-    bool haveabstolvec =
-      options.getfield ("haveabstolvec").bool_value ();
+    bool haveabstolvec = options.getfield ("haveabstolvec").bool_value ();
 
     if (haveabstolvec)
       {
-        ColumnVector abs_tol =
-          options.getfield("AbsTol").vector_value ();
+        ColumnVector abs_tol = options.getfield("AbsTol").vector_value ();
 
         dae.set_tolerance (abs_tol, rel_tol);
       }
     else
       {
-        realtype abs_tol =
-          options.getfield("AbsTol").double_value ();
+        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 ();
+    realtype maxstep = options.getfield("MaxStep").double_value ();
 
     dae.set_maxstep (maxstep);
 
     //Set initial step
-    if (!(options.getfield("InitialStep").is_empty ()))
+    if (! options.getfield("InitialStep").is_empty ())
       {
-        realtype initialstep =
-          options.getfield("InitialStep").double_value ();
+        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 ();
+    int maxorder = options.getfield("MaxOrder").int_value ();
 
     dae.set_maxorder (maxorder);
 
     //Set Refine
-    const int refine =
-      options.getfield("Refine").int_value ();
+    const int refine = options.getfield("Refine").int_value ();
 
-    bool haverefine =
-      (refine > 1);
+    bool haverefine = (refine > 1);
 
     octave_function *output_fcn = 0;
     ColumnVector outputsel;
 
     // OutputFcn
-    bool haveoutputfunction =
-      options.getfield("haveoutputfunction").bool_value ();
+    bool haveoutputfunction
+      = options.getfield("haveoutputfunction").bool_value ();
 
     if (haveoutputfunction)
       output_fcn = options.getfield("OutputFcn").function_value ();
 
     // OutputSel
-    bool haveoutputsel =
-      options.getfield("haveoutputselection").bool_value ();
+    bool haveoutputsel = options.getfield("haveoutputselection").bool_value ();
 
     if (haveoutputsel)
       outputsel = options.getfield("OutputSel").vector_value ();
@@ -1230,8 +1071,8 @@
     octave_function *event_fcn = 0;
 
     // Events
-    bool haveeventfunction =
-      options.getfield("haveeventfunction").bool_value ();
+    bool haveeventfunction
+      = options.getfield("haveeventfunction").bool_value ();
 
     if (haveeventfunction)
       event_fcn = options.getfield("Events").function_value ();
@@ -1246,8 +1087,7 @@
                             haveeventfunction, event_fcn);
 
     // Statistics
-    bool havestats =
-      options.getfield("havestats").bool_value ();
+    bool havestats = options.getfield("havestats").bool_value ();
 
     if (havestats)
       dae.print_stat ();
@@ -1278,32 +1118,30 @@
 
   octave_value f_arg = args(0);
 
-  if (f_arg.is_function_handle ())
-    ida_fcn = f_arg.function_value ();
-  else
+  if (! f_arg.is_function_handle ())
     error ("__ode15__: odefun must be a function handle");
 
+  ida_fcn = f_arg.function_value ();
+
   // Check input tspan
-  ColumnVector tspan =
-    args(1).xvector_value ("__ode15__: TRANGE must be a vector of numbers");
+  ColumnVector tspan
+    = args(1).xvector_value ("__ode15__: TRANGE must be a vector of numbers");
 
-  int numt =
-    tspan.numel ();
+  int numt = tspan.numel ();
 
-  realtype t0 =
-    tspan (0);
+  realtype t0 = tspan (0);
 
   if (numt < 2)
     error ("__ode15__: TRANGE must contain at least 2 elements");
-  else if (!(tspan.is_sorted ()) || (tspan(0) == tspan(numt - 1)))
+  else if (! tspan.is_sorted () || 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 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");
+  ColumnVector yp0
+    = args(3).xvector_value ("__ode15__: initial state yp0 must be a vector");
 
 
   if (y0.numel () != yp0.numel ())
@@ -1315,13 +1153,11 @@
   if (! args(4).is_map ())
     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");
+  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);
+  return octave::do_ode15 (ida_fcn, tspan, numt, t0, y0, yp0, options);
 
 
 #else