lazy contiguous subrange indexing

[Jaroslav Hajek-2]

Hello,

please consider the attached two patches.

Summary (1st patch):
If a dense indexing expression can be reduced to a contiguous subrange
(like x(2:end-1), A(:,:,3:10) etc.), a shared (i.e. cheap) copy is
made and the actual extraction is postponed until (if ever) the value
is written to. The obvious benefit is greater performance of various
operations with indexed expressions, such as y = x(1:n-1) + x(2:n),
which will avoid making two temporary copies. This is done by
modifying Array<T> to allow it work with contiguous subranges of its
(possibly shared) data. There is an obvious risk of increased memory
consumption, because a big block of memory may be kept allocated while
only part of it is actually used.
This is resolved by checking prior to storing a value to a "permanent"
location (same where checks for null matrices are done), and possibly
economizing if an orphaned block is detected. The detection (see
Array<T>::maybe_economize) is fairly simple and could be improved at
the cost of performance of the operation. Since most indexing
expressions do not live longer than their parent object I suppose this
will seldom be an issue, but it is still easily possible to construct
an example where memory is wastefully kept allocated.
However, the issue is easy to avoid if you know what you do, so I
think it's OK. Comments are welcome.

To get an idea of the performance speed up, here's a simplistic test
(set n to suitable value):

n = 2e7;
x = 1/n * [1:n]; y = x.^2;
# compute 2nd-order differences
tic; d2y = diff (y, 2); toc
# compute integral
tic; int = trapz (x, y); toc

with a recent tip, I get:

Elapsed time is 0.697262 seconds.
Elapsed time is 0.960276 seconds.

with the new patch, I got:

Elapsed time is 0.213054 seconds.
Elapsed time is 0.497459 seconds.

which means speedups by 228% and 93%.
P.S. I picked some benefitted functions because measuring directly the
performance speed-up of the affected indexing expressions is
meaningless.


Summary (2nd patch):
This is a long-postponed cleanup of Array<T>, DiagArray2<T> and
PermMatrix, to allow the private members of Array<T> to become really
private,
and reimplement DiagArray and PermMatrix without the dangerous abuse
of Array<T>::dimensions (that has caused me quite a pain with the
previous patch).
This patch should probably be applied in any case, I'm putting it here
because right now it's created on top of the 1st patch and I've only
ran make check on the result.

comments? OK to push?

cheers to Octave

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[lazy.diff]

---

# HG changeset patch
# User Jaroslav Hajek <highegg@...>
# Date 1231945605 -3600
# Node ID f3a54bae02c8e81cd79a8ceac45065f8af07f7fc
# Parent  0e0bd07e6ae2df29a519e2c26a5561a5f5d4a0be
implement non-copying contiguous range indexing

diff --git a/liboctave/Array.cc b/liboctave/Array.cc
--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -3,7 +3,7 @@
 
 Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004,
               2005, 2006, 2007 John W. Eaton
-Copyright (C) 2008 Jaroslav Hajek <highegg@...>
+Copyright (C) 2008, 2009 Jaroslav Hajek <highegg@...>
 
 This file is part of Octave.
 
@@ -47,7 +47,8 @@
 
 template <class T>
 Array<T>::Array (const Array<T>& a, const dim_vector& dv)
-  : rep (a.rep), dimensions (dv)
+  : rep (a.rep), dimensions (dv),
+    slice_data (a.slice_data), slice_len (a.slice_len)
 {
   rep->count++;
 
@@ -76,6 +77,8 @@
       rep->count++;
 
       dimensions = a.dimensions;
+      slice_data = a.slice_data;
+      slice_len = a.slice_len;
     }
 
   return *this;
@@ -680,6 +683,11 @@
   template <class T>
   void fill (const T& val, T *dest) const { do_fill (val, dest, top); }
 
+  bool is_cont_range (octave_idx_type& l,
+                            octave_idx_type& u) const
+    {
+      return top == 0 && idx[0].is_cont_range (dim[0], l, u);
+    }
 };
 
 // Helper class for multi-d recursive resizing
@@ -758,9 +766,7 @@
   if (i.is_colon ())
     {
       // A(:) produces a shallow copy as a column vector.
-      retval.dimensions = dim_vector (n, 1);
-      rep->count++;
-      retval.rep = rep;
+      retval = Array<T> (*this, dim_vector (n, 1));
     }
   else if (i.extent (n) != n)
     {
@@ -797,12 +803,19 @@
             rd = dim_vector (1, il);
         }
 
-      // Don't use resize here to avoid useless initialization for POD
-      // types.
-      retval = Array<T> (rd);
+      octave_idx_type l, u;
+      if (il != 0 && i.is_cont_range (n, l, u))
+        // If suitable, produce a shallow slice.
+        retval = Array<T> (*this, rd, l, u);
+      else
+        {
+          // Don't use resize here to avoid useless initialization for POD
+          // types.
+          retval = Array<T> (rd);
 
-      if (il != 0)
-        i.index (data (), n, retval.fortran_vec ());
+          if (il != 0)
+            i.index (data (), n, retval.fortran_vec ());
+        }
     }
 
   return retval;
@@ -830,18 +843,30 @@
     {
       octave_idx_type n = numel (), il = i.length (r), jl = j.length (c);
 
-      // Don't use resize here to avoid useless initialization for POD types.
-      retval = Array<T> (dim_vector (il, jl));
-
       idx_vector ii (i);
 
-      const T* src = data ();
-      T *dest = retval.fortran_vec ();
+      if (ii.maybe_reduce (r, j, c))
+        {
+          octave_idx_type l, u;
+          if (ii.length () > 0 && ii.is_cont_range (n, l, u))
+            // If suitable, produce a shallow slice.
+            retval = Array<T> (*this, dim_vector (il, jl), l, u);
+          else
+            {
+              // Don't use resize here to avoid useless initialization for POD types.
+              retval = Array<T> (dim_vector (il, jl));
 
-      if (ii.maybe_reduce (r, j, c))
-        ii.index (src, n, dest);
+              ii.index (data (), n, retval.fortran_vec ());
+            }
+        }
       else
         {
+          // Don't use resize here to avoid useless initialization for POD types.
+          retval = Array<T> (dim_vector (il, jl));
+
+          const T* src = data ();
+          T *dest = retval.fortran_vec ();
+
           for (octave_idx_type k = 0; k < jl; k++)
             dest += i.index (src + r * j.xelem (k), r, dest);
         }
@@ -898,14 +923,21 @@
           for (int i = 0; i < ial; i++) rdv(i) = ia(i).length (dv(i));
           rdv.chop_trailing_singletons ();
 
-          // Don't use resize here to avoid useless initialization for POD types.
-          retval = Array<T> (rdv);
-
           // Prepare for recursive indexing
           rec_index_helper rh (dv, ia);
 
-          // Do it.
-          rh.index (data (), retval.fortran_vec ());
+          octave_idx_type l, u;
+          if (rh.is_cont_range (l, u))
+            // If suitable, produce a shallow slice.
+            retval = Array<T> (*this, rdv, l, u);
+          else
+            {
+              // Don't use resize here to avoid useless initialization for POD types.
+              retval = Array<T> (rdv);
+
+              // Do it.
+              rh.index (data (), retval.fortran_vec ());
+            }
         }
     }
 
@@ -1842,7 +1874,7 @@
 {
   make_unique ();
 
-  return rep->data;
+  return slice_data;
 }
 
 template <class T>
@@ -2168,10 +2200,12 @@
 void
 Array<T>::print_info (std::ostream& os, const std::string& prefix) const
 {
-  os << prefix << "rep address: " << rep << "\n"
-     << prefix << "rep->len:    " << rep->len << "\n"
-     << prefix << "rep->data:   " << static_cast<void *> (rep->data) << "\n"
-     << prefix << "rep->count:  " << rep->count << "\n";
+  os << prefix << "rep address: " << rep << '\n'
+     << prefix << "rep->len:    " << rep->len << '\n'
+     << prefix << "rep->data:   " << static_cast<void *> (rep->data) << '\n'
+     << prefix << "rep->count:  " << rep->count << '\n'
+     << prefix << "slice_data:  " << static_cast<void *> (slice_data) << '\n'
+     << prefix << "slice_len:   " << slice_len << '\n';
 
   // 2D info:
   //
diff --git a/liboctave/Array.h b/liboctave/Array.h
--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -3,7 +3,7 @@
 
 Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003,
               2004, 2005, 2006, 2007 John W. Eaton
-Copyright (C) 2008 Jaroslav Hajek <highegg@...>
+Copyright (C) 2008, 2009 Jaroslav Hajek <highegg@...>
 
 This file is part of Octave.
 
@@ -30,6 +30,7 @@
 #include <cstddef>
 
 #include <iostream>
+#include <algorithm>
 
 #include "dim-vector.h"
 #include "idx-vector.h"
@@ -58,7 +59,12 @@
     octave_idx_type len;
     int count;
 
-    ArrayRep (T *d, octave_idx_type l) : data (d), len (l), count (1) { }
+    ArrayRep (T *d, octave_idx_type l, bool copy = false)
+      : data (copy ? new T [l] : d), len (l), count (1)
+        {
+          if (copy)
+            std::copy (d, d + l, data);
+        }
 
     ArrayRep (void) : data (0), len (0), count (1) { }
 
@@ -67,29 +73,18 @@
     explicit ArrayRep (octave_idx_type n, const T& val)
       : data (new T [n]), len (n), count (1)
       {
- fill (val);
+        std::fill (data, data + n, val);
       }
 
     ArrayRep (const ArrayRep& a)
       : data (new T [a.len]), len (a.len), count (1)
       {
-        for (octave_idx_type i = 0; i < len; i++)
-  data[i] = a.data[i];
+        std::copy (a.data, a.data + a.len, data);
       }
 
     ~ArrayRep (void) { delete [] data; }
 
     octave_idx_type length (void) const { return len; }
-
-    void fill (const T& val)
-      {
- for (octave_idx_type i = 0; i < len; i++)
-  data[i] = val;
-      }
-
-    T& elem (octave_idx_type n) { return data[n]; }
-
-    T elem (octave_idx_type n) const { return data[n]; }
 
   private:
 
@@ -110,19 +105,29 @@
       if (rep->count > 1)
  {
   --rep->count;
-  rep = new ArrayRep (*rep);
+  rep = new ArrayRep (slice_data, slice_len, true);
+          slice_data = rep->data;
  }
-    }
+      else if (slice_len != rep->len)
+        {
+          // Possibly economize here.
+          ArrayRep *new_rep = new ArrayRep (slice_data, slice_len, true);
+          delete rep;
+          rep = new_rep;
+          slice_data = rep->data;
+        }
+      }
 
   void make_unique (const T& val)
     {
       if (rep->count > 1)
  {
   --rep->count;
-  rep = new ArrayRep (rep->length (), val);
+  rep = new ArrayRep (slice_len, val);
+          slice_data = rep->data;
  }
       else
- rep->fill (val);
+        std::fill (slice_data, slice_data + slice_len, val);
     }
 
   typedef T element_type;
@@ -136,12 +141,33 @@
 
 protected:
 
+  T* slice_data;
+  octave_idx_type slice_len;
+
   Array (T *d, octave_idx_type n)
-    : rep (new typename Array<T>::ArrayRep (d, n)), dimensions (n) { }
+    : rep (new typename Array<T>::ArrayRep (d, n)), dimensions (n)
+    {
+      slice_data = rep->data;
+      slice_len = rep->len;
+    }
 
   Array (T *d, const dim_vector& dv)
     : rep (new typename Array<T>::ArrayRep (d, get_size (dv))),
-      dimensions (dv) { }
+      dimensions (dv)
+    {
+      slice_data = rep->data;
+      slice_len = rep->len;
+    }
+
+  // slice constructor
+  Array (const Array<T>& a, const dim_vector& dv,
+         octave_idx_type l, octave_idx_type u)
+    : rep(a.rep), dimensions (dv)
+    {
+      rep->count++;
+      slice_data = a.slice_data + l;
+      slice_len = std::min (u, a.slice_len) - l;
+    }
 
 private:
 
@@ -168,14 +194,25 @@
 public:
 
   Array (void)
-    : rep (nil_rep ()), dimensions () { rep->count++; }
+    : rep (nil_rep ()), dimensions ()
+    {
+      rep->count++;
+      slice_data = rep->data;
+      slice_len = rep->len;
+    }
 
   explicit Array (octave_idx_type n)
-    : rep (new typename Array<T>::ArrayRep (n)), dimensions (n) { }
+    : rep (new typename Array<T>::ArrayRep (n)), dimensions (n)
+    {
+      slice_data = rep->data;
+      slice_len = rep->len;
+    }
 
   explicit Array (octave_idx_type n, const T& val)
     : rep (new typename Array<T>::ArrayRep (n)), dimensions (n)
     {
+      slice_data = rep->data;
+      slice_len = rep->len;
       fill (val);
     }
 
@@ -185,6 +222,8 @@
     : rep (new typename Array<T>::ArrayRep (coerce (a.data (), a.length ()), a.length ())),
       dimensions (a.dimensions)
     {
+      slice_data = rep->data;
+      slice_len = rep->len;
     }
 
   // No type conversion case.
@@ -192,18 +231,26 @@
     : rep (a.rep), dimensions (a.dimensions)
     {
       rep->count++;
+      slice_data = a.slice_data;
+      slice_len = a.slice_len;
     }
 
 public:
 
   Array (const dim_vector& dv)
     : rep (new typename Array<T>::ArrayRep (get_size (dv))),
-      dimensions (dv) { }
+      dimensions (dv)
+    {
+      slice_data = rep->data;
+      slice_len = rep->len;
+    }
 
   Array (const dim_vector& dv, const T& val)
     : rep (new typename Array<T>::ArrayRep (get_size (dv))),
       dimensions (dv)
     {
+      slice_data = rep->data;
+      slice_len = rep->len;
       fill (val);
     }
 
@@ -215,7 +262,7 @@
 
   void fill (const T& val) { make_unique (val); }
 
-  octave_idx_type capacity (void) const { return rep->length (); }
+  octave_idx_type capacity (void) const { return slice_len; }
   octave_idx_type length (void) const { return capacity (); }
   octave_idx_type nelem (void) const { return capacity (); }
   octave_idx_type numel (void) const { return nelem (); }
@@ -258,8 +305,8 @@
 
   // No checking, even for multiple references, ever.
 
-  T& xelem (octave_idx_type n) { return rep->elem (n); }
-  T xelem (octave_idx_type n) const { return rep->elem (n); }
+  T& xelem (octave_idx_type n) { return slice_data [n]; }
+  T xelem (octave_idx_type n) const { return slice_data [n]; }
 
   T& xelem (octave_idx_type i, octave_idx_type j) { return xelem (dim1()*j+i); }
   T xelem (octave_idx_type i, octave_idx_type j) const { return xelem (dim1()*j+i); }
@@ -279,7 +326,7 @@
 
   T& checkelem (octave_idx_type n)
     {
-      if (n < 0 || n >= rep->length ())
+      if (n < 0 || n >= slice_len)
  return range_error ("T& Array<T>::checkelem", n);
       else
  {
@@ -341,7 +388,7 @@
 
   T checkelem (octave_idx_type n) const
     {
-      if (n < 0 || n >= rep->length ())
+      if (n < 0 || n >= slice_len)
  return range_error ("T Array<T>::checkelem", n);
       else
  return xelem (n);
@@ -407,7 +454,7 @@
   Array<T> transpose (void) const;
   Array<T> hermitian (T (*fcn) (const T&) = 0) const;
 
-  const T *data (void) const { return rep->data; }
+  const T *data (void) const { return slice_data; }
 
   const T *fortran_vec (void) const { return data (); }
 
@@ -513,6 +560,17 @@
 
   Array<T>& insert (const Array<T>& a, const Array<octave_idx_type>& idx);
 
+  void maybe_economize (void)
+    {
+      if (rep->count == 1 && slice_len != rep->len)
+        {
+          ArrayRep *new_rep = new ArrayRep (slice_data, slice_len, true);
+          delete rep;
+          rep = new_rep;
+          slice_data = rep->data;
+        }
+    }
+
   void print_info (std::ostream& os, const std::string& prefix) const;
 
   // Unsafe.  This function exists to support the MEX interface.
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -0,0 +1,14 @@
[phone]  Jaroslav Hajek  <highegg@...>
+
+ * Array.cc, Array.h (all Array<T> constructors): Handle slice_data and
+ slice_len.
+ (Array<T>::Array<T> (const Array<T>&, const dim_vector&,
+ octave_idx_type, octave_idx_type)): New constructor.
+ (Array<T>::index): Use shallow copy when index reduces to a contiguous
+ range.
+ (Array<T>::make_unique): Rewrite.
+ (Array<T>::ArrayRep): Delete redundant methods.
+ (rec_index_helper::is_cont_range): New method.
+ (Array<T>::maybe_economize): New method.
+ * DiagArray2.cc (DiagArray2<T>::resize): Fix the mess.
+
diff --git a/liboctave/DiagArray2.cc b/liboctave/DiagArray2.cc
--- a/liboctave/DiagArray2.cc
+++ b/liboctave/DiagArray2.cc
@@ -107,6 +107,7 @@
   if (r == dim1 () && c == dim2 ())
     return;
 
+  // FIXME: this is a mess. DiagArray2 really needs a rewrite.
   typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
   const T *old_data = this->data ();
   octave_idx_type old_len = this->length ();
@@ -114,6 +115,8 @@
   octave_idx_type new_len = r < c ? r : c;
 
   Array<T>::rep = new typename Array<T>::ArrayRep (new_len);
+  Array<T>::slice_data = Array<T>::rep->data;
+  Array<T>::slice_len = Array<T>::rep->len;
 
   this->dimensions = dim_vector (r, c);
 
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -0,0 +1,15 @@
[phone]  Jaroslav Hajek  <highegg@...>
+
+ * ov.cc (octave_value::maybe_economize): New method.
+ (octave_value::non_null_value): rename to storable_value.
+ (octave_value::make_non_null_value): rename to make_storable_value.
+ * ov-base.h (octave_base_value::maybe_economize): New method.
+ * ov-base-mat.h (octave_base_mat::maybe_economize): New override.
+ * oct-obj.cc (octave_value_list::normalize_null_values):
+ Rename to make_storable_values, use make_storable_value.
+ * oct-obj.h: Dtto.
+ * ov-builtin.cc: non_null_value -> storable_value.
+ * ov-cell.cc: Dtto.
+ * ov-struct.cc: Dtto.
+ * pt-decl.h: Dtto.
+
diff --git a/src/oct-obj.cc b/src/oct-obj.cc
--- a/src/oct-obj.cc
+++ b/src/oct-obj.cc
@@ -236,11 +236,11 @@
 }
 
 void
-octave_value_list::normalize_null_values (void)
+octave_value_list::make_storable_values (void)
 {
   octave_idx_type len = length ();
   for (octave_idx_type i = 0; i < len; i++)
-    data[i].make_non_null_value ();
+    data[i].make_storable_value ();
 }
 
 /*
diff --git a/src/oct-obj.h b/src/oct-obj.h
--- a/src/oct-obj.h
+++ b/src/oct-obj.h
@@ -121,7 +121,7 @@
 
   string_vector name_tags (void) const { return names; }
 
-  void normalize_null_values (void);
+  void make_storable_values (void);
 
 private:
 
diff --git a/src/ov-base-mat.h b/src/ov-base-mat.h
--- a/src/ov-base-mat.h
+++ b/src/ov-base-mat.h
@@ -73,6 +73,8 @@
   octave_value squeeze (void) const { return MT (matrix.squeeze ()); }
 
   octave_value full_value (void) const { return matrix; }
+
+  void maybe_economize (void) { matrix.maybe_economize (); }
 
   octave_value subsref (const std::string& type,
  const std::list<octave_value_list>& idx);
diff --git a/src/ov-base.h b/src/ov-base.h
--- a/src/ov-base.h
+++ b/src/ov-base.h
@@ -151,6 +151,8 @@
   virtual octave_value full_value (void) const;
 
   virtual octave_base_value *try_narrowing_conversion (void) { return 0; }
+
+  virtual void maybe_economize (void) { }
 
   virtual octave_value
   subsref (const std::string& type,
diff --git a/src/ov-builtin.cc b/src/ov-builtin.cc
--- a/src/ov-builtin.cc
+++ b/src/ov-builtin.cc
@@ -107,7 +107,7 @@
   retval = (*f) (args, nargout);
           // Do not allow null values to be returned from functions.
           // FIXME -- perhaps true builtins should be allowed?
-          retval.normalize_null_values ();
+          retval.make_storable_values ();
  }
       catch (octave_execution_exception)
  {
diff --git a/src/ov-cell.cc b/src/ov-cell.cc
--- a/src/ov-cell.cc
+++ b/src/ov-cell.cc
@@ -270,7 +270,7 @@
       }
     else if (i.all_scalars () || do_index_op (i, true).numel () == 1)
               // Regularize a null matrix if stored into a cell.
-              octave_base_matrix<Cell>::assign (i, Cell (t_rhs.non_null_value ()));
+              octave_base_matrix<Cell>::assign (i, Cell (t_rhs.storable_value ()));
             else if (! error_state)
               error ("scalar indices required for {} in assignment.");
 
diff --git a/src/ov-struct.cc b/src/ov-struct.cc
--- a/src/ov-struct.cc
+++ b/src/ov-struct.cc
@@ -412,7 +412,7 @@
       }
     else
               // Regularize a null matrix if stored into a struct component.
-      map.assign (key, t_rhs.non_null_value ());
+      map.assign (key, t_rhs.storable_value ());
 
     if (! error_state)
       {
diff --git a/src/ov.cc b/src/ov.cc
--- a/src/ov.cc
+++ b/src/ov.cc
@@ -1170,7 +1170,7 @@
 {
   if (op == op_asn_eq)
     // Regularize a null matrix if stored into a variable.
-    operator = (rhs.non_null_value ());
+    operator = (rhs.storable_value ());
   else
     {
       // FIXME -- only do the following stuff if we can't find
@@ -1552,18 +1552,23 @@
                                              type_name (), "complex vector"));
 }
 
+// FIXME: This is a good place for pre-storage hooks, but the functions should
+// probably be named differently. These functions will be called
 
 octave_value
-octave_value::non_null_value (void) const
+octave_value::storable_value (void) const
 {
+  octave_value retval = *this;
   if (is_null_value ())
-    return octave_value (rep->empty_clone ());
+    retval = octave_value (rep->empty_clone ());
   else
-    return *this;
+    retval.maybe_economize ();
+
+  return retval;
 }
 
 void
-octave_value::make_non_null_value (void)
+octave_value::make_storable_value (void)
 {
   if (is_null_value ())
     {
@@ -1572,6 +1577,8 @@
         delete rep;
       rep = rc;
     }
+  else
+    maybe_economize ();
 }
 
 int
diff --git a/src/ov.h b/src/ov.h
--- a/src/ov.h
+++ b/src/ov.h
@@ -847,13 +847,18 @@
   Array<FloatComplex> float_complex_vector_value (bool frc_str_conv = false,
        bool frc_vec_conv = false) const;
 
-  // Make a copy that is not a special null matrix
+  // Possibly economize a lazy-indexed value.
 
-  octave_value non_null_value (void) const;
+  void maybe_economize (void)
+    { rep->maybe_economize (); }
+
+  // Make a copy suitable for storing.
+
+  octave_value storable_value (void) const;
 
   // Ditto, but in place.
 
-  void make_non_null_value (void);
+  void make_storable_value (void);
 
   // Conversions.  These should probably be private.  If a user of this
   // class wants a certain kind of constant, he should simply ask for
diff --git a/src/pt-decl.h b/src/pt-decl.h
--- a/src/pt-decl.h
+++ b/src/pt-decl.h
@@ -66,7 +66,7 @@
   bool lvalue_ok (void) { return id ? id->lvalue_ok () : false; }
 
   // Do not allow functions return null values
-  octave_value rvalue (void) { return id ? id->rvalue ().non_null_value () : octave_value (); }
+  octave_value rvalue (void) { return id ? id->rvalue ().storable_value () : octave_value (); }
 
   octave_value_list rvalue (int nargout)
   {

---

[arclp.diff]

---

# HG changeset patch
# User Jaroslav Hajek <highegg@...>
# Date 1231946505 -3600
# Node ID a62a982163b6ec855cf6d360a9e106e2ac14d342
# Parent  f3a54bae02c8e81cd79a8ceac45065f8af07f7fc
clean up Array and DiagArray2

diff --git a/liboctave/Array.cc b/liboctave/Array.cc
--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -46,6 +46,26 @@
 // all the derived classes.
 
 template <class T>
+void
+Array<T>::make_unique (void)
+{
+  if (rep->count > 1)
+    {
+      --rep->count;
+      rep = new ArrayRep (slice_data, slice_len, true);
+      slice_data = rep->data;
+    }
+  else if (slice_len != rep->len)
+    {
+      // Possibly economize here.
+      ArrayRep *new_rep = new ArrayRep (slice_data, slice_len, true);
+      delete rep;
+      rep = new_rep;
+      slice_data = rep->data;
+    }
+}
+
+template <class T>
 Array<T>::Array (const Array<T>& a, const dim_vector& dv)
   : rep (a.rep), dimensions (dv),
     slice_data (a.slice_data), slice_len (a.slice_len)
@@ -82,6 +102,20 @@
     }
 
   return *this;
+}
+
+template <class T>
+void
+Array<T>::fill (const T& val)
+{
+  if (rep->count > 1)
+    {
+      --rep->count;
+      rep = new ArrayRep (length (), val);
+      slice_data = rep->data;
+    }
+  else
+    std::fill (slice_data, slice_data + slice_len, val);
 }
 
 template <class T>
@@ -131,13 +165,7 @@
     }
  }
 
-      // FIXME -- it would be better if we did not have to do
-      // this, so we could share the data while still having different
-      // dimension vectors.
-
-      retval.make_unique ();
-
-      retval.dimensions = new_dimensions;
+      retval = Array<T> (*this, new_dimensions);
     }
 
   return retval;
diff --git a/liboctave/Array.h b/liboctave/Array.h
--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -97,49 +97,15 @@
 
 public:
 
-  // !!! WARNING !!! -- these should be protected, not public.  You
-  // should not access these methods directly!
-
-  void make_unique (void)
-    {
-      if (rep->count > 1)
- {
-  --rep->count;
-  rep = new ArrayRep (slice_data, slice_len, true);
-          slice_data = rep->data;
- }
-      else if (slice_len != rep->len)
-        {
-          // Possibly economize here.
-          ArrayRep *new_rep = new ArrayRep (slice_data, slice_len, true);
-          delete rep;
-          rep = new_rep;
-          slice_data = rep->data;
-        }
-      }
-
-  void make_unique (const T& val)
-    {
-      if (rep->count > 1)
- {
-  --rep->count;
-  rep = new ArrayRep (slice_len, val);
-          slice_data = rep->data;
- }
-      else
-        std::fill (slice_data, slice_data + slice_len, val);
-    }
+  void make_unique (void);
 
   typedef T element_type;
 
-  // !!! WARNING !!! -- these should be protected, not public.  You
-  // should not access these data members directly!
+protected:
 
   typename Array<T>::ArrayRep *rep;
 
   dim_vector dimensions;
-
-protected:
 
   T* slice_data;
   octave_idx_type slice_len;
@@ -220,7 +186,7 @@
   template <class U>
   Array (const Array<U>& a)
     : rep (new typename Array<T>::ArrayRep (coerce (a.data (), a.length ()), a.length ())),
-      dimensions (a.dimensions)
+      dimensions (a.dims ())
     {
       slice_data = rep->data;
       slice_len = rep->len;
@@ -260,7 +226,7 @@
 
   Array<T>& operator = (const Array<T>& a);
 
-  void fill (const T& val) { make_unique (val); }
+  void fill (const T& val);
 
   octave_idx_type capacity (void) const { return slice_len; }
   octave_idx_type length (void) const { return capacity (); }
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -0,0 +1,14 @@
[phone]  Jaroslav Hajek  <highegg@...>
+
+ * Array.h (Array<T>::rep, Array<T>::dimensions): Make protected.
+ * Array.cc (Array<T>::make_unique): Move implementation here.
+ (Array<T>::fill): Dtto.
+ * DiagArray2.h (DiagArray2<T>): Reimplement without abusing
+ Array<T> internals.
+ (DiagArray2<T>::operator Array2<T>): New method.
+ * DiagArray2.cc (DiagArray2<T>): Update methods.
+ * MDiagArray2.h (MDiagArray2<T>::operator Array2<T>): Simplify.
+ * PermMatrix.h (PermMatrix): Reimplement without abusing
+ Array<T> internals.
+ * PermMatrix.cc (PermMatrix): Update methods.
+
diff --git a/liboctave/DiagArray2.cc b/liboctave/DiagArray2.cc
--- a/liboctave/DiagArray2.cc
+++ b/liboctave/DiagArray2.cc
@@ -37,17 +37,42 @@
 #include "lo-error.h"
 
 template <class T>
+const typename DiagArray2<T>::Proxy&
+DiagArray2<T>::Proxy::operator = (const T& val) const
+{
+  if (i == j)
+    {
+      if (object)
+        object->set (val, i);
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("invalid assignment to off-diagonal in diagonal array");
+
+  return *this;
+}
+
+template <class T>
+DiagArray2<T>::Proxy::operator T () const
+{
+  if (object && i == j)
+    return object->get (i);
+  else
+    {
+      static T foo;
+      return foo;
+    }
+}
+
+template <class T>
 Array<T>
 DiagArray2<T>::diag (octave_idx_type k) const
 {
   Array<T> d;
 
   if (k == 0)
-    {
-      // The main diagonal is shallow-copied.
-      d = *this;
-      d.dimensions = dim_vector (length ());
-    }
+    // The main diagonal is shallow-copied.
+    d = *this;
   else if (k > 0 && k < cols ())
     d = Array<T> (std::min (cols () - k, rows ()), T ());
   else if (k < 0 && -k < rows ())
@@ -64,7 +89,8 @@
 DiagArray2<T>::transpose (void) const
 {
   DiagArray2<T> retval (*this);
-  retval.dimensions = dim_vector (dim2 (), dim1 ());
+  retval.d1 = d2;
+  retval.d2 = d1;
   return retval;
 }
 
@@ -91,7 +117,20 @@
       (*current_liboctave_error_handler) ("range error in DiagArray2");
       return T ();
     }
-  return (r == c) ? Array<T>::xelem (r) : T (0);
+  return elem (r, c);
+}
+
+template <class T>
+typename DiagArray2<T>::Proxy
+DiagArray2<T>::checkelem (octave_idx_type r, octave_idx_type c)
+{
+  if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
+    {
+      (*current_liboctave_error_handler) ("range error in DiagArray2");
+      return Proxy (0, r, c);
+    }
+  else
+    return Proxy (this, r, c);
 }
 
 template <class T>
@@ -104,37 +143,16 @@
       return;
     }
 
-  if (r == dim1 () && c == dim2 ())
-    return;
-
-  // FIXME: this is a mess. DiagArray2 really needs a rewrite.
-  typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
-  const T *old_data = this->data ();
-  octave_idx_type old_len = this->length ();
-
-  octave_idx_type new_len = r < c ? r : c;
-
-  Array<T>::rep = new typename Array<T>::ArrayRep (new_len);
-  Array<T>::slice_data = Array<T>::rep->data;
-  Array<T>::slice_len = Array<T>::rep->len;
-
-  this->dimensions = dim_vector (r, c);
-
-  if (old_data && old_len > 0)
+  if (r != dim1 () || c != dim2 ())
     {
-      octave_idx_type min_len = old_len < new_len ? old_len : new_len;
-
-      for (octave_idx_type i = 0; i < min_len; i++)
- xelem (i, i) = old_data[i];
+      Array<T>::resize (std::min (r, c));
+      d1 = r; d2 = c;
     }
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
 }
 
 template <class T>
 void
-DiagArray2<T>::resize (octave_idx_type r, octave_idx_type c, const T& val)
+DiagArray2<T>::resize_fill (octave_idx_type r, octave_idx_type c, const T& val)
 {
   if (r < 0 || c < 0)
     {
@@ -142,32 +160,21 @@
       return;
     }
 
-  if (r == dim1 () && c == dim2 ())
-    return;
+  if (r != dim1 () || c != dim2 ())
+    {
+      Array<T>::resize_fill (std::min (r, c), val);
+      d1 = r; d2 = c;
+    }
+}
 
-  typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
-  const T *old_data = this->data ();
-  octave_idx_type old_len = this->length ();
+template <class T>
+DiagArray2<T>::operator Array2<T> (void) const
+{
+  Array2<T> result (dim1 (), dim2 ());
+  for (octave_idx_type i = 0, len = length (); i < len; i++)
+    result.xelem (i, i) = dgelem (i);
 
-  octave_idx_type new_len = r < c ? r : c;
-
-  Array<T>::rep = new typename Array<T>::ArrayRep (new_len);
-
-  this->dimensions = dim_vector (r, c);
-
-  octave_idx_type min_len = old_len < new_len ? old_len : new_len;
-
-  if (old_data && old_len > 0)
-    {
-      for (octave_idx_type i = 0; i < min_len; i++)
- xelem (i, i) = old_data[i];
-    }
-
-  for (octave_idx_type i = min_len; i < new_len; i++)
-    xelem (i, i) = val;
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  return result;
 }
 
 /*
diff --git a/liboctave/DiagArray2.h b/liboctave/DiagArray2.h
--- a/liboctave/DiagArray2.h
+++ b/liboctave/DiagArray2.h
@@ -3,7 +3,7 @@
 
 Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007
               John W. Eaton
-Copyright (C) 2008 Jaroslav Hajek
+Copyright (C) 2008, 2009 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -30,10 +30,10 @@
 #include <cstdlib>
 
 #include "Array.h"
+#include "Array2.h"
 #include "lo-error.h"
 
 // A two-dimensional array with diagonal elements only.
-//
 // Idea and example code for Proxy class and functions from:
 //
 // From: kanze@... (James Kanze)
@@ -45,13 +45,12 @@
 // James Kanze                             email: kanze@...
 // GABI Software, Sarl., 8 rue du Faisan, F-67000 Strasbourg, France
 
-// Array<T> is inherited privately because we abuse the dimensions variable
-// for true dimensions. Therefore, the inherited Array<T> object is not a valid
-// Array<T> object, and should not be publicly accessible.
+// Array<T> is inherited privately so that some methods, like index, don't
+// produce unexpected results.
 
 template <class T>
 class
-DiagArray2 : private Array<T>
+DiagArray2 : protected Array<T>
 {
 private:
 
@@ -66,30 +65,9 @@
     Proxy (DiagArray2<T> *ref, octave_idx_type r, octave_idx_type c)
       : i (r), j (c), object (ref) { }
 
-    const Proxy& operator = (const T& val) const
-      {
- if (i == j)
-  {
-    if (object)
-      object->set (val, i);
-  }
- else
-  (*current_liboctave_error_handler)
-    ("invalid assignment to off-diagonal in diagonal array");
+    const Proxy& operator = (const T& val) const;
 
- return *this;
-      }
-
-    operator T () const
-      {
- if (object && i == j)
-  return object->get (i);
- else
-  {
-    static T foo;
-    return foo;
-  }
-      }
+    operator T () const;
 
   private:
 
@@ -106,100 +84,92 @@
 
   };
 
-friend class Proxy;
+  friend class Proxy;
 
 protected:
+  octave_idx_type d1, d2;
 
-  DiagArray2 (T *d, octave_idx_type r, octave_idx_type c) : Array<T> (d, r < c ? r : c)
-    { Array<T>::dimensions = dim_vector (r, c); }
+  DiagArray2 (T *d, octave_idx_type r, octave_idx_type c)
+    : Array<T> (d, std::min (r, c)), d1 (r), d2 (c) { }
 
 public:
 
   typedef T element_type;
 
-  DiagArray2 (void) : Array<T> (dim_vector (0, 0)) { }
+  DiagArray2 (void)
+    : Array<T> (), d1 (0), d2 (0) { }
 
-  DiagArray2 (octave_idx_type r, octave_idx_type c) : Array<T> (r < c ? r : c)
-    { this->dimensions = dim_vector (r, c); }
+  DiagArray2 (octave_idx_type r, octave_idx_type c)
+    : Array<T> (std::min (r, c)), d1 (r), d2 (c) { }
 
-  DiagArray2 (octave_idx_type r, octave_idx_type c, const T& val) : Array<T> (r < c ? r : c)
-    {
-      this->dimensions = dim_vector (r, c);
+  DiagArray2 (octave_idx_type r, octave_idx_type c, const T& val)
+    : Array<T> (std::min (r, c), val), d1 (r), d2 (c) { }
 
-      Array<T>::fill (val);
-    }
+  DiagArray2 (const Array<T>& a)
+    : Array<T> (a), d1 (a.numel ()), d2 (a.numel ()) { }
 
-  DiagArray2 (const Array<T>& a) : Array<T> (a)
-    { this->dimensions = dim_vector (a.length (), a.length ()); }
-
-  DiagArray2 (const DiagArray2<T>& a) : Array<T> (a)
-    { this->dimensions = a.dims (); }
+  DiagArray2 (const DiagArray2<T>& a)
+    : Array<T> (a), d1 (a.d1), d2 (a.d2) { }
 
   template <class U>
-  DiagArray2 (const DiagArray2<U>& a) : Array<T> (a.diag ())
-    { this->dimensions = a.dims (); }
+  DiagArray2 (const DiagArray2<U>& a)
+  : Array<T> (a.diag ()), d1 (a.dim1 ()), d2 (a.dim2 ()) { }
 
   ~DiagArray2 (void) { }
 
   DiagArray2<T>& operator = (const DiagArray2<T>& a)
     {
       if (this != &a)
- Array<T>::operator = (a);
+        {
+          Array<T>::operator = (a);
+          d1 = a.d1;
+          d2 = a.d2;
+        }
 
       return *this;
     }
 
-
-  octave_idx_type dim1 (void) const { return Array<T>::dimensions(0); }
-  octave_idx_type dim2 (void) const { return Array<T>::dimensions(1); }
+  octave_idx_type dim1 (void) const { return d1; }
+  octave_idx_type dim2 (void) const { return d2; }
 
   octave_idx_type rows (void) const { return dim1 (); }
   octave_idx_type cols (void) const { return dim2 (); }
   octave_idx_type columns (void) const { return dim2 (); }
 
+  // FIXME: a dangerous ambiguity?
   octave_idx_type length (void) const { return Array<T>::length (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
   octave_idx_type numel (void) const { return nelem (); }
 
   size_t byte_size (void) const { return length () * sizeof (T); }
 
-  dim_vector dims (void) const { return Array<T>::dimensions; }
+  dim_vector dims (void) const { return dim_vector (d1, d2); }
 
   Array<T> diag (octave_idx_type k = 0) const;
 
-  Proxy elem (octave_idx_type r, octave_idx_type c)
-    {
-      return Proxy (this, r, c);
-    }
-
-  Proxy checkelem (octave_idx_type r, octave_idx_type c)
-    {
-      if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
- {
-  (*current_liboctave_error_handler) ("range error in DiagArray2");
-  return Proxy (0, r, c);
- }
-      else
- return Proxy (this, r, c);
-    }
-
-  Proxy operator () (octave_idx_type r, octave_idx_type c)
-    {
-      if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
- {
-  (*current_liboctave_error_handler) ("range error in DiagArray2");
-  return Proxy (0, r, c);
- }
-      else
- return Proxy (this, r, c);
-  }
+  // Warning: the non-const two-index versions will silently ignore assignments
+  // to off-diagonal elements.
 
   T elem (octave_idx_type r, octave_idx_type c) const
     {
-      return (r == c) ? Array<T>::xelem (r) : T (0);
+      return (r == c) ? Array<T>::elem (r) : T (0);
     }
 
+  T& elem (octave_idx_type r, octave_idx_type c)
+    {
+      static T zero (0);
+      return (r == c) ? Array<T>::elem (r) : zero;
+    }
+
+  T dgelem (octave_idx_type i) const
+    { return Array<T>::elem (i); }
+
+  T& dgelem (octave_idx_type i)
+    { return Array<T>::elem (i); }
+
   T checkelem (octave_idx_type r, octave_idx_type c) const;
+  Proxy checkelem (octave_idx_type r, octave_idx_type c);
+
   T operator () (octave_idx_type r, octave_idx_type c) const
     {
 #if defined (BOUNDS_CHECKING)
@@ -209,24 +179,39 @@
 #endif
     }
 
+  // FIXME: can this cause problems?
+#if defined (BOUNDS_CHECKING)
+  Proxy operator () (octave_idx_type r, octave_idx_type c)
+    {
+      return checkelem (r, c);
+    }
+#else
+  T& operator () (octave_idx_type r, octave_idx_type c)
+    {
+      return elem (r, c);
+    }
+#endif
+
   // No checking.
-
-  T& xelem (octave_idx_type r, octave_idx_type c)
-    {
-      static T foo (0);
-      return (r == c) ? Array<T>::xelem (r) : foo;
-    }
 
   T xelem (octave_idx_type r, octave_idx_type c) const
     {
       return (r == c) ? Array<T>::xelem (r) : T (0);
     }
 
+  T& dgxelem (octave_idx_type i)
+    { return Array<T>::xelem (i); }
+
+  T dgxelem (octave_idx_type i) const
+    { return Array<T>::xelem (i); }
+
   void resize (octave_idx_type n, octave_idx_type m);
-  void resize (octave_idx_type n, octave_idx_type m, const T& val);
+  void resize_fill (octave_idx_type n, octave_idx_type m, const T& val);
 
   DiagArray2<T> transpose (void) const;
   DiagArray2<T> hermitian (T (*fcn) (const T&) = 0) const;
+
+  operator Array2<T> (void) const;
 
   const T *data (void) const { return Array<T>::data (); }
 
diff --git a/liboctave/MDiagArray2.h b/liboctave/MDiagArray2.h
--- a/liboctave/MDiagArray2.h
+++ b/liboctave/MDiagArray2.h
@@ -72,17 +72,7 @@
 
   operator MArray2<T> () const
     {
-      octave_idx_type nr = DiagArray2<T>::dim1 ();
-      octave_idx_type nc = DiagArray2<T>::dim2 ();
-
-      MArray2<T> retval (nr, nc,  T (0));
-
-      octave_idx_type len = nr < nc ? nr : nc;
-
-      for (octave_idx_type i = 0; i < len; i++)
- retval.xelem (i, i) = this->xelem (i, i);
-
-      return retval;
+      return DiagArray2<T>::operator Array2<T> ();
     }
 
   octave_idx_type nnz (void) const
diff --git a/liboctave/PermMatrix.cc b/liboctave/PermMatrix.cc
--- a/liboctave/PermMatrix.cc
+++ b/liboctave/PermMatrix.cc
@@ -39,7 +39,6 @@
 PermMatrix::PermMatrix (const Array<octave_idx_type>& p, bool colp, bool check)
   : Array<octave_idx_type> (p), _colp(colp)
 {
-  this->dimensions = dim_vector (p.length (), p.length ());
   if (check)
     {
       if (! idx_vector (p).is_permutation (p.length ()))
@@ -61,14 +60,12 @@
       Array<octave_idx_type> idxa (len);
       for (octave_idx_type i = 0; i < len; i++) idxa(i) = idx(i);
       Array<octave_idx_type>::operator = (idxa);
-      this->dimensions = dim_vector (len, len);
     }
 }
 
 PermMatrix::PermMatrix (octave_idx_type n)
   : Array<octave_idx_type> (n), _colp (false)
 {
-  this->dimensions = dim_vector (n, n);
   for (octave_idx_type i = 0; i < n; i++) xelem (i) = i;
 }
 
diff --git a/liboctave/PermMatrix.h b/liboctave/PermMatrix.h
--- a/liboctave/PermMatrix.h
+++ b/liboctave/PermMatrix.h
@@ -26,15 +26,11 @@
 #include "Array.h"
 #include "mx-defs.h"
 
-// Array<T> is inherited privately because we abuse the dimensions variable
-// for true dimensions. Therefore, the inherited Array<T> object is not a valid
-// Array<T> object, and should not be publicly accessible.
+// Array<T> is inherited privately so that some methods, like index, don't
+// produce unexpected results.
 
-class PermMatrix : private Array<octave_idx_type>
+class PermMatrix : protected Array<octave_idx_type>
 {
-private:
-
-  octave_idx_type get (octave_idx_type i) const { return Array<octave_idx_type>::xelem (i); }
 
 public:
 
@@ -51,36 +47,34 @@
   PermMatrix (const idx_vector& idx, bool colp = false, octave_idx_type n = 0);
 
   octave_idx_type dim1 (void) const
-    { return Array<octave_idx_type>::dimensions(0); }
+    { return Array<octave_idx_type>::length (); }
   octave_idx_type dim2 (void) const
-    { return Array<octave_idx_type>::dimensions(1); }
+    { return Array<octave_idx_type>::length (); }
 
   octave_idx_type rows (void) const { return dim1 (); }
   octave_idx_type cols (void) const { return dim2 (); }
   octave_idx_type columns (void) const { return dim2 (); }
 
+  octave_idx_type perm_length (void) const
+    { return Array<octave_idx_type>::length (); }
   octave_idx_type length (void) const
-    { return Array<octave_idx_type>::length (); }
+    { return dim1 () * dim2 (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
   octave_idx_type numel (void) const { return nelem (); }
 
-  size_t byte_size (void) const { return length () * sizeof (octave_idx_type); }
+  size_t byte_size (void) const { return perm_length () * sizeof (octave_idx_type); }
 
-  dim_vector dims (void) const { return Array<octave_idx_type>::dimensions; }
+  dim_vector dims (void) const { return dim_vector (dim1 (), dim2 ()); }
 
   Array<octave_idx_type> pvec (void) const
-    {
-      Array<octave_idx_type> retval (*this);
-      retval.dimensions = dim_vector (length ());
-      return retval;
-    }
+    { return *this; }
 
   octave_idx_type
   elem (octave_idx_type i, octave_idx_type j) const
     {
       return (_colp
-              ? ((get(j) != i) ? 1 : 0)
-              : ((get(i) != j) ? 1 : 0));
+              ? ((Array<octave_idx_type>::elem (j) != i) ? 1 : 0)
+              : ((Array<octave_idx_type>::elem (i) != j) ? 1 : 0));
     }
 
   octave_idx_type
diff --git a/src/strfns.cc b/src/strfns.cc
--- a/src/strfns.cc
+++ b/src/strfns.cc
@@ -395,7 +395,7 @@
     {
       // Broadcast the string.
 
-      boolNDArray output (cell.dimensions, false);
+      boolNDArray output (cell.dims (), false);
 
       std::string s = r == 0 ? std::string () : str[0];
 
@@ -430,7 +430,7 @@
  {
   // Must match in all dimensions.
 
-  boolNDArray output (cell.dimensions, false);
+  boolNDArray output (cell.dims (), false);
 
   if (cell.length () == r)
     {
@@ -471,8 +471,8 @@
       r2 = cell2.length ();
     }
 
-  const dim_vector size1 = cell1.dimensions;
-  const dim_vector size2 = cell2.dimensions;
+  const dim_vector size1 = cell1.dims ();
+  const dim_vector size2 = cell2.dims ();
 
   boolNDArray output (size1, false);
 
@@ -671,7 +671,7 @@
     {
       // Broadcast the string.
 
-      boolNDArray output (cell.dimensions, false);
+      boolNDArray output (cell.dims (), false);
 
       if (c < n)
  {
@@ -724,7 +724,7 @@
  {
   // Must match in all dimensions.
 
-  boolNDArray output (cell.dimensions, false);
+  boolNDArray output (cell.dims (), false);
 
   if (cell.numel () == r)
     {
@@ -774,8 +774,8 @@
       r2 = cell2.length ();
     }
 
-  const dim_vector size1 = cell1.dimensions;
-  const dim_vector size2 = cell2.dimensions;
+  const dim_vector size1 = cell1.dims ();
+  const dim_vector size2 = cell2.dims ();
 
   boolNDArray output (size1, false);
 

---

[John W. Eaton-3]

On 14-Jan-2009, Jaroslav Hajek wrote:

| Summary (1st patch):
| If a dense indexing expression can be reduced to a contiguous subrange
| (like x(2:end-1), A(:,:,3:10) etc.), a shared (i.e. cheap) copy is
| made and the actual extraction is postponed until (if ever) the value
| is written to. The obvious benefit is greater performance of various
| operations with indexed expressions, such as y = x(1:n-1) + x(2:n),
| which will avoid making two temporary copies. This is done by
| modifying Array<T> to allow it work with contiguous subranges of its
| (possibly shared) data. There is an obvious risk of increased memory
| consumption, because a big block of memory may be kept allocated while
| only part of it is actually used.
| This is resolved by checking prior to storing a value to a "permanent"
| location (same where checks for null matrices are done), and possibly
| economizing if an orphaned block is detected. The detection (see
| Array<T>::maybe_economize) is fairly simple and could be improved at
| the cost of performance of the operation. Since most indexing
| expressions do not live longer than their parent object I suppose this
| will seldom be an issue, but it is still easily possible to construct
| an example where memory is wastefully kept allocated.
| However, the issue is easy to avoid if you know what you do, so I
| think it's OK. Comments are welcome.
|
| To get an idea of the performance speed up, here's a simplistic test
| (set n to suitable value):
|
| n = 2e7;
| x = 1/n * [1:n]; y = x.^2;
| # compute 2nd-order differences
| tic; d2y = diff (y, 2); toc
| # compute integral
| tic; int = trapz (x, y); toc
|
| with a recent tip, I get:
|
| Elapsed time is 0.697262 seconds.
| Elapsed time is 0.960276 seconds.
|
| with the new patch, I got:
|
| Elapsed time is 0.213054 seconds.
| Elapsed time is 0.497459 seconds.
|
| which means speedups by 228% and 93%.
| P.S. I picked some benefitted functions because measuring directly the
| performance speed-up of the affected indexing expressions is
| meaningless.
|
|
| Summary (2nd patch):
| This is a long-postponed cleanup of Array<T>, DiagArray2<T> and
| PermMatrix, to allow the private members of Array<T> to become really
| private,
| and reimplement DiagArray and PermMatrix without the dangerous abuse
| of Array<T>::dimensions (that has caused me quite a pain with the
| previous patch).
| This patch should probably be applied in any case, I'm putting it here
| because right now it's created on top of the 1st patch and I've only
| ran make check on the result.
|
| comments? OK to push?

I think these changes are fine.  I see that you already pushed them,
which is OK, but maybe in the future when you ask for comments you
could wait a few days?  :-)

I noticed the following incomplete comment in ov.cc:

  // FIXME: This is a good place for pre-storage hooks, but the functions should
  // probably be named differently. These functions will be called

Can you please explain (in a comment somewhere) what is meant by
"storable_value"?

Also, it might be good to add some explanation in Array.h about how
slice_data is supposed to be used.  I'd like to try to avoid confusion
for people who might want to make modifications to the Array class in
the future.

Thanks,

jwe

[John W. Eaton-3]

On 14-Jan-2009, Jaroslav Hajek wrote:

| with a recent tip, I get:
|
| Elapsed time is 0.697262 seconds.
| Elapsed time is 0.960276 seconds.
|
| with the new patch, I got:
|
| Elapsed time is 0.213054 seconds.
| Elapsed time is 0.497459 seconds.
|
| which means speedups by 228% and 93%.
| P.S. I picked some benefitted functions because measuring directly the
| performance speed-up of the affected indexing expressions is
| meaningless.

Oh, and I should mention that it is great to have the speedup.

Thanks!

jwe

[Jaroslav Hajek-2]

On Fri, Jan 16, 2009 at 10:24 PM, John W. Eaton <jwe@...> wrote:

> On 14-Jan-2009, Jaroslav Hajek wrote:
>
> | Summary (1st patch):
> | If a dense indexing expression can be reduced to a contiguous subrange
> | (like x(2:end-1), A(:,:,3:10) etc.), a shared (i.e. cheap) copy is
> | made and the actual extraction is postponed until (if ever) the value
> | is written to. The obvious benefit is greater performance of various
> | operations with indexed expressions, such as y = x(1:n-1) + x(2:n),
> | which will avoid making two temporary copies. This is done by
> | modifying Array<T> to allow it work with contiguous subranges of its
> | (possibly shared) data. There is an obvious risk of increased memory
> | consumption, because a big block of memory may be kept allocated while
> | only part of it is actually used.
> | This is resolved by checking prior to storing a value to a "permanent"
> | location (same where checks for null matrices are done), and possibly
> | economizing if an orphaned block is detected. The detection (see
> | Array<T>::maybe_economize) is fairly simple and could be improved at
> | the cost of performance of the operation. Since most indexing
> | expressions do not live longer than their parent object I suppose this
> | will seldom be an issue, but it is still easily possible to construct
> | an example where memory is wastefully kept allocated.
> | However, the issue is easy to avoid if you know what you do, so I
> | think it's OK. Comments are welcome.
> |
> | To get an idea of the performance speed up, here's a simplistic test
> | (set n to suitable value):
> |
> | n = 2e7;
> | x = 1/n * [1:n]; y = x.^2;
> | # compute 2nd-order differences
> | tic; d2y = diff (y, 2); toc
> | # compute integral
> | tic; int = trapz (x, y); toc
> |
> | with a recent tip, I get:
> |
> | Elapsed time is 0.697262 seconds.
> | Elapsed time is 0.960276 seconds.
> |
> | with the new patch, I got:
> |
> | Elapsed time is 0.213054 seconds.
> | Elapsed time is 0.497459 seconds.
> |
> | which means speedups by 228% and 93%.
> | P.S. I picked some benefitted functions because measuring directly the
> | performance speed-up of the affected indexing expressions is
> | meaningless.
> |
> |
> | Summary (2nd patch):
> | This is a long-postponed cleanup of Array<T>, DiagArray2<T> and
> | PermMatrix, to allow the private members of Array<T> to become really
> | private,
> | and reimplement DiagArray and PermMatrix without the dangerous abuse
> | of Array<T>::dimensions (that has caused me quite a pain with the
> | previous patch).
> | This patch should probably be applied in any case, I'm putting it here
> | because right now it's created on top of the 1st patch and I've only
> | ran make check on the result.
> |
> | comments? OK to push?
>
> I think these changes are fine.  I see that you already pushed them,
> which is OK, but maybe in the future when you ask for comments you
> could wait a few days?  :-)
>

Sorry. Basically the impulse to push it was that I wanted to refer to
the first changeset in a discussion (on this ML). A stupid reason, of
course.
My apologies, once more.

> I noticed the following incomplete comment in ov.cc:
>
>  // FIXME: This is a good place for pre-storage hooks, but the functions should
>  // probably be named differently. These functions will be called
>
> Can you please explain (in a comment somewhere) what is meant by
> "storable_value"?
>

storable_value and make_storable_value is a hook that is called on an
octave_value
prior to storing it to a "permanent" location, like a named variable,
a cell or a struct component, or a return value of a function.
Previously, this was named non_null_value because all it did was to
canonicalize the special null matrices (for element deletion). I
realized that the proper places to check for maybe_economize are the
same, so I renamed the hook function to reflect
*where it's being called* rather than what it does. If, later, some
other "lazy" mechanism is employed in Octave, then checks for possible
"unlazifying" will probably also go in the same place.
A suggestion for a better name or better scheme is welcome. Ill put in
an explaining comment.


> Also, it might be good to add some explanation in Array.h about how
> slice_data is supposed to be used.  I'd like to try to avoid confusion
> for people who might want to make modifications to the Array class in
> the future.

OK, I'll do it.

cheers

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[John W. Eaton-3]

On 17-Jan-2009, Jaroslav Hajek wrote:

| My apologies, once more.

There's no need to apologize.  Sometimes its difficult to tell whether
people are slow to comment or just not planning to comment...

| other "lazy" mechanism is employed in Octave, then checks for possible
| "unlazifying" will probably also go in the same place.
| A suggestion for a better name or better scheme is welcome. Ill put in
| an explaining comment.

OK.  I don't have a better name at the moment.

Thanks,

jwe