Mercurial > octave

changeset 28716:9beec32ba3d6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Overhaul usage of integer types in __ode15__.cc (bug #58795). * __ode15__.cc: Use octave_f77_int_type or octave_idx_type where appropriate and check for potentially narrowing conversions.
author Markus Mützel <markus.muetzel@gmx.de>
date Thu, 03 Sep 2020 13:36:05 +0200
parents 0fcbb0faf7de
children 1ffe315b1052
files libinterp/dldfcn/__ode15__.cc
diffstat 1 files changed, 50 insertions(+), 39 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/dldfcn/__ode15__.cc	Thu Sep 10 14:59:10 2020 -0700
+++ b/libinterp/dldfcn/__ode15__.cc	Thu Sep 03 13:36:05 2020 +0200
@@ -30,6 +30,7 @@
 #include "dColVector.h"
 #include "dMatrix.h"
 #include "dSparse.h"
+#include "f77-fcn.h"
 
 #include "Cell.h"
 #include "defun-dld.h"
@@ -166,8 +167,8 @@
       : 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_spdfdyp (nullptr), m_fun (nullptr), m_jacfun (nullptr),
+        m_jacspfun (nullptr), m_jacdcell (nullptr), m_jacspcell (nullptr),
         m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr)
     { }
 
@@ -177,8 +178,8 @@
       : 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_spdfdyp (nullptr), m_fun (daefun), m_jacfun (nullptr),
+        m_jacspfun (nullptr), m_jacdcell (nullptr), m_jacspcell (nullptr),
         m_sunJacMatrix (nullptr), m_sunLinearSolver (nullptr)
     { }
 
@@ -245,9 +246,9 @@
 
     void initialize (void);
 
-    static ColumnVector NVecToCol (N_Vector& v, long int n);
+    static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n);
 
-    static N_Vector ColToNVec (const ColumnVector& data, long int n);
+    static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
 
     void
     set_up (const ColumnVector& y);
@@ -266,9 +267,9 @@
     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)
+    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);
@@ -291,8 +292,8 @@
     }
 
     void
-    jacsparse_impl (realtype t, realtype cj, N_Vector& yy,
-                    N_Vector& yyp, SUNMatrix& Jac);
+    jacsparse_impl (realtype t, realtype cj,
+                    N_Vector& yy, N_Vector& yyp, SUNMatrix& Jac);
 #  endif
 
     void set_maxstep (realtype maxstep);
@@ -300,14 +301,14 @@
     void set_initialstep (realtype initialstep);
 
     bool
-    interpolate (int& cont, Matrix& output, ColumnVector& tout,
+    interpolate (octave_idx_type& 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);
+                 octave_idx_type& temp, ColumnVector& yold);
 
     bool
     outputfun (const octave_value& output_fcn, bool haveoutputsel,
@@ -321,12 +322,13 @@
            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);
+           octave_idx_type cont, octave_idx_type& temp, realtype told,
+           ColumnVector& yold);
 
     void set_maxorder (int maxorder);
 
     octave_value_list
-    integrate (const int numt, const ColumnVector& tt,
+    integrate (const octave_idx_type numt, const ColumnVector& tt,
                const ColumnVector& y0, const ColumnVector& yp0,
                const int refine, bool haverefine, bool haveoutputfcn,
                const octave_value& output_fcn, bool haveoutputsel,
@@ -344,7 +346,7 @@
     bool m_havejacfun;
     bool m_havejacsparse;
     void *m_mem;
-    int m_num;
+    octave_f77_int_type m_num;
     octave_value m_ida_fun;
     octave_value m_ida_jac;
     Matrix *m_dfdy;
@@ -396,6 +398,9 @@
         // 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.
+        // FIXME: m_num*m_num might be larger than the largest
+        // octave_f77_int_type (integer overflow).  Consider using a save
+        // calculation method.
         m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, m_num*m_num, CSC_MAT);
         if (! m_sunJacMatrix)
           error ("Unable to create sparse Jacobian for Sundials");
@@ -410,7 +415,9 @@
         IDASetJacFn (m_mem, IDA::jacsparse);
 
 #  else
-        error ("SUNDIALS SUNLINSOL KLU is not available in this version of Octave");
+        error ("SUNDIALS SUNLINSOL KLU was unavailable or disabled when "
+               "Octave was built");
+
 #  endif
 
       }
@@ -439,7 +446,7 @@
                       N_Vector& yy, N_Vector& yyp, SUNMatrix& JJ)
 
   {
-    long int Neq = NV_LENGTH_S (yy);
+    octave_f77_int_type Neq = NV_LENGTH_S (yy);
 
     ColumnVector y = NVecToCol (yy, Neq);
 
@@ -452,8 +459,9 @@
     else
       jac = (*m_jacdcell) (m_dfdy, m_dfdyp, cj);
 
+    octave_f77_int_type num_jac = to_f77_int (jac.numel ());
     std::copy (jac.fortran_vec (),
-               jac.fortran_vec () + jac.numel (),
+               jac.fortran_vec () + num_jac,
                SUNDenseMatrix_Data (JJ));
   }
 
@@ -475,41 +483,43 @@
       jac = (*m_jacspcell) (m_spdfdy, m_spdfdyp, cj);
 
     SUNMatZero_Sparse (Jac);
+    // We have to use "sunindextype *" here but still need to check that
+    // conversion of each element to "octave_f77_int_type" is save.
     sunindextype *colptrs = SUNSparseMatrix_IndexPointers (Jac);
     sunindextype *rowvals = SUNSparseMatrix_IndexValues (Jac);
 
-    for (int i = 0; i < m_num + 1; i++)
-      colptrs[i] = jac.cidx (i);
+    for (octave_f77_int_type i = 0; i < m_num + 1; i++)
+      colptrs[i] = to_f77_int (jac.cidx (i));
 
     double *d = SUNSparseMatrix_Data (Jac);
-    for (int i = 0; i < jac.nnz (); i++)
+    for (octave_f77_int_type i = 0; i < to_f77_int (jac.nnz ()); i++)
       {
-        rowvals[i] = jac.ridx (i);
+        rowvals[i] = to_f77_int (jac.ridx (i));
         d[i] = jac.data (i);
       }
   }
 #  endif
 
   ColumnVector
-  IDA::NVecToCol (N_Vector& v, long int n)
+  IDA::NVecToCol (N_Vector& v, octave_f77_int_type n)
   {
     ColumnVector data (n);
     realtype *punt = nv_data_s (v);
 
-    for (octave_idx_type i = 0; i < n; i++)
+    for (octave_f77_int_type i = 0; i < n; i++)
       data(i) = punt[i];
 
     return data;
   }
 
   N_Vector
-  IDA::ColToNVec (const ColumnVector& data, long int n)
+  IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n)
   {
     N_Vector v = N_VNew_Serial (n);
 
     realtype *punt = nv_data_s (v);
 
-    for (octave_idx_type i = 0; i < n; i++)
+    for (octave_f77_int_type i = 0; i < n; i++)
       punt[i] = data(i);
 
     return v;
@@ -527,7 +537,7 @@
   void
   IDA::initialize (void)
   {
-    m_num = m_y0.numel ();
+    m_num = to_f77_int (m_y0.numel ());
     m_mem = IDACreate ();
 
     N_Vector yy = ColToNVec (m_y0, m_num);
@@ -559,7 +569,7 @@
   }
 
   octave_value_list
-  IDA::integrate (const int numt, const ColumnVector& tspan,
+  IDA::integrate (const octave_idx_type numt, const ColumnVector& tspan,
                   const ColumnVector& y, const ColumnVector& yp,
                   const int refine, bool haverefine, bool haveoutputfcn,
                   const octave_value& output_fcn, bool haveoutputsel,
@@ -570,7 +580,7 @@
     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;
+    octave_idx_type cont = 0, temp = 0;
     bool status = 0;
     std::string string = "";
     ColumnVector yold = y;
@@ -726,11 +736,12 @@
 
   bool
   IDA::event (const octave_value& event_fcn,
-              ColumnVector& te, Matrix& ye, ColumnVector& ie,
-              realtype tsol, const ColumnVector& y, const std::string& flag,
+              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)
+              ColumnVector& oldisterminal, ColumnVector& olddir,
+              octave_idx_type cont, octave_idx_type& temp, realtype told,
+              ColumnVector& yold)
   {
     bool status = 0;
 
@@ -825,14 +836,14 @@
   }
 
   bool
-  IDA::interpolate (int& cont, Matrix& output, ColumnVector& tout,
+  IDA::interpolate (octave_idx_type& 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)
+                    octave_idx_type& temp, ColumnVector& yold)
   {
     realtype h = 0, tcur = 0;
     bool status = 0;
@@ -1052,7 +1063,7 @@
   octave_value_list
   do_ode15 (const octave_value& ida_fcn,
             const ColumnVector& tspan,
-            const int numt,
+            const octave_idx_type numt,
             const realtype t0,
             const ColumnVector& y0,
             const ColumnVector& yp0,
@@ -1206,7 +1217,7 @@
 #if defined (HAVE_SUNDIALS)
 
   // Check number of parameters
-  int nargin = args.length ();
+  octave_idx_type nargin = args.length ();
 
   if (nargin != 5)
     print_usage ();
@@ -1221,7 +1232,7 @@
   ColumnVector tspan
     = args(1).xvector_value ("__ode15__: TRANGE must be a vector of numbers");
 
-  int numt = tspan.numel ();
+  octave_idx_type numt = tspan.numel ();
 
   realtype t0 = tspan (0);