--- a/configure.ac	Sat Jun 21 13:13:05 2014 -0600
+++ b/configure.ac	Mon Jul 14 13:07:59 2014 -0600
@@ -231,6 +231,19 @@
   AC_DEFINE(BOUNDS_CHECKING, 1, [Define to 1 to use internal bounds checking.])
 fi
 
+### Enable uniqueness-checking on array modifications
+### As with bounds-checking, this slows down some things a bit so it is turned
+### off by default
+
+UNIQUENESS_CHECKING=no
+AC_ARG_ENABLE([uniqueness-check],
+  [AS_HELP_STRING([--enable-uniqueness-check],
+    [enable uniqueness checking for direct array modification])],
+  [if test "$enableval" = yes; then UNIQUENESS_CHECKING=yes; fi], [])
+if test $UNIQUENESS_CHECKING = yes; then
+  AC_DEFINE(UNIQUENESS_CHECKING, 1, [Define to 1 to use internal uniqueness checking.])
+fi
+
 ### Use Octave's built-in memory allocator rather than straightforward malloc.
 ### Disabled by default.
 

--- a/libinterp/corefcn/jit-typeinfo.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/libinterp/corefcn/jit-typeinfo.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -286,7 +286,7 @@
       Array<idx_vector> idx;
       make_indices (indicies, idx_count, idx);
 
-      Array<double> ret = mat->array->index (idx);
+      const Array<double> ret = mat->array->index (idx);
       return ret.xelem (0);
     }
   catch (const octave_execution_exception&)

--- a/liboctave/array/Array-util.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/array/Array-util.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -26,6 +26,7 @@
 #endif
 
 #include "Array-util.h"
+#include "Array.h"
 #include "dim-vector.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
@@ -174,60 +175,55 @@
   return retval;
 }
 
-octave_idx_type
-compute_index (octave_idx_type n, const dim_vector& dims)
+void
+check_index_bounds (int nd, int dim, octave_idx_type i, octave_idx_type size)
 {
-  if (n < 0)
+  if (i < 0)
     gripe_invalid_index ();
-  if (n >= dims.numel ())
-    gripe_index_out_of_range (1, 1, n+1, dims.numel ());
+  if (i >= size)
+    gripe_index_out_of_range (nd, dim + 1, i + 1, size);
+}
 
-  return n;
+void
+check_index (octave_idx_type n, const dim_vector& dims)
+{
+  check_index_bounds (1, 0, n, dims.numel ());
 }
 
-octave_idx_type
-compute_index (octave_idx_type i, octave_idx_type j, const dim_vector& dims)
+void
+check_index (octave_idx_type i, octave_idx_type j, const dim_vector& dims)
 {
-  if (i < 0 || j < 0)
-    gripe_invalid_index ();
-  if (i >= dims(0))
-    gripe_index_out_of_range (2, 1, i+1, dims(0));
-  if (j >= dims.numel (1))
-    gripe_index_out_of_range (2, 2, j+1, dims.numel (1));
-
-  return j*dims(0) + i;
+  check_index_bounds (2, 0, i, dims(0));
+  check_index_bounds (2, 1, j, dims.numel (1));
 }
 
-octave_idx_type
-compute_index (octave_idx_type i, octave_idx_type j, octave_idx_type k,
-               const dim_vector& dims)
+void
+check_index (octave_idx_type i, octave_idx_type j, octave_idx_type k,
+             const dim_vector& dims)
 {
-  if (i < 0 || j < 0 || k < 0)
-    gripe_invalid_index ();
-  if (i >= dims(0))
-    gripe_index_out_of_range (3, 1, i+1, dims(0));
-  if (j >= dims(1))
-    gripe_index_out_of_range (3, 2, j+1, dims(1));
-  if (k >= dims.numel (2))
-    gripe_index_out_of_range (3, 3, k+1, dims.numel (2));
-
-  return (k*dims(1) + j)*dims(0) + i;
+  check_index_bounds (3, 0, i, dims(0));
+  check_index_bounds (3, 1, j, dims(1));
+  check_index_bounds (3, 2, k, dims.numel (2));
 }
 
-octave_idx_type
-compute_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims)
+void
+check_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims)
 {
   int nd = ra_idx.length ();
   const dim_vector dv = dims.redim (nd);
   for (int d = 0; d < nd; d++)
     {
-      if (ra_idx(d) < 0)
-        gripe_invalid_index ();
-      if (ra_idx(d) >= dv(d))
-        gripe_index_out_of_range (nd, d+1, ra_idx(d)+1, dv(d));
+      check_index_bounds (nd, d, ra_idx(d), dv(d));
     }
+}
 
-  return dv.compute_index (ra_idx.data ());
+octave_idx_type
+compute_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims,
+               bool check)
+{
+  if (BOUNDS_CHECKING_DEFINED || check)
+    check_index (ra_idx, dims);
+  return dims.compute_index (ra_idx.data ());
 }
 
 Array<octave_idx_type>

--- a/liboctave/array/Array-util.h	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/array/Array-util.h	Mon Jul 14 13:07:59 2014 -0600
@@ -23,11 +23,13 @@
 #if !defined (octave_Array_util_h)
 #define octave_Array_util_h 1
 
-#include "Array.h"
 #include "dim-vector.h"
 #include "idx-vector.h"
 #include "lo-array-gripes.h"
 
+template <typename T>
+class Array;
+
 extern OCTAVE_API bool index_in_bounds (const Array<octave_idx_type>& ra_idx,
                                         const dim_vector& dimensions);
 
@@ -47,20 +49,59 @@
 
 extern OCTAVE_API bool any_ones (const Array<octave_idx_type>& arr);
 
+extern OCTAVE_API void
+check_index_bounds (int nd, int dim, octave_idx_type i, octave_idx_type size);
+
+extern OCTAVE_API void
+check_index (octave_idx_type n, const dim_vector& dims);
+
+extern OCTAVE_API void
+check_index (octave_idx_type i, octave_idx_type j, const dim_vector& dims);
+
+extern OCTAVE_API void
+check_index (octave_idx_type i, octave_idx_type j, octave_idx_type k,
+               const dim_vector& dims);
+
+extern OCTAVE_API void
+check_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims);
+
+#if defined(BOUNDS_CHECKING)
+#define BOUNDS_CHECKING_DEFINED true
+#else
+#define BOUNDS_CHECKING_DEFINED false
+#endif
+
 // These four compute a linear index for given dimensions, throwing
-// exceptions on invalid indices.
-extern OCTAVE_API octave_idx_type
-compute_index (octave_idx_type n, const dim_vector& dims);
+// exceptions on invalid indices if check is true or BOUNDS_CHECKING is defined.
+inline octave_idx_type
+compute_index (octave_idx_type n, const dim_vector& dims, bool check = true)
+{
+  if (BOUNDS_CHECKING_DEFINED || check)
+    check_index (n, dims);
+  return n;
+}
+
+inline octave_idx_type
+compute_index (octave_idx_type i, octave_idx_type j, const dim_vector& dims,
+               bool check = true)
+{
+  if (BOUNDS_CHECKING_DEFINED || check)
+    check_index (i, j, dims);
+  return j * dims(0) + i;
+}
+
+inline octave_idx_type
+compute_index (octave_idx_type i, octave_idx_type j, octave_idx_type k,
+               const dim_vector& dims, bool check = true)
+{
+  if (BOUNDS_CHECKING_DEFINED || check)
+    check_index (i, j, k, dims);
+  return (k * dims (1) + j) * dims(0) + i;
+}
 
 extern OCTAVE_API octave_idx_type
-compute_index (octave_idx_type i, octave_idx_type j, const dim_vector& dims);
-
-extern OCTAVE_API octave_idx_type
-compute_index (octave_idx_type i, octave_idx_type j, octave_idx_type k,
-               const dim_vector& dims);
-
-extern OCTAVE_API octave_idx_type
-compute_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims);
+compute_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims,
+               bool check = true);
 
 extern OCTAVE_API Array<octave_idx_type>
 conv_to_int_array (const Array<idx_vector>& a);

--- a/liboctave/array/Array.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/array/Array.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -36,7 +36,6 @@
 #include <new>
 
 #include "Array.h"
-#include "Array-util.h"
 #include "idx-vector.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
@@ -189,10 +188,7 @@
 Array<T>::checkelem (octave_idx_type n)
 {
   // Do checks directly to avoid recomputing slice_len.
-  if (n < 0)
-    gripe_invalid_index ();
-  if (n >= slice_len)
-    gripe_index_out_of_range (1, 1, n+1, slice_len);
+  check_index_bounds (1, 0, n, slice_len);
 
   return elem (n);
 }

--- a/liboctave/array/Array.h	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/array/Array.h	Mon Jul 14 13:07:59 2014 -0600
@@ -32,6 +32,7 @@
 #include <algorithm>
 #include <iosfwd>
 
+#include "Array-util.h"
 #include "dim-vector.h"
 #include "idx-vector.h"
 #include "lo-traits.h"
@@ -341,26 +342,46 @@
   const
   { return dimensions.compute_index (ra_idx.data (), ra_idx.length ()); }
 
-  // No checking, even for multiple references, ever.
+  // Check for multiple references only if uniqueness-checking is enabled
+  // Check for out-of-range only if bounds-checking is enabled
+  // Otherwise no checking
 
-  T& xelem (octave_idx_type n) { return slice_data[n]; }
-  crefT xelem (octave_idx_type n) const { return slice_data[n]; }
+  T&
+  xelem (octave_idx_type n)
+  {
+#if defined (UNIQUENESS_CHECKING)
+    if (!is_unique ())
+      gripe_modifying_nonunique ();
+#endif
+#if defined (BOUNDS_CHECKING)
+    check_index_bounds (1, 0, n, slice_len);
+#endif
+    return slice_data[n];
+  }
+  crefT
+  xelem (octave_idx_type n) const
+  {
+#if defined (BOUNDS_CHECKING)
+    check_index_bounds (1, 0, n, slice_len);
+#endif
+    return slice_data[n];
+  }
 
   T& xelem (octave_idx_type i, octave_idx_type j)
-  { return xelem (dim1 ()*j+i); }
+  { return xelem (::compute_index(i, j, dims (), false)); }
   crefT xelem (octave_idx_type i, octave_idx_type j) const
-  { return xelem (dim1 ()*j+i); }
+  { return xelem (::compute_index(i, j, dims (), false)); }
 
   T& xelem (octave_idx_type i, octave_idx_type j, octave_idx_type k)
-  { return xelem (i, dim2 ()*k+j); }
+  { return xelem (::compute_index(i, j, k, dims (), false)); }
   crefT xelem (octave_idx_type i, octave_idx_type j, octave_idx_type k) const
-  { return xelem (i, dim2 ()*k+j); }
+  { return xelem (::compute_index(i, j, k, dims (), false)); }
 
   T& xelem (const Array<octave_idx_type>& ra_idx)
-  { return xelem (compute_index_unchecked (ra_idx)); }
+  { return xelem (::compute_index(ra_idx, dims (), false)); }
 
   crefT xelem (const Array<octave_idx_type>& ra_idx) const
-  { return xelem (compute_index_unchecked (ra_idx)); }
+  { return xelem (::compute_index(ra_idx, dims (), false)); }
 
   // FIXME: would be nice to fix this so that we don't unnecessarily force
   //        a copy, but that is not so easy, and I see no clean way to do it.
@@ -721,6 +742,8 @@
 
 private:
 
+  bool is_unique (void) const { return rep->count == 1; }
+
   void resize2 (octave_idx_type nr, octave_idx_type nc, const T& rfv);
   void resize2 (octave_idx_type nr, octave_idx_type nc)
   {

--- a/liboctave/numeric/CmplxQR.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/numeric/CmplxQR.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -159,13 +159,14 @@
         {
           // afact will become q.
           q = afact;
+          const ComplexMatrix& c_afact = afact;
           octave_idx_type k = qr_type == qr_type_economy ? n : m;
           r = ComplexMatrix (k, n);
           for (octave_idx_type j = 0; j < n; j++)
             {
               octave_idx_type i = 0;
               for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
+                r.xelem (i, j) = c_afact.xelem (i, j);
               for (; i < k; i++)
                 r.xelem (i, j) = 0;
             }

--- a/liboctave/numeric/dbleQR.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/numeric/dbleQR.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -160,13 +160,14 @@
         {
           // afact will become q.
           q = afact;
+          const Matrix& c_afact = afact;
           octave_idx_type k = qr_type == qr_type_economy ? n : m;
           r = Matrix (k, n);
           for (octave_idx_type j = 0; j < n; j++)
             {
               octave_idx_type i = 0;
               for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
+                r.xelem (i, j) = c_afact.xelem (i, j);
               for (; i < k; i++)
                 r.xelem (i, j) = 0;
             }

--- a/liboctave/numeric/fCmplxQR.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/numeric/fCmplxQR.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -162,13 +162,14 @@
         {
           // afact will become q.
           q = afact;
+          const FloatComplexMatrix& c_afact = afact;
           octave_idx_type k = qr_type == qr_type_economy ? n : m;
           r = FloatComplexMatrix (k, n);
           for (octave_idx_type j = 0; j < n; j++)
             {
               octave_idx_type i = 0;
               for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
+                r.xelem (i, j) = c_afact.xelem (i, j);
               for (; i < k; i++)
                 r.xelem (i, j) = 0;
             }

--- a/liboctave/numeric/floatQR.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/numeric/floatQR.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -158,13 +158,14 @@
         {
           // afact will become q.
           q = afact;
+          const FloatMatrix& c_afact = afact;
           octave_idx_type k = qr_type == qr_type_economy ? n : m;
           r = FloatMatrix (k, n);
           for (octave_idx_type j = 0; j < n; j++)
             {
               octave_idx_type i = 0;
               for (; i <= j; i++)
-                r.xelem (i, j) = afact.xelem (i, j);
+                r.xelem (i, j) = c_afact.xelem (i, j);
               for (; i < k; i++)
                 r.xelem (i, j) = 0;
             }

--- a/liboctave/util/interp-idx.h	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/util/interp-idx.h	Mon Jul 14 13:07:59 2014 -0600
@@ -31,14 +31,17 @@
 }
 
 // Simple method for taking a row-column pair together with the matrix
-// height and returning the corresponding index as an octave data-type
+// dimensions and returning the corresponding index as an octave data-type
 // (note that for large heights, there's a risk of losing precision.
 // This method will not overflow or throw a bad alloc, it will simply
 // choose the nearest possible double-value to the proper index).
 inline double
-to_interp_idx (octave_idx_type row, octave_idx_type col, octave_idx_type height)
+to_interp_idx (octave_idx_type row, octave_idx_type col, const dim_vector& dims)
 {
-  return col * static_cast<double> (height) + row + 1;
+#if defined(BOUNDS_CHECKING)
+  check_index (row, col, dims);
+#endif
+  return col * static_cast<double> (dims(0)) + row + 1;
 }
 
 #endif

--- a/liboctave/util/lo-array-gripes.cc	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/util/lo-array-gripes.cc	Mon Jul 14 13:07:59 2014 -0600
@@ -166,3 +166,9 @@
     ("A(I,J,...) = X: dimensions mismatch");
 }
 
+void
+gripe_modifying_nonunique (void)
+{
+  (*current_liboctave_error_handler)
+    ("Attempting to directly modify array with multiple references");
+}

--- a/liboctave/util/lo-array-gripes.h	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/util/lo-array-gripes.h	Mon Jul 14 13:07:59 2014 -0600
@@ -71,4 +71,7 @@
 extern void OCTAVE_API
 gripe_assignment_dimension_mismatch (void);
 
+extern void OCTAVE_API
+gripe_modifying_nonunique (void);
+
 #endif

--- a/liboctave/util/nz-iterators.h	Sat Jun 21 13:13:05 2014 -0600
+++ b/liboctave/util/nz-iterators.h	Mon Jul 14 13:07:59 2014 -0600
@@ -245,7 +245,7 @@
   double
   interp_idx (void) const
   {
-    return to_interp_idx (row (), col (), mat.rows ());
+    return to_interp_idx (row (), col (), mat.dims ());
   }
   octave_idx_type
   col (void) const
@@ -313,7 +313,7 @@
   double
   interp_idx (void) const
   {
-    return to_interp_idx (row (), col (), mat.rows ());
+    return to_interp_idx (row (), col (), mat.dims ());
   }
   octave_idx_type
   col (void) const
@@ -369,7 +369,7 @@
   octave_idx_type
   interp_idx (void) const
   {
-    return to_interp_idx (row (), col (), mat.rows ());
+    return to_interp_idx (row (), col (), mat.dims ());
   }
   octave_idx_type
   col (void) const