--- a/liboctave/Array-C.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-C.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -96,10 +96,10 @@
 
 INSTANTIATE_ARRAY_AND_ASSIGN (Complex, OCTAVE_API);
 
-INSTANTIATE_ARRAY_ASSIGN (Complex, double, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (Complex, int, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (Complex, short, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (Complex, char, OCTAVE_API);
+INSTANTIATE_ARRAY_ASSIGN (Complex, double, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (Complex, int, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (Complex, short, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (Complex, char, OCTAVE_API)
 
 #include "Array2.h"
 

--- a/liboctave/Array-d.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-d.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -388,9 +388,9 @@
 
 INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API);
 
-INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (double, char, OCTAVE_API);
+INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (double, char, OCTAVE_API)
 
 #include "Array2.h"
 

--- a/liboctave/Array-f.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-f.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -388,9 +388,9 @@
 
 INSTANTIATE_ARRAY_AND_ASSIGN (float, OCTAVE_API);
 
-INSTANTIATE_ARRAY_ASSIGN (float, int, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (float, short, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (float, char, OCTAVE_API);
+INSTANTIATE_ARRAY_ASSIGN (float, int, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (float, short, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (float, char, OCTAVE_API)
 
 #include "Array2.h"
 

--- a/liboctave/Array-fC.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-fC.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -96,10 +96,10 @@
 
 INSTANTIATE_ARRAY_AND_ASSIGN (FloatComplex, OCTAVE_API);
 
-INSTANTIATE_ARRAY_ASSIGN (FloatComplex, float, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (FloatComplex, int, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (FloatComplex, short, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (FloatComplex, char, OCTAVE_API);
+INSTANTIATE_ARRAY_ASSIGN (FloatComplex, float, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (FloatComplex, int, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (FloatComplex, short, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (FloatComplex, char, OCTAVE_API)
 
 #include "Array2.h"
 

--- a/liboctave/Array-i.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-i.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -45,8 +45,8 @@
 INSTANTIATE_ARRAY_AND_ASSIGN (long long, OCTAVE_API);
 #endif
 
-INSTANTIATE_ARRAY_ASSIGN (int, short, OCTAVE_API);
-INSTANTIATE_ARRAY_ASSIGN (int, char, OCTAVE_API);
+INSTANTIATE_ARRAY_ASSIGN (int, short, OCTAVE_API)
+INSTANTIATE_ARRAY_ASSIGN (int, char, OCTAVE_API)
 
 INSTANTIATE_ARRAY_SORT (octave_int8);
 INSTANTIATE_ARRAY_SORT (octave_int16);

--- a/liboctave/Array-s.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-s.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -35,7 +35,7 @@
 
 INSTANTIATE_ARRAY_AND_ASSIGN (short, OCTAVE_API);
 
-INSTANTIATE_ARRAY_ASSIGN (short, char, OCTAVE_API);
+INSTANTIATE_ARRAY_ASSIGN (short, char, OCTAVE_API)
 
 #include "Array2.h"
 

--- a/liboctave/Array-util.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-util.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -322,34 +322,6 @@
 }
 
 bool
-is_in (octave_idx_type num, const idx_vector& idx)
-{
-  octave_idx_type n = idx.capacity ();
-
-  for (octave_idx_type i = 0; i < n; i++)
-    if (idx.elem (i) == num)
-      return true;
-
-  return false;
-}
-
-octave_idx_type
-how_many_lgt (const octave_idx_type num, idx_vector& idxv)
-{
-  octave_idx_type retval = 0;
-
-  octave_idx_type n = idxv.capacity ();
-
-  for (octave_idx_type i = 0; i < n; i++)
-    {
-      if (num > idxv.elem (i))
-	retval++;
-    }
-
-  return retval;
-}
-
-bool
 all_ones (const Array<octave_idx_type>& arr)
 {
   bool retval = true;
@@ -413,62 +385,90 @@
   return retval;
 }
 
-dim_vector
-short_freeze (Array<idx_vector>& ra_idx, const dim_vector& dimensions,
-	      int resize_ok)
+dim_vector zero_dims_inquire (const Array<idx_vector>& ia,
+                              const dim_vector& rhdv)
 {
-  dim_vector retval;
-
-  int n = ra_idx.length ();
-
-  int n_dims = dimensions.length ();
-
-  if (n == n_dims)
-    {
-      retval = freeze (ra_idx, dimensions, resize_ok);
-    }
-  else if (n < n_dims)
+  int ial = ia.length (), rhdvl = rhdv.length ();
+  dim_vector rdv;
+  rdv.resize (ial);
+  bool *scalar = new bool[ial], *colon = new bool[ial];
+  // Mark scalars and colons, count non-scalar indices.
+  int nonsc = 0; 
+  bool all_colons = true;
+  for (int i = 0; i < ial; i++)
     {
-      retval.resize (n);
-      
-      for (int i = 0; i < n - 1; i++)
-        retval(i) = ra_idx(i).freeze (dimensions(i), "dimension", resize_ok);
-
-      int size_left = 1;
-
-      for (int i = n - 1; i < n_dims; i++)
-        size_left *= dimensions(i); 
- 
-      if (ra_idx(n-1).is_colon())
-        {
-	  retval(n-1) = size_left;
-	}
-      else
-	{
-	  octave_idx_type last_ra_idx = ra_idx(n-1)(0);
-	  for (octave_idx_type i = 1; i < ra_idx (n - 1).capacity (); i++)
-	    last_ra_idx = (last_ra_idx > ra_idx(n-1)(i) ? last_ra_idx : 
-			   ra_idx(n-1)(i));
-
-	  if (last_ra_idx < size_left)
-            {
-              retval(n-1) = ra_idx(n-1).freeze (size_left,
-						"dimension", resize_ok);
-            }
-          else
-            {
-	      // Make it larger than it should be to get an error
-	      // later.
-
-	      retval.resize (n_dims+1);
-
-	      (*current_liboctave_error_handler)
-		("index exceeds N-d array dimensions");
-	    }
-	}
+      // FIXME: should we check for length() instead?
+      scalar[i] = ia(i).is_scalar ();
+      colon[i] = ia(i).is_colon ();
+      if (! scalar[i]) nonsc++;
+      if (! colon[i]) rdv(i) = ia(i).extent (0);
+      all_colons = all_colons && colon[i];
     }
 
-  return retval;
+  bool match = false;
+  // If the number of nonscalar indices matches the dimensionality of RHS,
+  // we try an exact match, inquiring even singleton dimensions.
+  if (all_colons)
+    {
+      rdv = rhdv;
+      rdv.resize(ial, 1);
+    }
+  else if (nonsc == rhdvl)
+    {
+      for (int i = 0, j = 0; i < ial; i++)
+        {
+          if (scalar[i]) continue;
+          if (colon[i])
+            rdv(i) = rhdv(j++);
+        }
+    }
+  else
+    {
+      dim_vector rhdv0 = rhdv;
+      rhdv0.chop_all_singletons ();
+      int rhdv0l = rhdv0.length ();
+      for (int i = 0, j = 0; i < ial; i++)
+        {
+          if (scalar[i]) continue;
+          if (colon[i])
+            rdv(i) =  (j < rhdv0l) ? rhdv0(j++) : 1;
+        }
+    }
+
+  delete [] scalar;
+  delete [] colon;
+
+  return rdv;
+}
+
+dim_vector zero_dims_inquire (const idx_vector& i, const idx_vector& j,
+                              const dim_vector& rhdv)
+{
+  bool icol = i.is_colon (), jcol = j.is_colon ();
+  dim_vector rdv;
+  if (icol && jcol && rhdv.length () == 2)
+    {
+      rdv(0) = rhdv(0);
+      rdv(1) = rhdv(1);
+    }
+  else
+    {
+      dim_vector rhdv0 = rhdv;
+      rhdv0.chop_all_singletons ();
+      int k = 0;
+      rdv(0) = i.extent (0);
+      if (icol)
+        rdv(0) = rhdv0(k++);
+      else if (! i.is_scalar ())
+        k++;
+      rdv(1) = j.extent (0);
+      if (jcol)
+        rdv(1) = rhdv0(k++);
+      else if (! j.is_scalar ())
+        k++;
+    }
+
+  return rdv;
 }
 
 int

--- a/liboctave/Array-util.h	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array-util.h	Mon Oct 20 16:54:28 2008 +0200
@@ -65,10 +65,6 @@
 extern OCTAVE_API bool all_colon_equiv (const Array<idx_vector>& ra_idx,
 			     const dim_vector& frozen_lengths);
 
-extern OCTAVE_API bool is_in (octave_idx_type num, const idx_vector& idx);
-
-extern OCTAVE_API octave_idx_type how_many_lgt (const octave_idx_type num, idx_vector& idxv);
-
 extern OCTAVE_API bool all_ones (const Array<octave_idx_type>& arr);
 
 extern OCTAVE_API Array<octave_idx_type> get_elt_idx (const Array<idx_vector>& ra_idx,
@@ -76,9 +72,11 @@
 
 extern OCTAVE_API Array<octave_idx_type> get_ra_idx (octave_idx_type idx, const dim_vector& dims);
 
-extern OCTAVE_API dim_vector short_freeze (Array<idx_vector>& ra_idx,
-				const dim_vector& dimensions,
-				int resize_ok);
+extern OCTAVE_API dim_vector zero_dims_inquire (const Array<idx_vector>& ia,
+                                                const dim_vector& rhdv);
+
+extern OCTAVE_API dim_vector zero_dims_inquire (const idx_vector& i, const idx_vector& j,
+                                                const dim_vector& rhdv);
 
 struct
 permute_vector

--- a/liboctave/Array.cc	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array.cc	Mon Oct 20 16:54:28 2008 +0200
@@ -3,6 +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@gmail.com>
 
 This file is part of Octave.
 
@@ -32,6 +33,7 @@
 #include <iostream>
 #include <sstream>
 #include <vector>
+#include <algorithm>
 #include <new>
 
 #include "Array.h"
@@ -44,7 +46,7 @@
 
 template <class T>
 Array<T>::Array (const Array<T>& a, const dim_vector& dv)
-  : rep (a.rep), dimensions (dv), idx (0), idx_count (0)
+  : rep (a.rep), dimensions (dv)
 {
   rep->count++;
 
@@ -58,8 +60,6 @@
 {
   if (--rep->count <= 0)
     delete rep;
-
-  delete [] idx;
 }
 
 template <class T>
@@ -75,10 +75,6 @@
       rep->count++;
 
       dimensions = a.dimensions;
-
-      delete [] idx;
-      idx_count = 0;
-      idx = 0;
     }
 
   return *this;
@@ -559,458 +555,896 @@
   return retval;
 }
 
-template <class T>
-void
-Array<T>::resize_no_fill (octave_idx_type n)
+// Helper class for multi-d index reduction and recursive indexing/indexed assignment.
+// Rationale: we could avoid recursion using a state machine instead. However, using
+// recursion is much more amenable to possible parallelization in the future.
+// Also, the recursion solution is cleaner and more understandable.
+class rec_index_helper
 {
-  if (n < 0)
+  octave_idx_type *dim, *cdim;
+  idx_vector *idx;
+  int top;
+
+public:
+  rec_index_helper (const dim_vector& dv, const Array<idx_vector>& ia)
     {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return;
+      int n = ia.length ();
+      assert (n > 0 && (dv.length () == std::max (n, 2)));
+
+      dim = new octave_idx_type [2*n];
+      // A hack to avoid double allocation
+      cdim = dim + n;
+      idx = new idx_vector [n];
+      top = 0;
+
+      dim[0] = dv(0);
+      cdim[0] = 1;
+      idx[0] = ia(0);
+
+      for (int i = 1; i < n; i++)
+        {
+          // Try reduction...
+          if (idx[top].maybe_reduce (dim[top], ia(i), dv(i)))
+            {
+              // Reduction successful, fold dimensions.
+              dim[top] *= dv(i);
+            }
+          else
+            {
+              // Unsuccessful, store index & cumulative dim.
+              top++;
+              idx[top] = ia(i);
+              dim[top] = dv(i);
+              cdim[top] = cdim[top-1] * dim[top-1];
+            } 
+        }
+    }
+
+  ~rec_index_helper (void) { delete [] idx; delete [] dim; }
+
+private:
+
+  // Recursive N-d indexing
+  template <class T>
+  T *do_index (const T *src, T *dest, int lev) const
+    {
+      if (lev == 0)
+        dest += idx[0].index (src, dim[0], dest);
+      else
+        {
+          octave_idx_type n = idx[lev].length (dim[lev]), d = cdim[lev];
+          for (octave_idx_type i = 0; i < n; i++)
+            dest = do_index (src + d*idx[lev].xelem (i), dest, lev-1);
+        }
+
+      return dest;
+    }
+  
+  // Recursive N-d indexed assignment
+  template <class T>
+  const T *do_assign (const T *src, T *dest, int lev) const
+    {
+      if (lev == 0)
+        src += idx[0].assign (src, dim[0], dest);
+      else
+        {
+          octave_idx_type n = idx[lev].length (dim[lev]), d = cdim[lev];
+          for (octave_idx_type i = 0; i < n; i++)
+            src = do_assign (src, dest + d*idx[lev].xelem (i), lev-1);
+        }
+
+      return src;
     }
 
-  if (n == length ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = rep;
-  const T *old_data = data ();
-  octave_idx_type old_len = length ();
-
-  rep = new typename Array<T>::ArrayRep (n);
-
-  dimensions = dim_vector (n);
-
-  if (n > 0 && old_data && old_len > 0)
+  // Recursive N-d indexed assignment
+  template <class T>
+  void do_fill (const T& val, T *dest, int lev) const
     {
-      octave_idx_type min_len = old_len < n ? old_len : n;
-
-      for (octave_idx_type i = 0; i < min_len; i++)
-	xelem (i) = old_data[i];
+      if (lev == 0)
+        idx[0].fill (val, dim[0], dest);
+      else
+        {
+          octave_idx_type n = idx[lev].length (dim[lev]), d = cdim[lev];
+          for (octave_idx_type i = 0; i < n; i++)
+            do_fill (val, dest + d*idx[lev].xelem (i), lev-1);
+        }
     }
 
-  if (--old_rep->count <= 0)
-    delete old_rep;
+public:
+
+  template <class T>
+  void index (const T *src, T *dest) const { do_index (src, dest, top); }
+
+  template <class T>
+  void assign (const T *src, T *dest) const { do_assign (src, dest, top); }
+
+  template <class T>
+  void fill (const T& val, T *dest) const { do_fill (val, dest, top); }
+
+};
+
+// Helper class for multi-d recursive resizing
+// This handles resize () in an efficient manner, touching memory only
+// once (apart from reinitialization)
+class rec_resize_helper
+{
+  octave_idx_type *cext, *sext, *dext;
+  int n;
+
+public:
+  rec_resize_helper (const dim_vector& ndv, const dim_vector& odv)
+    {
+      int l = ndv.length ();
+      assert (odv.length () == l);
+      octave_idx_type ld = 1;
+      int i = 0;
+      for (; i < l-1 && ndv(i) == odv(i); i++) ld *= ndv(i);
+      n = l - i;
+      cext = new octave_idx_type[3*n];
+      // Trick to avoid three allocations
+      sext = cext + n;
+      dext = sext + n;
+
+      octave_idx_type sld = ld, dld = ld;
+      for (int j = 0; j < n; j++)
+        {
+          cext[j] = std::min (ndv(i+j), odv(i+j));
+          sext[j] = sld *= odv(i+j);
+          dext[j] = dld *= ndv(i+j);
+        }
+      cext[0] *= ld;
+    }
+
+  ~rec_resize_helper (void) { delete [] cext; }
+
+private:
+  // recursive resizing
+  template <class T>
+  void do_resize_fill (const T* src, T *dest, const T& rfv, int lev) const
+    {
+      if (lev == 0)
+        {
+          T* destc = std::copy (src, src + cext[0], dest);
+          std::fill (destc, dest + dext[0], rfv);
+        }
+      else
+        {
+          octave_idx_type sd = sext[lev-1], dd = dext[lev-1], k;
+          for (k = 0; k < cext[lev]; k++)
+            do_resize_fill (src + k * sd, dest + k * dd, rfv, lev - 1);
+
+          std::fill (dest + k * dd, dest + dext[lev], rfv);
+        }
+    }
+public:
+  template <class T>
+  void resize_fill (const T* src, T *dest, const T& rfv) const 
+    { do_resize_fill (src, dest, rfv, n-1); }
+
+};
+
+static void gripe_index_out_of_range (void)
+{
+  (*current_liboctave_error_handler)
+    ("A(I): Index exceeds matrix dimension.");
 }
 
 template <class T>
-void
-Array<T>::resize_no_fill (const dim_vector& dv)
+Array<T>
+Array<T>::index (const idx_vector& i) const
 {
-  octave_idx_type n = dv.length ();
+  octave_idx_type n = numel ();
+  Array<T> retval;
 
-  for (octave_idx_type i = 0; i < n; i++)
+  if (i.is_colon ())
     {
-      if (dv(i) < 0)
-	{
-	  (*current_liboctave_error_handler)
-	    ("can't resize to negative dimension");
-	  return;
-	}
+      // A(:) produces a shallow copy as a column vector.
+      retval.dimensions = dim_vector (n, 1);
+      rep->count++;
+      retval.rep = rep;
     }
-
-  bool same_size = true;
-
-  if (dimensions.length () != n)
+  else if (i.extent (n) != n)
     {
-      same_size = false;
+      gripe_index_out_of_range ();
     }
   else
     {
-      for (octave_idx_type i = 0; i < n; i++)
-	{
-	  if (dv(i) != dimensions(i))
-	    {
-	      same_size = false;
-	      break;
-	    }
-	}
+      // FIXME: This is the only place where orig_dimensions are used.
+      dim_vector rd = i.orig_dimensions ();
+      octave_idx_type il = i.length (n);
+
+      // FIXME: This is for Matlab compatibility. Matlab 2007 given `b = ones(3,1)'
+      // yields the following:
+      // b(zeros(0,0)) gives []
+      // b(zeros(1,0)) gives zeros(0,1)
+      // b(zeros(0,1)) gives zeros(0,1)
+      // b(zeros(0,m)) gives zeros(0,m)
+      // b(zeros(m,0)) gives zeros(m,0)
+      // b(1:2) gives ones(2,1)
+      // b(ones(2)) gives ones(2) etc.
+      // As you can see, the behaviour is weird, but the tests end up pretty
+      // simple. Nah, I don't want to suggest that this is ad hoc :) Neither do
+      // I want to say that Matlab is a lousy piece of s...oftware.
+      if (ndims () == 2 && n != 1)
+        {
+          if (columns () == 1 && rd(0) == 1)
+            rd = dim_vector (il, 1);
+          else if (rows () == 1 && rd(1) == 1)
+            rd = dim_vector (1, il);
+        }
+
+      // 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 ());
+    }
+
+  return retval;
+}
+
+template <class T>
+Array<T>
+Array<T>::index (const idx_vector& i, const idx_vector& j) const
+{
+  // Get dimensions, allowing Fortran indexing in the 2nd dim.
+  dim_vector dv = dimensions.redim (2);
+  octave_idx_type r = dv(0), c = dv(1);
+  Array<T> retval;
+
+  if (i.is_colon () && j.is_colon ())
+    {
+      // A(:,:) produces a shallow copy.
+      retval = Array<T> (*this, dv);
+    }
+  else if (i.extent (r) != r || j.extent (c) != c)
+    {
+      gripe_index_out_of_range ();
+    }
+  else
+    {
+      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))
+        ii.index (src, n, dest);
+      else
+        {
+          for (octave_idx_type k = 0; k < jl; k++)
+            dest += i.index (src + r * j.xelem (k), r, dest);
+        }
     }
 
-  if (same_size)
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = rep;
-  const T *old_data = data ();
-
-  octave_idx_type ts = get_size (dv);
+  return retval;
+}
 
-  rep = new typename Array<T>::ArrayRep (ts);
+template <class T>
+Array<T>
+Array<T>::index (const Array<idx_vector>& ia) const
+{
+  int ial = ia.length ();
+  Array<T> retval;
 
-  dim_vector dv_old = dimensions;
-  octave_idx_type  dv_old_orig_len = dv_old.length ();
-  dimensions = dv;
-  octave_idx_type ts_old = get_size (dv_old);
-
-  if (ts > 0 && ts_old > 0 && dv_old_orig_len > 0)
+  // FIXME: is this dispatching necessary?
+  if (ial == 1)
+    retval = index (ia(0));
+  else if (ial == 2)
+    retval = index (ia(0), ia(1));
+  else if (ial > 0)
     {
-      Array<octave_idx_type> ra_idx (dimensions.length (), 0);
+      // Get dimensions, allowing Fortran indexing in the last dim.
+      dim_vector dv = dimensions.redim (ial);
 
-      if (n > dv_old_orig_len)
-	{
-	  dv_old.resize (n);
+      // Check for out of bounds conditions.
+      bool all_colons = true, mismatch = false;
+      for (int i = 0; i < ial; i++)
+        {
+          if (ia(i).extent (dv(i)) != dv(i))
+            {
+              mismatch = true;
+              break;
+            }
+          else
+            all_colons = all_colons && ia(i).is_colon ();
+        }
+
 
-	  for (octave_idx_type i = dv_old_orig_len; i < n; i++)
-	    dv_old.elem (i) = 1;
-	}
+      if (mismatch)
+        {
+          gripe_index_out_of_range ();
+        }
+      else if (all_colons)
+        {
+          // A(:,:,...,:) produces a shallow copy.
+          retval = Array<T> (*this, dv);
+        }
+      else 
+        {
+          // Form result dimensions.
+          dim_vector rdv;
+          rdv.resize (ial);
+          for (int i = 0; i < ial; i++) rdv(i) = ia(i).length (dv(i));
+          rdv.chop_trailing_singletons ();
 
-      for (octave_idx_type i = 0; i < ts; i++)
-	{
-	  if (index_in_bounds (ra_idx, dv_old))
-	    rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)];
+          // Don't use resize here to avoid useless initialization for POD types.
+          retval = Array<T> (rdv);
 
-	  increment_index (ra_idx, dimensions);
-	}
+          // Prepare for recursive indexing
+          rec_index_helper rh (dv, ia);
+
+          // Do it.
+          rh.index (data (), retval.fortran_vec ());
+        }
     }
 
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  return retval;
 }
 
+// FIXME: the following is a common error message to resize, regardless of whether it's
+// called from assign or elsewhere. It seems OK to me, but eventually the gripe can be
+// specialized. Anyway, propagating various error messages into procedure is, IMHO, a
+// nonsense. If anything, we should change error handling here (and throughout liboctave)
+// to allow custom handling of errors
+static void gripe_invalid_resize (void)
+{
+  (*current_liboctave_error_handler)
+    ("resize: Invalid resizing operation or ambiguous assignment to an out-of-bounds array element.");
+}
+
+// The default fill value. Override if you want a different one.
+
+template <class T>
+T Array<T>::resize_fill_value ()
+{
+  return T ();
+}
+
+// Yes, we could do resize using index & assign.  However, that would possibly
+// involve a lot more memory traffic than we actually need.
+
 template <class T>
 void
-Array<T>::resize_no_fill (octave_idx_type r, octave_idx_type c)
+Array<T>::resize_fill (octave_idx_type n, const T& rfv)
 {
-  if (r < 0 || c < 0)
+  if (n >= 0 && ndims () == 2)
     {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return;
-    }
-
-  int n = ndims ();
-
-  if (n == 0)
-    dimensions = dim_vector (0, 0);
-
-  assert (ndims () == 2);
-
-  if (r == dim1 () && c == dim2 ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
-  const T *old_data = data ();
-
-  octave_idx_type old_d1 = dim1 ();
-  octave_idx_type old_d2 = dim2 ();
-  octave_idx_type old_len = length ();
-
-  octave_idx_type ts = get_size (r, c);
-
-  rep = new typename Array<T>::ArrayRep (ts);
-
-  dimensions = dim_vector (r, c);
-
-  if (ts > 0 && old_data && old_len > 0)
-    {
-      octave_idx_type min_r = old_d1 < r ? old_d1 : r;
-      octave_idx_type min_c = old_d2 < c ? old_d2 : c;
-
-      for (octave_idx_type j = 0; j < min_c; j++)
-	for (octave_idx_type i = 0; i < min_r; i++)
-	  xelem (i, j) = old_data[old_d1*j+i];
-    }
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
-}
+      dim_vector dv;
+      // This is driven by Matlab's behaviour of giving a *row* vector on
+      // some out-of-bounds assignments. Specifically, matlab allows a(i) with
+      // out-of-bouds i when a is either of 0x0, 1x0, 1x1, 0xN, and gives
+      // a row vector in all cases (yes, even the last one, search me why).
+      // Giving a column vector would make much more sense (given the way
+      // trailing singleton dims are treated), but hey, Matlab is not here to
+      // make sense, it's here to make money ;)
+      bool invalid = false;
+      if (rows () == 0 || rows () == 1)
+        dv = dim_vector (1, n);          
+      else if (columns () == 1)
+        dv = dim_vector (n, 1);
+      else
+         invalid = true;
+        
+      if (invalid)
+        gripe_invalid_resize ();
+      else
+        {
+          octave_idx_type nx = numel ();
+          if (n != nx)
+            {
+              Array<T> tmp = Array<T> (dv);
+              T *dest = tmp.fortran_vec ();
 
-template <class T>
-void
-Array<T>::resize_no_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p)
-{
-  if (r < 0 || c < 0 || p < 0)
-    {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return;
-    }
-
-  int n = ndims ();
-
-  if (n == 0)
-    dimensions = dim_vector (0, 0, 0);
-
-  assert (ndims () == 3);
-
-  if (r == dim1 () && c == dim2 () && p == dim3 ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = rep;
-  const T *old_data = data ();
+              octave_idx_type n0 = std::min (n, nx), n1 = n - n0;
+              dest = std::copy (data (), data () + n0, dest);
+              std::fill (dest, dest + n1, rfv);
 
-  octave_idx_type old_d1 = dim1 ();
-  octave_idx_type old_d2 = dim2 ();
-  octave_idx_type old_d3 = dim3 ();
-  octave_idx_type old_len = length ();
-
-  octave_idx_type ts = get_size (get_size (r, c), p);
-
-  rep = new typename Array<T>::ArrayRep (ts);
-
-  dimensions = dim_vector (r, c, p);
-
-  if (ts > 0 && old_data && old_len > 0)
-    {
-      octave_idx_type min_r = old_d1 < r ? old_d1 : r;
-      octave_idx_type min_c = old_d2 < c ? old_d2 : c;
-      octave_idx_type min_p = old_d3 < p ? old_d3 : p;
-
-      for (octave_idx_type k = 0; k < min_p; k++)
-	for (octave_idx_type j = 0; j < min_c; j++)
-	  for (octave_idx_type i = 0; i < min_r; i++)
-	    xelem (i, j, k) = old_data[old_d1*(old_d2*k+j)+i];
+              *this = tmp;
+            }
+        }
     }
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  else
+    gripe_invalid_resize ();
 }
 
 template <class T>
 void
-Array<T>::resize_and_fill (octave_idx_type n, const T& val)
+Array<T>::resize_fill (octave_idx_type r, octave_idx_type c, const T& rfv)
 {
-  if (n < 0)
+  if (r >= 0 && c >= 0 && ndims () == 2)
     {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return;
+      octave_idx_type rx = rows (), cx = columns ();
+      if (r != rx || c != cx)
+        {
+          Array<T> tmp = Array<T> (dim_vector (r, c));
+          T *dest = tmp.fortran_vec ();
+
+          octave_idx_type r0 = std::min (r, rx), r1 = r - r0;
+          octave_idx_type c0 = std::min (c, cx), c1 = c - c0;
+          const T *src = data ();
+          if (r == rx)
+            dest = std::copy (src, src + r * c0, dest);
+          else
+            {
+              for (octave_idx_type k = 0; k < c0; k++)
+                {
+                  dest = std::copy (src, src + r0, dest);
+                  src += rx;
+                  std::fill (dest, dest + r1, rfv);
+                  dest += r1;
+                }
+            }
+
+          std::fill (dest, dest + r * c1, rfv);
+
+          *this = tmp;
+        }
+    }
+  else
+    gripe_invalid_resize ();
+
+}
+
+template<class T>
+void
+Array<T>::resize_fill (const dim_vector& dv, const T& rfv)
+{
+  int dvl = dv.length ();
+  if (dvl == 1)
+    resize (dv(0), rfv);
+  else if (dvl == 2)
+    resize (dv(0), dv(1), rfv);
+  else if (dimensions != dv)
+    {
+      if (dimensions.length () <= dvl)
+        {
+          Array<T> tmp (dv);
+          // Prepare for recursive resizing.
+          rec_resize_helper rh (dv, dimensions.redim (dvl));
+
+          // Do it.
+          rh.resize_fill (data (), tmp.fortran_vec (), rfv);   
+          *this = tmp;
+        }
+      else
+        gripe_invalid_resize ();
+    }
+}
+
+template <class T>
+Array<T> 
+Array<T>::index (const idx_vector& i, bool resize_ok, const T& rfv) const
+{
+  Array<T> tmp = *this;
+  if (resize_ok)
+    {
+      octave_idx_type n = numel (), nx = i.extent (n);
+      if (n != nx)
+        tmp.resize_fill (nx, rfv);
+
+      if (tmp.numel () != nx)
+        return Array<T> ();
     }
 
-  if (n == length ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = rep;
-  const T *old_data = data ();
-  octave_idx_type old_len = length ();
-
-  rep = new typename Array<T>::ArrayRep (n);
-
-  dimensions = dim_vector (n);
+  return tmp.index (i);
+}
 
-  if (n > 0)
+template <class T>
+Array<T> 
+Array<T>::index (const idx_vector& i, const idx_vector& j, 
+                 bool resize_ok, const T& rfv) const
+{
+  Array<T> tmp = *this;
+  if (resize_ok)
     {
-      octave_idx_type min_len = old_len < n ? old_len : n;
+      dim_vector dv = dimensions.redim (2);
+      octave_idx_type r = dv(0), c = dv(1);
+      octave_idx_type rx = i.extent (r), cx = j.extent (c);
+      if (r != rx || c != cx)
+        tmp.resize_fill (rx, cx, rfv);
 
-      if (old_data && old_len > 0)
-	{
-	  for (octave_idx_type i = 0; i < min_len; i++)
-	    xelem (i) = old_data[i];
-	}
-
-      for (octave_idx_type i = old_len; i < n; i++)
-	xelem (i) = val;
+      if (tmp.rows () != rx || tmp.columns () != cx)
+        return Array<T> ();
     }
 
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  return tmp.index (i, j);  
+}
+
+template <class T>
+Array<T> 
+Array<T>::index (const Array<idx_vector>& ia,
+                 bool resize_ok, const T& rfv) const
+{
+  Array<T> tmp = *this;
+  if (resize_ok)
+    {
+      int ial = ia.length ();
+      dim_vector dv = dimensions.redim (ial);
+      dim_vector dvx; dvx.resize (ial);
+      for (int i = 0; i < ial; i++) dvx(i) = ia(i).extent (dv (i));
+      if (! (dvx == dv))
+        tmp.resize_fill (dvx, rfv);
+
+      if (tmp.dimensions != dvx)
+        return Array<T> ();
+    }
+
+  return tmp.index (ia);  
+}
+
+
+static void 
+gripe_invalid_assignment_size (void)
+{
+  (*current_liboctave_error_handler)
+    ("A(I) = X: X must have the same size as I");
+}
+
+static void
+gripe_assignment_dimension_mismatch (void)
+{
+  (*current_liboctave_error_handler)
+    ("A(I,J,...) = X: dimensions mismatch");
 }
 
 template <class T>
 void
-Array<T>::resize_and_fill (octave_idx_type r, octave_idx_type c, const T& val)
+Array<T>::assign (const idx_vector& i, const Array<T>& rhs, const T& rfv)
 {
-  if (r < 0 || c < 0)
-    {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return;
-    }
-
-  if (ndims () == 0)
-    dimensions = dim_vector (0, 0);
+  octave_idx_type n = numel (), rhl = rhs.numel ();
 
-  assert (ndims () == 2);
-
-  if (r == dim1 () && c == dim2 ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
-  const T *old_data = data ();
-
-  octave_idx_type old_d1 = dim1 ();
-  octave_idx_type old_d2 = dim2 ();
-  octave_idx_type old_len = length ();
-
-  octave_idx_type ts = get_size (r, c);
+  if (rhl == 1 || i.length (n) == rhl)
+    {
+      octave_idx_type nx = i.extent (n);
+      // Try to resize first if necessary. 
+      if (nx != n)
+        {
+          resize_fill (nx, rfv);      
+          n = numel ();
+        }
 
-  rep = new typename Array<T>::ArrayRep (ts);
-
-  dimensions = dim_vector (r, c);
-
-  if (ts > 0)
-    {
-      octave_idx_type min_r = old_d1 < r ? old_d1 : r;
-      octave_idx_type min_c = old_d2 < c ? old_d2 : c;
-
-      if (old_data && old_len > 0)
-	{
-	  for (octave_idx_type j = 0; j < min_c; j++)
-	    for (octave_idx_type i = 0; i < min_r; i++)
-	      xelem (i, j) = old_data[old_d1*j+i];
-	}
-
-      for (octave_idx_type j = 0; j < min_c; j++)
-	for (octave_idx_type i = min_r; i < r; i++)
-	  xelem (i, j) = val;
-
-      for (octave_idx_type j = min_c; j < c; j++)
-	for (octave_idx_type i = 0; i < r; i++)
-	  xelem (i, j) = val;
+      // If the resizing was ambiguous, resize has already griped.
+      if (nx == n)
+        {
+          if (i.is_colon ())
+            {
+              // A(:) = X makes a full fill or a shallow copy.
+              if (rhl == 1)
+                fill (rhs(0));
+              else
+                *this = rhs.reshape (dimensions);
+            }
+          else
+            {
+              if (rhl == 1)
+                i.fill (rhs(0), n, fortran_vec ());
+              else
+                i.assign (rhs.data (), n, fortran_vec ());
+            }
+        }
     }
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  else
+    gripe_invalid_assignment_size ();
 }
 
 template <class T>
 void
-Array<T>::resize_and_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val)
+Array<T>::assign (const idx_vector& i, const idx_vector& j,
+                  const Array<T>& rhs, const T& rfv)
 {
-  if (r < 0 || c < 0 || p < 0)
+  // Get RHS extents, discarding singletons.
+  dim_vector rhdv = rhs.dims (); 
+  // Get LHS extents, allowing Fortran indexing in the second dim.
+  dim_vector dv = dimensions.redim (2);
+  // Check for out-of-bounds and form resizing dimensions.
+  dim_vector rdv; 
+  // In the special when all dimensions are zero, colons are allowed to inquire
+  // the shape of RHS. The rules are more obscure, so we solve that elsewhere.
+  if (dv.all_zero ())
+    rdv = zero_dims_inquire (i, j, rhdv);
+  else
     {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return;
+      rdv(0) = i.extent (dv(0));
+      rdv(1) = j.extent (dv(1));
     }
 
-  if (ndims () == 0)
-    dimensions = dim_vector (0, 0, 0);
-
-  assert (ndims () == 3);
+  bool isfill = rhs.numel () == 1;
+  octave_idx_type il = i.length (rdv(0)), jl = j.length (rdv(1));
+  rhdv.chop_all_singletons ();
+  bool match = isfill || (rhdv.length () == 2 && il == rhdv(0) && jl == rhdv(1));
+  match = match || (il == 1 && jl == rhdv(0) && rhdv(1) == 1);
 
-  if (r == dim1 () && c == dim2 () && p == dim3 ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = rep;
-  const T *old_data = data ();
+  if (match)
+    {
+      // Resize if requested.
+      if (rdv != dv)
+        {
+          resize (rdv, rfv);
+          dv = dimensions;
+        }
 
-  octave_idx_type old_d1 = dim1 ();
-  octave_idx_type old_d2 = dim2 ();
-  octave_idx_type old_d3 = dim3 ();
-
-  octave_idx_type old_len = length ();
-
-  octave_idx_type ts = get_size (get_size (r, c), p);
-
-  rep = new typename Array<T>::ArrayRep (ts);
-
-  dimensions = dim_vector (r, c, p);
-
-  if (ts > 0)
-    {
-      octave_idx_type min_r = old_d1 < r ? old_d1 : r;
-      octave_idx_type min_c = old_d2 < c ? old_d2 : c;
-      octave_idx_type min_p = old_d3 < p ? old_d3 : p;
+      // If the resizing was invalid, resize has already griped.
+      if (dv == rdv)
+        {
+          if (i.is_colon () && j.is_colon ())
+            {
+              // A(:,:) = X makes a full fill or a shallow copy
+              if (isfill)
+                fill (rhs(0));
+              else
+                *this = rhs.reshape (dimensions);
+            }
+          else
+            {
+              // The actual work.
+              octave_idx_type n = numel (), r = rows (), c = columns ();
+              idx_vector ii (i);
 
-      if (old_data && old_len > 0)
-	for (octave_idx_type k = 0; k < min_p; k++)
-	  for (octave_idx_type j = 0; j < min_c; j++)
-	    for (octave_idx_type i = 0; i < min_r; i++)
-	      xelem (i, j, k) = old_data[old_d1*(old_d2*k+j)+i];
-
-      // FIXME -- if the copy constructor is expensive, this
-      // may win.  Otherwise, it may make more sense to just copy the
-      // value everywhere when making the new ArrayRep.
+              const T* src = rhs.data ();
+              T *dest = fortran_vec ();
 
-      for (octave_idx_type k = 0; k < min_p; k++)
-	for (octave_idx_type j = min_c; j < c; j++)
-	  for (octave_idx_type i = 0; i < min_r; i++)
-	    xelem (i, j, k) = val;
+              // Try reduction first.
+              if (ii.maybe_reduce (r, j, c))
+                {
+                  if (isfill)
+                    ii.fill (*src, n, dest);
+                  else
+                    ii.assign (src, n, dest);
+                }
+              else
+                {
+                  if (isfill)
+                    {
+                      for (octave_idx_type k = 0; k < jl; k++)
+                        i.fill (*src, r, dest + r * j.xelem (k));
+                    }
+                  else
+                    {
+                      for (octave_idx_type k = 0; k < jl; k++)
+                        src += i.assign (src, r, dest + r * j.xelem (k));
+                    }
+                }
+            }
+          
+        }
 
-      for (octave_idx_type k = 0; k < min_p; k++)
-	for (octave_idx_type j = 0; j < c; j++)
-	  for (octave_idx_type i = min_r; i < r; i++)
-	    xelem (i, j, k) = val;
-
-      for (octave_idx_type k = min_p; k < p; k++)
-	for (octave_idx_type j = 0; j < c; j++)
-	  for (octave_idx_type i = 0; i < r; i++)
-	    xelem (i, j, k) = val;
     }
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  else
+    gripe_assignment_dimension_mismatch ();
 }
 
 template <class T>
 void
-Array<T>::resize_and_fill (const dim_vector& dv, const T& val)
+Array<T>::assign (const Array<idx_vector>& ia,
+                  const Array<T>& rhs, const T& rfv)
 {
-  octave_idx_type n = dv.length ();
+  int ial = ia.length ();
 
-  for (octave_idx_type i = 0; i < n; i++)
+  // FIXME: is this dispatching necessary / desirable?
+  if (ial == 1)
+    assign (ia(0), rhs, rfv);
+  else if (ial == 2)
+    assign (ia(0), ia(1), rhs, rfv);
+  else if (ial > 0)
     {
-      if (dv(i) < 0)
-	{
-	  (*current_liboctave_error_handler)
-	    ("can't resize to negative dimension");
-	  return;
-	}
-    }
+      // Get RHS extents, discarding singletons.
+      dim_vector rhdv = rhs.dims ();
+
+      // Get LHS extents, allowing Fortran indexing in the second dim.
+      dim_vector dv = dimensions.redim (ial);
+      
+      // Get the extents forced by indexing. 
+      dim_vector rdv;
+
+      // In the special when all dimensions are zero, colons are allowed to
+      // inquire the shape of RHS. The rules are more obscure, so we solve that
+      // elsewhere.
+      if (dv.all_zero ())
+        rdv = zero_dims_inquire (ia, rhdv);
+      else
+        {
+          rdv.resize (ial);
+          for (int i = 0; i < ial; i++)
+            rdv(i) = ia(i).extent (dv(i));
+        }
+
+      // Check whether LHS and RHS match, up to singleton dims.
+      bool match = true, all_colons = true, isfill = rhs.numel () == 1;
+
+      rhdv.chop_all_singletons ();
+      int j = 0, rhdvl = rhdv.length ();
+      for (int i = 0; i < ial; i++)
+        {
+          all_colons = all_colons && ia(i).is_colon ();
+          octave_idx_type l = ia(i).length (rdv(i));
+          if (l == 1) continue;
+          match = match && j < rhdvl && l == rhdv(j++);
+        }
+
+      match = match && (j == rhdvl || rhdv(j) == 1);
+      match = match || isfill;
+            
+      if (match)
+        {
+          // Resize first if necessary.
+          if (rdv != dv)
+            {
+              resize_fill (rdv, rfv);
+              dv = dimensions;
+              chop_trailing_singletons ();
+            }
 
-  bool same_size = true;
+          // If the resizing was invalid, resize has already griped.
+          if (dv == rdv)
+            {
+              if (all_colons)
+                {
+                  // A(:,:,...,:) = X makes a full fill or a shallow copy.
+                  if (isfill)
+                    fill (rhs(0));
+                  else
+                    *this = rhs.reshape (dimensions);
+                }
+              else
+                {
+                  // Do the actual work.
+
+                  // Prepare for recursive indexing
+                  rec_index_helper rh (dv, ia);
 
-  if (dimensions.length () != n)
-    {
-      same_size = false;
+                  // Do it.
+                  if (isfill)
+                    rh.fill (rhs(0), fortran_vec ());
+                  else
+                    rh.assign (rhs.data (), fortran_vec ());
+                }
+            }
+        }
+      else 
+        gripe_assignment_dimension_mismatch ();
     }
-  else
+}
+
+template <class T>
+void 
+Array<T>::delete_elements (const idx_vector& i)
+{
+  octave_idx_type n = numel ();
+  if (i.is_colon ())
+    { 
+      *this = Array<T> ();
+    }
+  else if (i.extent (n) != n)
     {
-      for (octave_idx_type i = 0; i < n; i++)
-	{
-	  if (dv(i) != dimensions(i))
-	    {
-	      same_size = false;
-	      break;
-	    }
-	}
+      gripe_index_out_of_range ();
+    }
+  else if (i.length (n) != 0)
+    {
+      octave_idx_type l, u;
+
+      bool col_vec = ndims () == 2 && columns () == 1 && rows () != 1;
+      if (i.is_cont_range (n, l, u))
+        {
+          // Special case deleting a contiguous range.
+          octave_idx_type m = n + l - u;
+          Array<T> tmp (dim_vector (col_vec ? m : 1, !col_vec ? m : 1));
+          const T *src = data ();
+          T *dest = tmp.fortran_vec ();
+          dest = std::copy (src, src + l, dest);
+          dest = std::copy (src + u, src + n, dest);
+          *this = tmp;
+        }
+      else
+        {
+          // Use index.
+          *this = index (i.complement (n));
+        }
+    }
+}
+
+template <class T>
+void 
+Array<T>::delete_elements (int dim, const idx_vector& i)
+{
+  if (dim > ndims ())
+    {
+      (*current_liboctave_error_handler)
+        ("invalid dimension in delete_elements");
+      return;
     }
 
-  if (same_size)
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = rep;
-  const T *old_data = data ();
-
-  octave_idx_type len = get_size (dv);
-
-  rep = new typename Array<T>::ArrayRep (len);
+  octave_idx_type n = dimensions (dim);
+  if (i.is_colon ())
+    { 
+      *this = Array<T> ();
+    }
+  else if (i.extent (n) != n)
+    {
+      gripe_index_out_of_range ();
+    }
+  else if (i.length (n) != 0)
+    {
+      octave_idx_type l, u;
 
-  dim_vector dv_old = dimensions;
-  octave_idx_type dv_old_orig_len = dv_old.length ();
-  dimensions = dv;
+      if (i.is_cont_range (n, l, u))
+        {
+          // Special case deleting a contiguous range.
+          octave_idx_type nd = n + l - u, dl = 1, du = 1;
+          dim_vector rdv = dimensions;
+          rdv(dim) = nd;
+          for (int k = 0; k < dim; k++) dl *= dimensions(k);
+          for (int k = dim + 1; k < ndims (); k++) du *= dimensions(k);
 
-  if (len > 0 && dv_old_orig_len > 0)
-    {
-      Array<octave_idx_type> ra_idx (dimensions.length (), 0);
-      
-      if (n > dv_old_orig_len)
-	{
-	  dv_old.resize (n);
+          // Special case deleting a contiguous range.
+          Array<T> tmp = Array<T> (rdv);
+          const T *src = data ();
+          T *dest = tmp.fortran_vec ();
+          l *= dl; u *= dl; n *= dl;
+          for (octave_idx_type k = 0; k < du; k++)
+            {
+              dest = std::copy (src, src + l, dest);
+              dest = std::copy (src + u, src + n, dest);
+              src += n;
+            }
 
-	  for (octave_idx_type i = dv_old_orig_len; i < n; i++)
-	    dv_old.elem (i) = 1;
-	}
+          *this = tmp;
+        }
+      else
+        {
+          // Use index.
+          Array<idx_vector> ia (ndims (), idx_vector::colon);
+          ia (dim) = i.complement (n);
+          *this = index (ia);
+        }
+    }
+}
 
-      for (octave_idx_type i = 0; i < len; i++)
-	{
-	  if (index_in_bounds (ra_idx, dv_old))
-	    rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)];
-	  else
-	    rep->elem (i) = val;
-	  
-	  increment_index (ra_idx, dimensions);
-	}
+template <class T>
+void 
+Array<T>::delete_elements (const Array<idx_vector>& ia)
+{
+  if (ia.length () == 1)
+    delete_elements (ia(0));
+  else
+    {
+      int len = ia.length (), k, dim = -1;
+      for (k = 0; k < len; k++)
+        {
+          if (! ia(k).is_colon ())
+            {
+              if (dim < 0)
+                dim = k;
+              else
+                break;
+            }
+        }
+      if (dim < 0)
+        {
+          dim_vector dv = dimensions;
+          dv(0) = 0;
+          *this = Array<T> (dv);
+        }
+      else if (k == len)
+        {
+          delete_elements (dim, ia(dim));
+        }
+      else
+        {
+          (*current_liboctave_error_handler)
+            ("A null assignment can only have one non-colon index.");
+        }
     }
-  else
-    for (octave_idx_type i = 0; i < len; i++)
-      rep->elem (i) = val;
 
-  if (--old_rep->count <= 0)
-    delete old_rep;
 }
 
+// FIXME: Remove these methods or implement them using assign.
+
 template <class T>
 Array<T>&
 Array<T>::insert (const Array<T>& a, octave_idx_type r, octave_idx_type c)
@@ -1194,6 +1628,7 @@
   return *this;
 }
 
+
 template <class T>
 Array<T>
 Array<T>::transpose (void) const
@@ -1406,1178 +1841,6 @@
 }
 
 template <class T>
-void
-Array<T>::clear_index (void) const
-{
-  delete [] idx;
-  idx = 0;
-  idx_count = 0;
-}
-
-template <class T>
-void
-Array<T>::set_index (const idx_vector& idx_arg) const
-{
-  int nd = ndims ();
-
-  if (! idx && nd > 0)
-    idx = new idx_vector [nd];
-
-  if (idx_count < nd)
-    {
-      idx[idx_count++] = idx_arg;
-    }
-  else
-    {
-      idx_vector *new_idx = new idx_vector [idx_count+1];
-
-      for (int i = 0; i < idx_count; i++)
-	new_idx[i] = idx[i];
-
-      new_idx[idx_count++] = idx_arg;
-
-      delete [] idx;
-
-      idx = new_idx;
-    }
-}
-
-template <class T>
-void
-Array<T>::maybe_delete_elements (idx_vector& idx_arg)
-{
-  switch (ndims ())
-    {
-    case 1:
-      maybe_delete_elements_1 (idx_arg);
-      break;
-
-    case 2:
-      maybe_delete_elements_2 (idx_arg);
-      break;
-
-    default:
-      (*current_liboctave_error_handler)
-	("Array<T>::maybe_delete_elements: invalid operation");
-      break;
-    }
-}
-
-template <class T>
-void
-Array<T>::maybe_delete_elements_1 (idx_vector& idx_arg)
-{
-  octave_idx_type len = length ();
-
-  if (len == 0)
-    return;
-
-  if (idx_arg.is_colon_equiv (len, 1))
-    resize_no_fill (0);
-  else
-    {
-      int num_to_delete = idx_arg.length (len);
-
-      if (num_to_delete != 0)
-	{
-	  octave_idx_type new_len = len;
-
-	  octave_idx_type iidx = 0;
-
-	  for (octave_idx_type i = 0; i < len; i++)
-	    if (i == idx_arg.elem (iidx))
-	      {
-		iidx++;
-		new_len--;
-
-		if (iidx == num_to_delete)
-		  break;
-	      }
-
-	  if (new_len > 0)
-	    {
-	      T *new_data = new T [new_len];
-
-	      octave_idx_type ii = 0;
-	      iidx = 0;
-	      for (octave_idx_type i = 0; i < len; i++)
-		{
-		  if (iidx < num_to_delete && i == idx_arg.elem (iidx))
-		    iidx++;
-		  else
-		    {
-		      new_data[ii] = xelem (i);
-		      ii++;
-		    }
-		}
-
-	      if (--rep->count <= 0)
-		delete rep;
-
-	      rep = new typename Array<T>::ArrayRep (new_data, new_len);
-
-	      dimensions.resize (1);
-	      dimensions(0) = new_len;
-	    }
-	  else
-	    (*current_liboctave_error_handler)
-	      ("A(idx) = []: index out of range");
-	}
-    }
-}
-
-template <class T>
-void
-Array<T>::maybe_delete_elements_2 (idx_vector& idx_arg)
-{
-  assert (ndims () == 2);
-
-  octave_idx_type nr = dim1 ();
-  octave_idx_type nc = dim2 ();
-
-  if (idx_arg.is_colon ())
-    {
-      // A(:) = [] always gives 0-by-0 matrix, even if A was empty.
-      resize_no_fill (0, 0);
-      return;
-    }
-
-  octave_idx_type n;
-  if (nr == 1)
-    n = nc;
-  else if (nc == 1)
-    n = nr;
-  else if (! idx_arg.orig_empty ())
-    {
-      // Reshape to row vector for Matlab compatibility.
-
-      n = nr * nc;
-      nr = 1;
-      nc = n;
-    }
-  else
-    return;
-
-  idx_arg.sort (true);
-
-  if (idx_arg.is_colon_equiv (n, 1))
-    {
-      if (nr == 1)
-        resize_no_fill (1, 0);
-      else if (nc == 1)
-        resize_no_fill (0, 1);
-      return;
-    }
-
-  octave_idx_type num_to_delete = idx_arg.length (n);
-
-  if (num_to_delete != 0)
-    {
-      octave_idx_type new_n = n;
-
-      octave_idx_type iidx = 0;
-
-      for (octave_idx_type i = 0; i < n; i++)
-	if (i == idx_arg.elem (iidx))
-	  {
-	    iidx++;
-	    new_n--;
-
-	    if (iidx == num_to_delete)
-	      break;
-	  }
-
-      if (new_n > 0)
-	{
-	  T *new_data = new T [new_n];
-
-	  octave_idx_type ii = 0;
-	  iidx = 0;
-	  for (octave_idx_type i = 0; i < n; i++)
-	    {
-	      if (iidx < num_to_delete && i == idx_arg.elem (iidx))
-		iidx++;
-	      else
-		{
-		  new_data[ii] = xelem (i);
-
-		  ii++;
-		}
-	    }
-
-	  if (--(Array<T>::rep)->count <= 0)
-	    delete Array<T>::rep;
-
-	  Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_n);
-
-	  dimensions.resize (2);
-
-	  if (nr == 1)
-	    {
-	      dimensions(0) = 1;
-	      dimensions(1) = new_n;
-	    }
-	  else
-	    {
-	      dimensions(0) = new_n;
-	      dimensions(1) = 1;
-	    }
-	}
-      else
-	(*current_liboctave_error_handler)
-	  ("A(idx) = []: index out of range");
-    }
-}
-
-template <class T>
-void
-Array<T>::maybe_delete_elements (idx_vector& idx_i, idx_vector& idx_j)
-{
-  assert (ndims () == 2);
-
-  octave_idx_type nr = dim1 ();
-  octave_idx_type nc = dim2 ();
-
-  if (idx_i.is_colon () && idx_j.is_colon ())
-    {
-      // A special case: A(:,:). Matlab gives 0-by-nc here, but perhaps we
-      // should not?
-      resize_no_fill (0, nc);
-    }
-  else if (idx_i.is_colon ())
-    {
-      idx_j.sort (true); // sort in advance to speed-up the following check
-
-      if (idx_j.is_colon_equiv (nc, 1))
-	resize_no_fill (nr, 0);
-      else
-	{
-	  octave_idx_type num_to_delete = idx_j.length (nc);
-
-	  if (num_to_delete != 0)
-            {
-              octave_idx_type new_nc = nc;
-
-              octave_idx_type iidx = 0;
-
-              for (octave_idx_type j = 0; j < nc; j++)
-                if (j == idx_j.elem (iidx))
-                  {
-                    iidx++;
-                    new_nc--;
-
-                    if (iidx == num_to_delete)
-                      break;
-                  }
-
-              if (new_nc > 0)
-                {
-                  T *new_data = new T [nr * new_nc];
-
-                  octave_idx_type jj = 0;
-                  iidx = 0;
-                  for (octave_idx_type j = 0; j < nc; j++)
-                    {
-                      if (iidx < num_to_delete && j == idx_j.elem (iidx))
-                        iidx++;
-                      else
-                        {
-                          for (octave_idx_type i = 0; i < nr; i++)
-                            new_data[nr*jj+i] = xelem (i, j);
-                          jj++;
-                        }
-                    }
-
-                  if (--(Array<T>::rep)->count <= 0)
-                    delete Array<T>::rep;
-
-                  Array<T>::rep = new typename Array<T>::ArrayRep (new_data, nr * new_nc);
-
-                  dimensions.resize (2);
-                  dimensions(1) = new_nc;
-                }
-              else
-                (*current_liboctave_error_handler)
-                  ("A(idx) = []: index out of range");
-            }
-	}
-    }
-  else if (idx_j.is_colon ())
-    {
-      idx_i.sort (true); // sort in advance to speed-up the following check
-
-      if (idx_i.is_colon_equiv (nr, 1))
-	resize_no_fill (0, nc);
-      else
-	{
-	  octave_idx_type num_to_delete = idx_i.length (nr);
-
-	  if (num_to_delete != 0)
-            {
-              octave_idx_type new_nr = nr;
-
-              octave_idx_type iidx = 0;
-
-              for (octave_idx_type i = 0; i < nr; i++)
-                if (i == idx_i.elem (iidx))
-                  {
-                    iidx++;
-                    new_nr--;
-
-                    if (iidx == num_to_delete)
-                      break;
-                  }
-
-              if (new_nr > 0)
-                {
-                  T *new_data = new T [new_nr * nc];
-
-                  octave_idx_type ii = 0;
-                  iidx = 0;
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    {
-                      if (iidx < num_to_delete && i == idx_i.elem (iidx))
-                        iidx++;
-                      else
-                        {
-                          for (octave_idx_type j = 0; j < nc; j++)
-                            new_data[new_nr*j+ii] = xelem (i, j);
-                          ii++;
-                        }
-                    }
-
-                  if (--(Array<T>::rep)->count <= 0)
-                    delete Array<T>::rep;
-
-                  Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_nr * nc);
-
-                  dimensions.resize (2);
-                  dimensions(0) = new_nr;
-                }
-              else
-                (*current_liboctave_error_handler)
-                  ("A(idx) = []: index out of range");
-            }
-	}
-    }
-  else if (! (idx_i.orig_empty () || idx_j.orig_empty ()))
-    {
-      (*current_liboctave_error_handler)
-        ("a null assignment can have only one non-colon index");
-    }
-}
-
-template <class T>
-void
-Array<T>::maybe_delete_elements (idx_vector&, idx_vector&, idx_vector&)
-{
-  assert (0);
-}
-
-template <class T>
-void
-Array<T>::maybe_delete_elements (Array<idx_vector>& ra_idx)
-{
-  octave_idx_type n_idx = ra_idx.length ();
-
-  // Special case matrices
-  if (ndims () == 2)
-    {
-      if (n_idx == 1)
-        {
-          maybe_delete_elements (ra_idx (0));
-          return;
-        }
-      else if (n_idx == 2)
-        {
-          maybe_delete_elements (ra_idx (0), ra_idx (1));
-          return;
-        }
-    }
-
-  dim_vector lhs_dims = dims ();
-
-  int n_lhs_dims = lhs_dims.length ();
-
-  if (lhs_dims.all_zero ())
-    return;
-
-  if (n_idx == 1 && ra_idx(0).is_colon ())
-    {
-      resize (dim_vector (0, 0));
-      return;
-    }
-
-  if (n_idx > n_lhs_dims)
-    {
-      for (int i = n_idx; i < n_lhs_dims; i++)
-	{
-	  // Ensure that extra indices are either colon or 1.
-
-	  if (! ra_idx(i).is_colon_equiv (1, 1))
-	    {
-	      (*current_liboctave_error_handler)
-		("index exceeds array dimensions");
-	      return;
-	    }
-	}
-
-      ra_idx.resize (n_lhs_dims);
-
-      n_idx = n_lhs_dims;
-    }
-
-  Array<int> idx_is_colon (n_idx, 0);
-
-  Array<int> idx_is_colon_equiv (n_idx, 0);
-
-  // Initialization of colon arrays.
-
-  for (octave_idx_type i = 0; i < n_idx; i++)
-    {
-      if (ra_idx(i).orig_empty ()) return;
-      idx_is_colon_equiv(i) = ra_idx(i).is_colon_equiv (lhs_dims(i), 1);
-
-      idx_is_colon(i) = ra_idx(i).is_colon ();
-    }
-
-  bool idx_ok = true;
-
-  // Check for index out of bounds.
-
-  for (octave_idx_type i = 0 ; i < n_idx - 1; i++)
-    {
-      if (! (idx_is_colon(i) || idx_is_colon_equiv(i)))
-	{
-	  ra_idx(i).sort (true);
-
-	  if (ra_idx(i).max () > lhs_dims(i))
-	    {
-	      (*current_liboctave_error_handler)
-		("index exceeds array dimensions");
-
-	      idx_ok = false;
-	      break;
-	    }
-	  else if (ra_idx(i).min () < 0) // I believe this is checked elsewhere
-	    {
-	      (*current_liboctave_error_handler)
-		("index must be one or larger");
-
-	      idx_ok = false;
-	      break;
-	    }
-	}
-    }
-
-  if (n_idx <= n_lhs_dims)
-    {
-      octave_idx_type last_idx = ra_idx(n_idx-1).max ();
-
-      octave_idx_type sum_el = lhs_dims(n_idx-1);
-
-      for (octave_idx_type i = n_idx; i < n_lhs_dims; i++)
-	  sum_el *= lhs_dims(i);
-
-      if (last_idx > sum_el - 1)
-	{
-	  (*current_liboctave_error_handler)
-	    ("index exceeds array dimensions");
-
-	  idx_ok = false;
-	}
-    }
-
-  if (idx_ok)
-    {
-      if (n_idx > 1
-	  && (all_ones (idx_is_colon)))
-	{
-	  // A(:,:,:) -- we are deleting elements in all dimensions, so
-	  // the result is [](0x0x0).
-
-	  dim_vector newdim = dims ();
-          newdim(0) = 0;
-
-	  resize (newdim);
-	}
-
-      else if (n_idx > 1
-	       && num_ones (idx_is_colon) == n_idx - 1
-	       && num_ones (idx_is_colon_equiv) == n_idx)
-	{
-	  // A(:,:,j) -- we are deleting elements in one dimension by
-	  // enumerating them.
-	  //
-	  // If we enumerate all of the elements, we should have zero
-	  // elements in that dimension with the same number of elements
-	  // in the other dimensions that we started with.
-
-	  dim_vector temp_dims;
-	  temp_dims.resize (n_idx);
-
-	  for (octave_idx_type i = 0; i < n_idx; i++)
-	    {
-	      if (idx_is_colon (i))
-		temp_dims(i) =  lhs_dims(i);
-	      else
-		temp_dims(i) = 0;
-	    }
-
-	  resize (temp_dims);
-	}
-      else if (n_idx > 1 && num_ones (idx_is_colon) == n_idx - 1)
-	{
-	  // We have colons in all indices except for one.
-	  // This index tells us which slice to delete
-
-	  if (n_idx < n_lhs_dims)
-	    {
-	      // Collapse dimensions beyond last index.
-
-	      if (! (ra_idx(n_idx-1).is_colon ()))
-		(*current_liboctave_warning_with_id_handler)
-		  ("Octave:fortran-indexing",
-		   "fewer indices than dimensions for N-d array");
-
-	      for (octave_idx_type i = n_idx; i < n_lhs_dims; i++)
-		lhs_dims(n_idx-1) *= lhs_dims(i);
-
-	      lhs_dims.resize (n_idx);
-
-	      // Reshape *this.
-	      dimensions = lhs_dims;
-	    }
-
-	  int non_col = 0;
-
-	  // Find the non-colon column.
-
-	  for (octave_idx_type i = 0; i < n_idx; i++)
-	    {
-	      if (! idx_is_colon(i))
-		non_col = i;
-	    }
-
-	  // The length of the non-colon dimension.
-
-	  octave_idx_type non_col_dim = lhs_dims (non_col);
-
-	  octave_idx_type num_to_delete = ra_idx(non_col).length (lhs_dims (non_col));
-
-	  if (num_to_delete > 0)
-	    {
-	      int temp = lhs_dims.num_ones ();
-
-	      if (non_col_dim == 1)
-		temp--;
-
-	      if (temp == n_idx - 1 && num_to_delete == non_col_dim)
-		{
-		  // We have A with (1x1x4), where A(1,:,1:4)
-		  // Delete all (0x0x0)
-
-		  dim_vector zero_dims (n_idx, 0);
-
-		  resize (zero_dims);
-		}
-	      else
-		{
-		  // New length of non-colon dimension
-		  // (calculated in the next for loop)
-
-		  octave_idx_type new_dim = non_col_dim;
-
-		  octave_idx_type iidx = 0;
-
-		  for (octave_idx_type j = 0; j < non_col_dim; j++)
-		    if (j == ra_idx(non_col).elem (iidx))
-		      {
-			iidx++;
-
-			new_dim--;
-
-			if (iidx == num_to_delete)
-			  break;
-		      }
-
-		  // Creating the new nd array after deletions.
-
-		  if (new_dim > 0)
-		    {
-		      // Calculate number of elements in new array.
-
-		      octave_idx_type num_new_elem=1;
-
-		      for (int i = 0; i < n_idx; i++)
-			{
-			  if (i == non_col)
-			    num_new_elem *= new_dim;
-
-			  else
-			    num_new_elem *= lhs_dims(i);
-			}
-
-		      T *new_data = new T [num_new_elem];
-
-		      Array<octave_idx_type> result_idx (n_lhs_dims, 0);
-
-		      dim_vector new_lhs_dim = lhs_dims;
-
-		      new_lhs_dim(non_col) = new_dim;
-
-		      octave_idx_type num_elem = 1;
-
-		      octave_idx_type numidx = 0;
-
-		      octave_idx_type n = length ();
-
-		      for (int i = 0; i < n_lhs_dims; i++)
-			if (i != non_col)
-			  num_elem *= lhs_dims(i);
-
-		      num_elem *= ra_idx(non_col).capacity ();
-
-		      for (octave_idx_type i = 0; i < n; i++)
-			{
-			  if (numidx < num_elem
-			      && is_in (result_idx(non_col), ra_idx(non_col)))
-			    numidx++;
-
-			  else
-			    {
-			      Array<octave_idx_type> temp_result_idx = result_idx;
-
-			      octave_idx_type num_lgt = how_many_lgt (result_idx(non_col),
-							  ra_idx(non_col));
-
-			      temp_result_idx(non_col) -= num_lgt;
-
-			      octave_idx_type kidx
-				= ::compute_index (temp_result_idx, new_lhs_dim);
-
-			      new_data[kidx] = xelem (result_idx);
-			    }
-
-			  increment_index (result_idx, lhs_dims);
-			}
-
-		      if (--rep->count <= 0)
-			delete rep;
-
-		      rep = new typename Array<T>::ArrayRep (new_data,
-							     num_new_elem);
-
-		      dimensions = new_lhs_dim;
-		    }
-		}
-	    }
-	}
-      else if (n_idx == 1)
-	{
-	  // This handle cases where we only have one index (not
-	  // colon).  The index denotes which elements we should
-	  // delete in the array which can be of any dimension. We
-	  // return a column vector, except for the case where we are
-	  // operating on a row vector. The elements are numerated
-	  // column by column.
-	  //
-	  // A(3,3,3)=2;
-	  // A(3:5) = []; A(6)=[]
-
-	  octave_idx_type lhs_numel = numel ();
-
-	  idx_vector idx_vec = ra_idx(0);
-
-	  idx_vec.freeze (lhs_numel, 0, true);
-      
-	  idx_vec.sort (true);
-
-	  octave_idx_type num_to_delete = idx_vec.length (lhs_numel);
-
-	  if (num_to_delete > 0)
-	    {
-	      octave_idx_type new_numel = lhs_numel - num_to_delete;
-
-	      T *new_data = new T[new_numel];
-
-	      Array<octave_idx_type> lhs_ra_idx (ndims (), 0);
-
-	      octave_idx_type ii = 0;
-	      octave_idx_type iidx = 0;
-
-	      for (octave_idx_type i = 0; i < lhs_numel; i++)
-		{
-		  if (iidx < num_to_delete && i == idx_vec.elem (iidx))
-		    {
-		      iidx++;
-		    }
-		  else
-		    {
-		      new_data[ii++] = xelem (lhs_ra_idx);
-		    }
-
-		  increment_index (lhs_ra_idx, lhs_dims);
-		}
-
-	      if (--(Array<T>::rep)->count <= 0)
-		delete Array<T>::rep;
-
-	      Array<T>::rep = new typename Array<T>::ArrayRep (new_data, new_numel);
-
-	      dimensions.resize (2);
-
-	      if (lhs_dims.length () == 2 && lhs_dims(1) == 1)
-		{
-		  dimensions(0) = new_numel;
-		  dimensions(1) = 1;
-		}
-	      else
-		{
-		  dimensions(0) = 1;
-		  dimensions(1) = new_numel;
-		}
-	    }
-	}
-      else if (num_ones (idx_is_colon) < n_idx)
-	{
-	  (*current_liboctave_error_handler)
-	    ("a null assignment can have only one non-colon index");
-	}
-    }
-}
-
-template <class T>
-Array<T>
-Array<T>::value (void) const
-{
-  Array<T> retval;
-
-  int n_idx = index_count ();
-
-  if (n_idx == 2)
-    {
-      idx_vector *tmp = get_idx ();
-
-      idx_vector idx_i = tmp[0];
-      idx_vector idx_j = tmp[1];
-
-      retval = index (idx_i, idx_j);
-    }
-  else if (n_idx == 1)
-    {
-      retval = index (idx[0]);
-    }
-  else
-    {
-      clear_index ();
-
-      (*current_liboctave_error_handler)
-	("Array<T>::value: invalid number of indices specified");
-    }
-
-  clear_index ();
-
-  return retval;
-}
-
-template <class T>
-Array<T>
-Array<T>::index (idx_vector& idx_arg, int resize_ok, const T& rfv) const
-{
-  Array<T> retval;
-
-  dim_vector dv = idx_arg.orig_dimensions ();
-
-  if (dv.length () > 2 || ndims () > 2)
-    retval = indexN (idx_arg, resize_ok, rfv);
-  else
-    {
-      switch (ndims ())
-	{
-	case 1:
-	  retval = index1 (idx_arg, resize_ok, rfv);
-	  break;
-
-	case 2:
-	  retval = index2 (idx_arg, resize_ok, rfv);
-	  break;
-
-	default:
-	  (*current_liboctave_error_handler)
-	    ("invalid array (internal error)");
-	  break;
-	}
-    }
-
-  return retval;
-}
-
-template <class T>
-Array<T>
-Array<T>::index1 (idx_vector& idx_arg, int resize_ok, const T& rfv) const
-{
-  Array<T> retval;
-
-  octave_idx_type len = length ();
-
-  octave_idx_type n = idx_arg.freeze (len, "vector", resize_ok);
-
-  if (idx_arg)
-    {
-      if (idx_arg.is_colon_equiv (len))
-	{
-	  retval = *this;
-	}
-      else if (n == 0)
-	{
-	  retval.resize_no_fill (0);
-	}
-      else
-	{
-	  retval.resize_no_fill (n);
-
-	  for (octave_idx_type i = 0; i < n; i++)
-	    {
-	      octave_idx_type ii = idx_arg.elem (i);
-	      if (ii >= len)
-		retval.elem (i) = rfv;
-	      else
-		retval.elem (i) = elem (ii);
-	    }
-	}
-    }
-
-  // idx_vector::freeze() printed an error message for us.
-
-  return retval;
-}
-
-template <class T>
-Array<T>
-Array<T>::index2 (idx_vector& idx_arg, int resize_ok, const T& rfv) const
-{
-  Array<T> retval;
-
-  assert (ndims () == 2);
-
-  octave_idx_type nr = dim1 ();
-  octave_idx_type nc = dim2 ();
-
-  octave_idx_type orig_len = nr * nc;
-
-  dim_vector idx_orig_dims = idx_arg.orig_dimensions ();
-
-  octave_idx_type idx_orig_rows = idx_arg.orig_rows ();
-  octave_idx_type idx_orig_columns = idx_arg.orig_columns ();
-
-  if (idx_arg.is_colon ())
-    {
-      // Fast magic colon processing.
-
-      octave_idx_type result_nr = nr * nc;
-      octave_idx_type result_nc = 1;
-
-      retval = Array<T> (*this, dim_vector (result_nr, result_nc));
-    }
-  else if (nr == 1 && nc == 1)
-    {
-      Array<T> tmp = Array<T>::index1 (idx_arg, resize_ok);
-
-      octave_idx_type len = tmp.length ();
-
-      if (len >= idx_orig_dims.numel ())
-	retval = Array<T> (tmp, idx_orig_dims);
-    }
-  else if (nr == 1 || nc == 1)
-    {
-      // If indexing a vector with a matrix, return value has same
-      // shape as the index.  Otherwise, it has same orientation as
-      // indexed object.
-
-      Array<T> tmp = Array<T>::index1 (idx_arg, resize_ok);
-
-      octave_idx_type len = tmp.length ();
-
-      if (idx_orig_rows == 1 || idx_orig_columns == 1)
-	{
-	  if (nr == 1)
-	    retval = Array<T> (tmp, dim_vector (1, len));
-	  else
-	    retval = Array<T> (tmp, dim_vector (len, 1));
-	}
-      else if (len >= idx_orig_dims.numel ())
-	retval = Array<T> (tmp, idx_orig_dims);
-    }
-  else
-    {
-      (*current_liboctave_warning_with_id_handler)
-	("Octave:fortran-indexing", "single index used for matrix");
-
-      // This code is only for indexing matrices.  The vector
-      // cases are handled above.
-
-      idx_arg.freeze (nr * nc, "matrix", resize_ok);
-
-      if (idx_arg)
-	{
-	  octave_idx_type result_nr = idx_orig_rows;
-	  octave_idx_type result_nc = idx_orig_columns;
-
-	  retval.resize_no_fill (result_nr, result_nc);
-
-	  octave_idx_type k = 0;
-	  for (octave_idx_type j = 0; j < result_nc; j++)
-	    {
-	      for (octave_idx_type i = 0; i < result_nr; i++)
-		{
-		  octave_idx_type ii = idx_arg.elem (k++);
-		  if (ii >= orig_len)
-		    retval.elem (i, j) = rfv;
-		  else
-		    {
-		      octave_idx_type fr = ii % nr;
-		      octave_idx_type fc = (ii - fr) / nr;
-		      retval.elem (i, j) = elem (fr, fc);
-		    }
-		}
-	    }
-	}
-      // idx_vector::freeze() printed an error message for us.
-    }
-
-  return retval;
-}
-
-template <class T>
-Array<T>
-Array<T>::indexN (idx_vector& ra_idx, int resize_ok, const T& rfv) const
-{
-  Array<T> retval;
-
-  dim_vector dv = dims ();
-
-  int n_dims = dv.length ();
-
-  octave_idx_type orig_len = dv.numel ();
-
-  dim_vector idx_orig_dims = ra_idx.orig_dimensions ();
-
-  if (ra_idx.is_colon ())
-    {
-      // Fast magic colon processing.
-
-      retval = Array<T> (*this, dim_vector (orig_len, 1));
-    }
-  else
-    {
-      bool vec_equiv = vector_equivalent (dv);
-
-      if (! (vec_equiv || ra_idx.is_colon ()))
-	(*current_liboctave_warning_with_id_handler)
-	  ("Octave:fortran-indexing", "single index used for N-d array");
-
-      octave_idx_type frozen_len
-	= ra_idx.freeze (orig_len, "nd-array", resize_ok);
-
-      if (ra_idx)
-	{
-	  dim_vector result_dims;
-
-	  if (vec_equiv && ! orig_len == 1)
-	    {
-	      result_dims = dv;
-
-	      for (int i = 0; i < n_dims; i++)
-		{
-		  if (result_dims(i) != 1)
-		    {
-		      // All but this dim should be one.
-		      result_dims(i) = frozen_len;
-		      break;
-		    }
-		}
-	    }
-	  else
-	    result_dims = idx_orig_dims;
-
-	  result_dims.chop_trailing_singletons ();
-
-	  retval.resize (result_dims);
-
-	  octave_idx_type n = result_dims.numel ();
-
-	  octave_idx_type k = 0;
-
-	  for (octave_idx_type i = 0; i < n; i++)
-	    {
-	      octave_idx_type ii = ra_idx.elem (k++);
-
-	      if (ii >= orig_len)
-	        retval.elem (i) = rfv;
-	      else
-		retval.elem (i) = elem (ii);
-	    }
-	}
-    }
-
-  return retval;
-}
-
-template <class T>
-Array<T>
-Array<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok,
-		 const T& rfv) const
-{
-  Array<T> retval;
-
-  if (ndims () != 2)
-    {
-      Array<idx_vector> ra_idx (2);
-      ra_idx(0) = idx_i;
-      ra_idx(1) = idx_j;
-      return index (ra_idx, resize_ok, rfv);
-    }
-
-  octave_idx_type nr = dim1 ();
-  octave_idx_type nc = dim2 ();
-
-  octave_idx_type n = idx_i.freeze (nr, "row", resize_ok);
-  octave_idx_type m = idx_j.freeze (nc, "column", resize_ok);
-
-  if (idx_i && idx_j)
-    {
-      if (idx_i.orig_empty () || idx_j.orig_empty () || n == 0 || m == 0)
-	{
-	  retval.resize_no_fill (n, m);
-	}
-      else if (idx_i.is_colon_equiv (nr) && idx_j.is_colon_equiv (nc))
-	{
-	  retval = *this;
-	}
-      else
-	{
-	  retval.resize_no_fill (n, m);
-
-	  for (octave_idx_type j = 0; j < m; j++)
-	    {
-	      octave_idx_type jj = idx_j.elem (j);
-	      for (octave_idx_type i = 0; i < n; i++)
-		{
-		  octave_idx_type ii = idx_i.elem (i);
-		  if (ii >= nr || jj >= nc)
-		    retval.elem (i, j) = rfv;
-		  else
-		    retval.elem (i, j) = elem (ii, jj);
-		}
-	    }
-	}
-    }
-
-  // idx_vector::freeze() printed an error message for us.
-
-  return retval;
-}
-
-template <class T>
-Array<T>
-Array<T>::index (Array<idx_vector>& ra_idx, int resize_ok, const T& rfv) const
-{
-  // This function handles all calls with more than one idx.
-  // For (3x3x3), the call can be A(2,5), A(2,:,:), A(3,2,3) etc.
-
-  Array<T> retval;
-
-  int n_dims = dimensions.length ();
-
-  // Remove trailing singletons in ra_idx, but leave at least ndims
-  // elements.
-
-  octave_idx_type ra_idx_len = ra_idx.length ();
-
-  bool trim_trailing_singletons = true;
-  for (octave_idx_type j = ra_idx_len; j > n_dims; j--)
-    {
-      idx_vector iidx = ra_idx (ra_idx_len-1);
-      if (iidx.capacity () == 1 && trim_trailing_singletons)
-	ra_idx_len--;
-      else
-	trim_trailing_singletons = false;
-
-      if (! resize_ok)
-	{
-	  for (octave_idx_type i = 0; i < iidx.capacity (); i++)
-	    if (iidx (i) != 0)
-	      {
-		(*current_liboctave_error_handler)
-		  ("index exceeds N-d array dimensions");
-
-		return retval;
-	      }
-	}
-    }
-
-  ra_idx.resize (ra_idx_len);
-
-  dim_vector new_dims = dims ();
-  dim_vector frozen_lengths;
-
-  if (!ra_idx (ra_idx_len - 1).orig_empty () && ra_idx_len < n_dims)
-    frozen_lengths = short_freeze (ra_idx, dimensions, resize_ok);
-  else
-    {
-      new_dims.resize (ra_idx_len, 1);
-      frozen_lengths = freeze (ra_idx, new_dims, resize_ok);
-    }
-
-  if (all_ok (ra_idx))
-    {
-      if (any_orig_empty (ra_idx) || frozen_lengths.any_zero ())
-	{
-	  frozen_lengths.chop_trailing_singletons ();
-
-	  retval.resize (frozen_lengths);
-	}
-      else if (frozen_lengths.length () == n_dims
-	       && all_colon_equiv (ra_idx, dimensions))
-	{
-	  retval = *this;
-	}
-      else
-	{
-	  dim_vector frozen_lengths_for_resize = frozen_lengths;
-
-	  frozen_lengths_for_resize.chop_trailing_singletons ();
-
-	  retval.resize (frozen_lengths_for_resize);
-
-	  octave_idx_type n = retval.length ();
-
-	  Array<octave_idx_type> result_idx (ra_idx.length (), 0);
-
-	  Array<octave_idx_type> elt_idx;
-
-	  for (octave_idx_type i = 0; i < n; i++)
-	    {
-	      elt_idx = get_elt_idx (ra_idx, result_idx);
-
-	      octave_idx_type numelem_elt = get_scalar_idx (elt_idx, new_dims);
-
-	      if (numelem_elt >= length () || numelem_elt < 0)
-		retval.elem (i) = rfv;
-	      else
-		retval.elem (i) = elem (numelem_elt);
-
-	      increment_index (result_idx, frozen_lengths);
-
-	    }
-	}
-    }
-
-  return retval;
-}
-
-template <class T>
 bool 
 ascending_compare (T a, T b)
 {
@@ -2851,7 +2114,7 @@
 	  if (nnr == 1)
 	    {
 	      octave_idx_type n = nnc + std::abs (k);
-	      d = Array<T> (dim_vector (n, n), resize_fill_value (T ()));
+	      d = Array<T> (dim_vector (n, n), resize_fill_value ());
 
 	      for (octave_idx_type i = 0; i < nnc; i++)
 		d.xelem (i+roff, i+coff) = elem (0, i);
@@ -2859,7 +2122,7 @@
 	  else
 	    {
 	      octave_idx_type n = nnr + std::abs (k);
-	      d = Array<T> (dim_vector (n, n), resize_fill_value (T ()));
+	      d = Array<T> (dim_vector (n, n), resize_fill_value ());
 
 	      for (octave_idx_type i = 0; i < nnr; i++)
 		d.xelem (i+roff, i+coff) = elem (i, 0);
@@ -2870,1011 +2133,6 @@
   return d;
 }
 
-// FIXME -- this is a mess.
-
-template <class LT, class RT>
-int
-assign (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
-{
-  int n_idx = lhs.index_count ();
-
-  // kluge...
-  if (lhs.ndims () == 0)
-    lhs.resize_no_fill (0, 0);
-
-  return (lhs.ndims () == 2
-	  && (n_idx == 1
-	      || (n_idx < 3
-		  && rhs.ndims () == 2
-		  && rhs.rows () == 0 && rhs.columns () == 0)))
-    ? assign2 (lhs, rhs, rfv) : assignN (lhs, rhs, rfv);
-}
-
-template <class LT, class RT>
-int
-assign1 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
-{
-  int retval = 1;
-
-  idx_vector *tmp = lhs.get_idx ();
-
-  idx_vector lhs_idx = tmp[0];
-
-  octave_idx_type lhs_len = lhs.length ();
-  octave_idx_type rhs_len = rhs.length ();
-
-  octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true);
-
-  if (n != 0)
-    {
-      dim_vector lhs_dims = lhs.dims ();
-
-      if (lhs_len != 0
-	  || lhs_dims.all_zero ()
-	  || (lhs_dims.length () == 2 && lhs_dims(0) < 2))
-	{
-	  if (rhs_len == n || rhs_len == 1)
-	    {
-	      lhs.make_unique ();
-
-	      octave_idx_type max_idx = lhs_idx.max () + 1;
-	      if (max_idx > lhs_len)
-		lhs.resize_and_fill (max_idx, rfv);
-	    }
-
-	  if (rhs_len == n)
-	    {
-	      lhs.make_unique ();
-
-	      if (lhs_idx.is_colon ())
-		{
-		  for (octave_idx_type i = 0; i < n; i++)
-		    lhs.xelem (i) = rhs.elem (i);
-		}
-	      else
-		{
-		  for (octave_idx_type i = 0; i < n; i++)
-		    {
-		      octave_idx_type ii = lhs_idx.elem (i);
-		      lhs.xelem (ii) = rhs.elem (i);
-		    }
-		}
-	    }
-	  else if (rhs_len == 1)
-	    {
-	      lhs.make_unique ();
-
-	      RT scalar = rhs.elem (0);
-
-	      if (lhs_idx.is_colon ())
-		{
-		  for (octave_idx_type i = 0; i < n; i++)
-		    lhs.xelem (i) = scalar;
-		}
-	      else
-		{
-		  for (octave_idx_type i = 0; i < n; i++)
-		    {
-		      octave_idx_type ii = lhs_idx.elem (i);
-		      lhs.xelem (ii) = scalar;
-		    }
-		}
-	    }
-	  else
-	    {
-	      lhs.clear_index ();
-
-	      (*current_liboctave_error_handler)
-		("A(I) = X: X must be a scalar or a vector with same length as I");
-
-	      retval = 0;
-	    }
-	}
-      else
-	{
-	  lhs.clear_index ();
-
-	  (*current_liboctave_error_handler)
-	    ("A(I) = X: unable to resize A");
-
-	  retval = 0;
-	}
-    }
-  else if (lhs_idx.is_colon ())
-    {
-      dim_vector lhs_dims = lhs.dims ();
-
-      if (lhs_dims.all_zero ())
-	{
-	  lhs.make_unique ();
-
-	  lhs.resize_no_fill (rhs_len);
-
-	  for (octave_idx_type i = 0; i < rhs_len; i++)
-	    lhs.xelem (i) = rhs.elem (i);
-	}
-      else if (rhs_len != lhs_len)
-	{
-	  lhs.clear_index ();
-
-	  (*current_liboctave_error_handler)
-	    ("A(:) = X: A must be the same size as X");
-	}
-    }
-  else if (! (rhs_len == 1 || rhs_len == 0))
-    {
-      lhs.clear_index ();
-
-      (*current_liboctave_error_handler)
-	("A([]) = X: X must also be an empty matrix or a scalar");
-
-      retval = 0;
-    }
-
-  lhs.clear_index ();
-
-  return retval;
-}
-
-#define MAYBE_RESIZE_LHS \
-  do \
-    { \
-      octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : idx_i.max () + 1; \
-      octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : idx_j.max () + 1; \
- \
-      octave_idx_type new_nr = max_row_idx > lhs_nr ? max_row_idx : lhs_nr; \
-      octave_idx_type new_nc = max_col_idx > lhs_nc ? max_col_idx : lhs_nc; \
- \
-      lhs.resize_and_fill (new_nr, new_nc, rfv); \
-    } \
-  while (0)
-
-template <class LT, class RT>
-int
-assign2 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
-{
-  int retval = 1;
-
-  int n_idx = lhs.index_count ();
-
-  octave_idx_type lhs_nr = lhs.rows ();
-  octave_idx_type lhs_nc = lhs.cols ();
-
-  Array<RT> xrhs = rhs;
-
-  octave_idx_type rhs_nr = xrhs.rows ();
-  octave_idx_type rhs_nc = xrhs.cols ();
-
-  if (xrhs.ndims () > 2)
-    {
-      xrhs = xrhs.squeeze ();
-
-      dim_vector dv_tmp = xrhs.dims ();
-
-      switch (dv_tmp.length ())
-	{
-	case 1:
-	  // FIXME -- this case should be unnecessary, because
-	  // squeeze should always return an object with 2 dimensions.
-	  if (rhs_nr == 1)
-	    rhs_nc = dv_tmp.elem (0);
-	  break;
-
-	case 2:
-	  rhs_nr = dv_tmp.elem (0);
-	  rhs_nc = dv_tmp.elem (1);
-	  break;
-
-	default:
-	  lhs.clear_index ();
-	  (*current_liboctave_error_handler)
-	    ("Array<T>::assign2: Dimension mismatch");
-	  return 0;
-	}
-    }
-
-  bool rhs_is_scalar = rhs_nr == 1 && rhs_nc == 1;
-
-  idx_vector *tmp = lhs.get_idx ();
-
-  idx_vector idx_i;
-  idx_vector idx_j;
-
-  if (n_idx > 1)
-    idx_j = tmp[1];
-
-  if (n_idx > 0)
-    idx_i = tmp[0];
-
-  if (n_idx == 2)
-    {
-      octave_idx_type n = idx_i.freeze (lhs_nr, "row", true);
-      octave_idx_type m = idx_j.freeze (lhs_nc, "column", true);
-
-      int idx_i_is_colon = idx_i.is_colon ();
-      int idx_j_is_colon = idx_j.is_colon ();
-
-      if (lhs_nr == 0 && lhs_nc == 0)
-	{
-	  if (idx_i_is_colon)
-	    n = rhs_nr;
-
-	  if (idx_j_is_colon)
-	    m = rhs_nc;
-	}
-
-      if (idx_i && idx_j)
-        {
-          if (rhs_is_scalar && n >= 0 && m >= 0)
-            {
-              // No need to do anything if either of the indices
-              // are empty.
-
-              if (n > 0 && m > 0)
-                {
-                  lhs.make_unique ();
-
-                  MAYBE_RESIZE_LHS;
-
-                  RT scalar = xrhs.elem (0, 0);
-
-                  for (octave_idx_type j = 0; j < m; j++)
-                    {
-                      octave_idx_type jj = idx_j.elem (j);
-                      for (octave_idx_type i = 0; i < n; i++)
-                        {
-                          octave_idx_type ii = idx_i.elem (i);
-                          lhs.xelem (ii, jj) = scalar;
-                        }
-                    }
-                }
-            }
-          else if ((n == 1 || m == 1)
-                   && (rhs_nr == 1 || rhs_nc == 1)
-                   && n * m == rhs_nr * rhs_nc)
-            {
-              lhs.make_unique ();
-
-              MAYBE_RESIZE_LHS;
-
-              if (n > 0 && m > 0)
-                {
-                  octave_idx_type k = 0;
-
-                  for (octave_idx_type j = 0; j < m; j++)
-                    {
-                      octave_idx_type jj = idx_j.elem (j);
-                      for (octave_idx_type i = 0; i < n; i++)
-                        {
-                          octave_idx_type ii = idx_i.elem (i);
-                          lhs.xelem (ii, jj) = xrhs.elem (k++);
-                        }
-                    }
-                }
-            }
-          else if (n == rhs_nr && m == rhs_nc)
-            {
-              lhs.make_unique ();
-
-              MAYBE_RESIZE_LHS;
-
-              if (n > 0 && m > 0)
-                {
-                  for (octave_idx_type j = 0; j < m; j++)
-                    {
-                      octave_idx_type jj = idx_j.elem (j);
-                      for (octave_idx_type i = 0; i < n; i++)
-                        {
-                          octave_idx_type ii = idx_i.elem (i);
-                          lhs.xelem (ii, jj) = xrhs.elem (i, j);
-                        }
-                    }
-                }
-            }
-          else if (n == 0 && m == 0)
-            {
-              if (! (rhs_is_scalar || (rhs_nr == 0 || rhs_nc == 0)))
-                {
-                  lhs.clear_index ();
-
-                  (*current_liboctave_error_handler)
-                    ("A([], []) = X: X must be an empty matrix or a scalar");
-
-                  retval = 0;
-                }
-            }
-          else
-            {
-              lhs.clear_index ();
-
-              (*current_liboctave_error_handler)
-                ("A(I, J) = X: X must be a scalar or the number of elements in I must match the number of rows in X and the number of elements in J must match the number of columns in X");
-
-              retval = 0;
-            }
-        }
-      // idx_vector::freeze() printed an error message for us.
-    }
-  else if (n_idx == 1)
-    {
-      int lhs_is_empty = lhs_nr == 0 || lhs_nc == 0;
-
-      if (lhs_is_empty || (lhs_nr == 1 && lhs_nc == 1))
-	{
-	  octave_idx_type lhs_len = lhs.length ();
-
-	  idx_i.freeze (lhs_len, 0, true);
-
-	  if (idx_i)
-            {
-              if (lhs_is_empty
-                  && idx_i.is_colon ()
-                  && ! (rhs_nr == 1 || rhs_nc == 1))
-                {
-                  (*current_liboctave_warning_with_id_handler)
-                    ("Octave:fortran-indexing",
-                     "A(:) = X: X is not a vector or scalar");
-                }
-              else
-                {
-                  octave_idx_type idx_nr = idx_i.orig_rows ();
-                  octave_idx_type idx_nc = idx_i.orig_columns ();
-
-                  if (! (rhs_nr == idx_nr && rhs_nc == idx_nc))
-                    (*current_liboctave_warning_with_id_handler)
-                      ("Octave:fortran-indexing",
-                       "A(I) = X: X does not have same shape as I");
-                }
-
-              if (assign1 (lhs, xrhs, rfv))
-                {
-                  octave_idx_type len = lhs.length ();
-
-                  if (len > 0)
-                    {
-                      // The following behavior is much simplified
-                      // over previous versions of Octave.  It
-                      // seems to be compatible with Matlab.
-
-                      lhs.dimensions = dim_vector (1, lhs.length ());
-                    }
-                }
-              else
-                retval = 0;
-            }
-	  // idx_vector::freeze() printed an error message for us.
-	}
-      else if (lhs_nr == 1)
-	{
-	  idx_i.freeze (lhs_nc, "vector", true);
-
-	  if (idx_i)
-            {
-              if (assign1 (lhs, xrhs, rfv))
-                lhs.dimensions = dim_vector (1, lhs.length ());
-              else
-                retval = 0;
-            }
-	  // idx_vector::freeze() printed an error message for us.
-	}
-      else if (lhs_nc == 1)
-	{
-	  idx_i.freeze (lhs_nr, "vector", true);
-
-	  if (idx_i)
-	    {
-              if (assign1 (lhs, xrhs, rfv))
-                lhs.dimensions = dim_vector (lhs.length (), 1);
-              else
-                retval = 0;
-	    }
-	  // idx_vector::freeze() printed an error message for us.
-	}
-      else
-	{
-	  if (! idx_i.is_colon ())
-	    (*current_liboctave_warning_with_id_handler)
-	      ("Octave:fortran-indexing", "single index used for matrix");
-
-	  octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix");
-
-	  if (idx_i)
-	    {
-	      if (len == 0)
-		{
-		  if (! (rhs_is_scalar || (rhs_nr == 0 || rhs_nc == 0)))
-		    {
-		      lhs.clear_index ();
-
-		      (*current_liboctave_error_handler)
-			("A([]) = X: X must be an empty matrix or scalar");
-
-		      retval = 0;
-		    }
-		}
-	      else if (len == rhs_nr * rhs_nc)
-		{
-		  lhs.make_unique ();
-
-		  if (idx_i.is_colon ())
-		    {
-		      for (octave_idx_type i = 0; i < len; i++)
-			lhs.xelem (i) = xrhs.elem (i);
-		    }
-		  else
-		    {
-		      for (octave_idx_type i = 0; i < len; i++)
-			{
-			  octave_idx_type ii = idx_i.elem (i);
-			  lhs.xelem (ii) = xrhs.elem (i);
-			}
-		    }
-		}
-	      else if (rhs_is_scalar)
-		{
-		  lhs.make_unique ();
-
-		  RT scalar = rhs.elem (0, 0);
-
-		  if (idx_i.is_colon ())
-		    {
-		      for (octave_idx_type i = 0; i < len; i++)
-			lhs.xelem (i) = scalar;
-		    }
-		  else
-		    {
-		      for (octave_idx_type i = 0; i < len; i++)
-			{
-			  octave_idx_type ii = idx_i.elem (i);
-			  lhs.xelem (ii) = scalar;
-			}
-		    }
-		}
-	      else
-		{
-		  lhs.clear_index ();
-
-		  (*current_liboctave_error_handler)
-      ("A(I) = X: X must be a scalar or a matrix with the same size as I");
-
-		  retval = 0;
-		}
-	    }
-	  // idx_vector::freeze() printed an error message for us.
-	}
-    }
-  else
-    {
-      (*current_liboctave_error_handler)
-	("invalid number of indices for matrix expression");
-
-      retval = 0;
-    }
-
-  lhs.clear_index ();
-
-  return retval;
-}
-
-template <class LT, class RT>
-int
-assignN (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv)
-{
-  int retval = 1;
-
-  dim_vector rhs_dims = rhs.dims ();
-
-  octave_idx_type rhs_dims_len = rhs_dims.length ();
-
-  bool rhs_is_scalar = is_scalar (rhs_dims);
-
-  int n_idx = lhs.index_count ();
-
-  idx_vector *idx_vex = lhs.get_idx ();
-
-  Array<idx_vector> idx = conv_to_array (idx_vex, n_idx);
-
-  if (n_idx == 0)
-    {
-      lhs.clear_index ();
-
-      (*current_liboctave_error_handler)
-	("invalid number of indices for matrix expression");
-
-      retval = 0;
-    }
-  else if (n_idx == 1)
-    {
-      idx_vector iidx = idx(0);
-      int iidx_is_colon = iidx.is_colon ();
-
-      if (! iidx_is_colon)
-	(*current_liboctave_warning_with_id_handler)
-	  ("Octave:fortran-indexing", "single index used for N-d array");
-
-      octave_idx_type lhs_len = lhs.length ();
-
-      octave_idx_type len = iidx.freeze (lhs_len, "N-d arrray");
-
-      if (iidx)
-	{
-	  if (len == 0)
-	    {
-	      if (! (rhs_dims.all_ones () || rhs_dims.any_zero ()))
-		{
-		  lhs.clear_index ();
-
-		  (*current_liboctave_error_handler)
-		    ("A([]) = X: X must be an empty matrix or scalar");
-
-		  retval = 0;
-		}
-	    }
-	  else if (len == rhs.length ())
-	    {
-	      lhs.make_unique ();
-
-	      if (iidx_is_colon)
-		{
-		  for (octave_idx_type i = 0; i < len; i++)
-		    lhs.xelem (i) = rhs.elem (i);
-		}
-	      else
-		{
-		  for (octave_idx_type i = 0; i < len; i++)
-		    {
-		      octave_idx_type ii = iidx.elem (i);
-
-		      lhs.xelem (ii) = rhs.elem (i);
-		    }
-		}
-	    }
-	  else if (rhs_is_scalar)
-	    {
-	      RT scalar = rhs.elem (0);
-
-	      lhs.make_unique ();
-
-	      if (iidx_is_colon)
-		{
-		  for (octave_idx_type i = 0; i < len; i++)
-		    lhs.xelem (i) = scalar;
-		}
-	      else
-		{
-		  for (octave_idx_type i = 0; i < len; i++)
-		    {
-		      octave_idx_type ii = iidx.elem (i);
-
-		      lhs.xelem (ii) = scalar;
-		    }
-		}
-	    }
-	  else
-	    {
-	      lhs.clear_index ();
-
-	      (*current_liboctave_error_handler)
-		("A(I) = X: X must be a scalar or a matrix with the same size as I");
-
-	      retval = 0;
-	    }
-
-	  // idx_vector::freeze() printed an error message for us.
-	}
-    }
-  else
-    {
-      // Maybe expand to more dimensions.
-
-      dim_vector lhs_dims = lhs.dims ();
-
-      octave_idx_type lhs_dims_len = lhs_dims.length ();
-
-      dim_vector final_lhs_dims = lhs_dims;
-
-      dim_vector frozen_len;
-
-      octave_idx_type orig_lhs_dims_len = lhs_dims_len;
-
-      bool orig_empty = lhs_dims.all_zero ();
-
-      if (n_idx < lhs_dims_len)
-	{
-	  // Collapse dimensions beyond last index.  Note that we
-	  // delay resizing LHS until we know that the assignment will
-	  // succeed.
-
-	  if (! (idx(n_idx-1).is_colon ()))
-	    (*current_liboctave_warning_with_id_handler)
-	      ("Octave:fortran-indexing",
-	       "fewer indices than dimensions for N-d array");
-
-	  for (int i = n_idx; i < lhs_dims_len; i++)
-	    lhs_dims(n_idx-1) *= lhs_dims(i);
-
-	  lhs_dims.resize (n_idx);
-
-	  lhs_dims_len = lhs_dims.length ();
-	}
-
-      // Resize.
-
-      dim_vector new_dims;
-      new_dims.resize (n_idx);
-
-      if (orig_empty)
-	{
-	  if (rhs_is_scalar)
-	    {
-	      for (int i = 0; i < n_idx; i++)
-		{
-		  if (idx(i).is_colon ())
-		    new_dims(i) = 1;
-		  else
-		    new_dims(i) = idx(i).orig_empty () ? 0 : idx(i).max () + 1;
-		}
-	    }
-	  else if (is_vector (rhs_dims))
-	    {
-	      int ncolon = 0;
-	      int fcolon = 0;
-	      octave_idx_type new_dims_numel = 1;
-	      int new_dims_vec = 0;
-	      for (int i = 0; i < n_idx; i++)
-		if (idx(i).is_colon ())
-		  {
-		    ncolon ++;
-		    if (ncolon == 1)
-		      fcolon = i;
-		  } 
-		else
-		  {
-		    octave_idx_type cap = idx(i).capacity ();
-		    new_dims_numel *= cap;
-		    if (cap != 1)
-		      new_dims_vec ++;
-		  }
-
-	      if (ncolon == n_idx)
-		{
-		  new_dims = rhs_dims;
-		  new_dims.resize (n_idx);
-		  for (int i = rhs_dims_len; i < n_idx; i++)
-		    new_dims (i) = 1;
-		}
-	      else
-		{
-		  octave_idx_type rhs_dims_numel = rhs_dims.numel ();
-	      	      
-		  for (int i = 0; i < n_idx; i++)
-		    new_dims(i) = idx(i).orig_empty () ? 0 : idx(i).max () + 1;
-
-		  if (new_dims_numel != rhs_dims_numel && 
-		      ncolon > 0 && new_dims_numel == 1)
-		    {
-		      if (ncolon == rhs_dims_len)
-			{
-			  int k = 0;
-			  for (int i = 0; i < n_idx; i++)
-			    if (idx(i).is_colon ())
-			      new_dims (i) = rhs_dims (k++);
-			}
-		      else
-			new_dims (fcolon) = rhs_dims_numel;
-		    }
-		  else if (new_dims_numel != rhs_dims_numel || new_dims_vec > 1)
-		    {
-		      lhs.clear_index ();
-
-		      (*current_liboctave_error_handler)
-			("A(IDX-LIST) = RHS: mismatched index and RHS dimension");
-		      return retval;
-		    }
-		}
-	    }
-	  else
-	    {
-	      int k = 0;
-	      for (int i = 0; i < n_idx; i++)
-		{
-		  if (idx(i).is_colon ())
-		    {
-		      if (k < rhs_dims_len)
-			new_dims(i) = rhs_dims(k++);
-		      else
-			new_dims(i) = 1;
-		    }
-		  else
-		    {
-		      octave_idx_type nelem = idx(i).capacity ();
-
-		      if (nelem >= 1 
-			  && (k < rhs_dims_len && nelem == rhs_dims(k)))
-			k++;
-		      else if (nelem != 1)
-			{
-			  lhs.clear_index ();
-
-			  (*current_liboctave_error_handler)
-			    ("A(IDX-LIST) = RHS: mismatched index and RHS dimension");
-			  return retval;
-			}
-		      new_dims(i) = idx(i).orig_empty () ? 0 : 
-			idx(i).max () + 1;
-		    }
-		}
-	    }
-	}
-      else
-	{
-	  for (int i = 0; i < n_idx; i++)
-	    {
-	      // We didn't start out with all zero dimensions, so if
-	      // index is a colon, it refers to the current LHS
-	      // dimension.  Otherwise, it is OK to enlarge to a
-	      // dimension given by the largest index, but if that
-	      // index is a colon the new dimension is singleton.
-
-	      if (i < lhs_dims_len
-		  && (idx(i).is_colon ()
-		      || idx(i).orig_empty ()
-		      || idx(i).max () < lhs_dims(i)))
-		new_dims(i) = lhs_dims(i);
-	      else if (! idx(i).is_colon ())
-		new_dims(i) = idx(i).max () + 1;
-	      else
-		new_dims(i) = 1;
-	    }
-	}
-
-      if (retval != 0)
-	{
-	  if (! orig_empty
-	      && n_idx < orig_lhs_dims_len
-	      && new_dims(n_idx-1) != lhs_dims(n_idx-1))
-	    {
-	      // We reshaped and the last dimension changed.  This has to
-	      // be an error, because we don't know how to undo that
-	      // later...
-
-	      lhs.clear_index ();
-
-	      (*current_liboctave_error_handler)
-		("array index %d (= %d) for assignment requires invalid resizing operation",
-		 n_idx, new_dims(n_idx-1));
-
-	      retval = 0;
-	    }
-	  else
-	    {
-	      // Determine final dimensions for LHS and reset the
-	      // current size of the LHS.  Note that we delay actually
-	      // resizing LHS until we know that the assignment will
-	      // succeed.
-
-	      if (n_idx < orig_lhs_dims_len)
-		{
-		  for (int i = 0; i < n_idx-1; i++)
-		    final_lhs_dims(i) = new_dims(i);
-		}
-	      else
-		final_lhs_dims = new_dims;
-
-	      lhs_dims_len = new_dims.length ();
-
-	      frozen_len = freeze (idx, new_dims, true);
-
-	      for (int i = 0; i < idx.length (); i++)
-		{
-		  if (! idx(i))
-		    {
-		      retval = 0;
-		      lhs.clear_index ();
-		      return retval;
-		    }
-		}
-
-	      if (rhs_is_scalar)
-		{
-		  lhs.make_unique ();
-
-		  if (n_idx < orig_lhs_dims_len)
-		    lhs = lhs.reshape (lhs_dims);
-
-		  lhs.resize_and_fill (new_dims, rfv);
-
-		  if  (! final_lhs_dims.any_zero ())
-		    {
-		      RT scalar = rhs.elem (0);
-
-		      if (n_idx == 1)
-			{
-			  idx_vector iidx = idx(0);
-
-			  octave_idx_type len = frozen_len(0);
-
-			  if (iidx.is_colon ())
-			    {
-			      for (octave_idx_type i = 0; i < len; i++)
-				lhs.xelem (i) = scalar;
-			    }
-			  else
-			    {
-			      for (octave_idx_type i = 0; i < len; i++)
-				{
-				  octave_idx_type ii = iidx.elem (i);
-
-				  lhs.xelem (ii) = scalar;
-				}
-			    }
-			}
-		      else if (lhs_dims_len == 2 && n_idx == 2)
-			{
-			  idx_vector idx_i = idx(0);
-			  idx_vector idx_j = idx(1);
-
-			  octave_idx_type i_len = frozen_len(0);
-			  octave_idx_type j_len = frozen_len(1);
-
-			  if (idx_i.is_colon())
-			    {
-			      for (octave_idx_type j = 0; j < j_len; j++)
-				{
-				  octave_idx_type off = new_dims (0) *
-				    idx_j.elem (j);
-				  for (octave_idx_type i = 0; i < i_len; i++)
-				    lhs.xelem (i + off) = scalar;
-				}
-			    }
-			  else
-			    {
-			      for (octave_idx_type j = 0; j < j_len; j++)
-				{
-				  octave_idx_type off = new_dims (0) *
-				    idx_j.elem (j);
-				  for (octave_idx_type i = 0; i < i_len; i++)
-				    {
-				      octave_idx_type ii = idx_i.elem (i);
-				      lhs.xelem (ii + off) = scalar;
-				    }
-				}
-			    }
-			}
-		      else
-			{
-			  octave_idx_type n = Array<LT>::get_size (frozen_len);
-
-			  Array<octave_idx_type> result_idx (lhs_dims_len, 0);
-
-			  for (octave_idx_type i = 0; i < n; i++)
-			    {
-			      Array<octave_idx_type> elt_idx = get_elt_idx (idx, result_idx);
-
-			      lhs.xelem (elt_idx) = scalar;
-
-			      increment_index (result_idx, frozen_len);
-			    }
-			}
-		    }
-		}
-	      else
-		{
-		  // RHS is matrix or higher dimension.
-
-		  octave_idx_type n = Array<LT>::get_size (frozen_len);
-
-		  if (n != rhs.numel ())
-		    {
-		      lhs.clear_index ();
-
-		      (*current_liboctave_error_handler)
-			("A(IDX-LIST) = X: X must be a scalar or size of X must equal number of elements indexed by IDX-LIST");
-
-			  retval = 0;
-		    }
-		  else
-		    {
-		      lhs.make_unique ();
-
-		      if (n_idx < orig_lhs_dims_len)
-			lhs = lhs.reshape (lhs_dims);
-
-		      lhs.resize_and_fill (new_dims, rfv);
-
-		      if  (! final_lhs_dims.any_zero ())
-			{
-			  if (n_idx == 1)
-			    {
-			      idx_vector iidx = idx(0);
-
-			      octave_idx_type len = frozen_len(0);
-
-			      if (iidx.is_colon ())
-				{
-				  for (octave_idx_type i = 0; i < len; i++)
-				    lhs.xelem (i) = rhs.elem (i);
-				}
-			      else
-				{
-				  for (octave_idx_type i = 0; i < len; i++)
-				    {
-				      octave_idx_type ii = iidx.elem (i);
-
-				      lhs.xelem (ii) = rhs.elem (i);
-				    }
-				}
-			    }
-			  else if (lhs_dims_len == 2 && n_idx == 2)
-			    {
-			      idx_vector idx_i = idx(0);
-			      idx_vector idx_j = idx(1);
-
-			      octave_idx_type i_len = frozen_len(0);
-			      octave_idx_type j_len = frozen_len(1);
-			      octave_idx_type k = 0;
-
-			      if (idx_i.is_colon())
-				{
-				  for (octave_idx_type j = 0; j < j_len; j++)
-				    {
-				      octave_idx_type off = new_dims (0) * 
-					idx_j.elem (j);
-				      for (octave_idx_type i = 0; 
-					   i < i_len; i++)
-					lhs.xelem (i + off) = rhs.elem (k++);
-				    }
-				}
-			      else
-				{
-				  for (octave_idx_type j = 0; j < j_len; j++)
-				    {
-				      octave_idx_type off = new_dims (0) * 
-					idx_j.elem (j);
-				      for (octave_idx_type i = 0; i < i_len; i++)
-					{
-					  octave_idx_type ii = idx_i.elem (i);
-					  lhs.xelem (ii + off) = rhs.elem (k++);
-					}
-				    }
-				}
-
-			    }
-			  else
-			    {
-			      n = Array<LT>::get_size (frozen_len);
-
-			      Array<octave_idx_type> result_idx (lhs_dims_len, 0);
-
-			      for (octave_idx_type i = 0; i < n; i++)
-				{
-				  Array<octave_idx_type> elt_idx = get_elt_idx (idx, result_idx);
-
-				  lhs.xelem (elt_idx) = rhs.elem (i);
-
-				  increment_index (result_idx, frozen_len);
-				}
-			    }
-			}
-		    }
-		}
-	    }
-	}
-
-      lhs.clear_index ();
-
-      if (retval != 0)
-	lhs = lhs.reshape (final_lhs_dims);
-    }
-
-  if (retval != 0)
-    lhs.chop_trailing_singletons ();
-
-  lhs.clear_index ();
-
-  return retval;
-}
-
 template <class T>
 void
 Array<T>::print_info (std::ostream& os, const std::string& prefix) const

--- a/liboctave/Array.h	Thu Oct 30 00:13:05 2008 +0100
+++ b/liboctave/Array.h	Mon Oct 20 16:54:28 2008 +0200
@@ -3,6 +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@gmail.com>
 
 This file is part of Octave.
 
@@ -31,23 +32,15 @@
 #include <iostream>
 
 #include "dim-vector.h"
+#include "idx-vector.h"
 #include "lo-utils.h"
 #include "oct-sort.h"
 #include "quit.h"
 
-class idx_vector;
-
 // One dimensional array class.  Handles the reference counting for
 // all the derived classes.
 
 template <class T>
-T
-resize_fill_value (const T& x)
-{
-  return x;
-}
-
-template <class T>
 class
 Array
 {
@@ -148,16 +141,12 @@
 
 protected:
 
-  mutable idx_vector *idx;
-  mutable int idx_count;
-
   Array (T *d, octave_idx_type n)
-    : rep (new typename Array<T>::ArrayRep (d, n)), dimensions (n),
-      idx (0), idx_count (0) { }
+    : rep (new typename Array<T>::ArrayRep (d, n)), dimensions (n) { }
 
   Array (T *d, const dim_vector& dv)
     : rep (new typename Array<T>::ArrayRep (d, get_size (dv))),
-      dimensions (dv), idx (0), idx_count (0) { }
+      dimensions (dv) { }
 
 private:
 
@@ -184,16 +173,13 @@
 public:
 
   Array (void)
-    : rep (nil_rep ()), dimensions (),
-      idx (0), idx_count (0) { rep->count++; }
+    : rep (nil_rep ()), dimensions () { rep->count++; }
 
   explicit Array (octave_idx_type n)
-    : rep (new typename Array<T>::ArrayRep (n)), dimensions (n),
-      idx (0), idx_count (0) { }
+    : rep (new typename Array<T>::ArrayRep (n)), dimensions (n) { }
 
   explicit Array (octave_idx_type n, const T& val)
-    : rep (new typename Array<T>::ArrayRep (n)), dimensions (n),
-      idx (0), idx_count (0)
+    : rep (new typename Array<T>::ArrayRep (n)), dimensions (n)
     {
       fill (val);
     }
@@ -202,13 +188,13 @@
   template <class U>
   Array (const Array<U>& a)
     : rep (new typename Array<T>::ArrayRep (coerce (a.data (), a.length ()), a.length ())),
-      dimensions (a.dimensions), idx (0), idx_count (0)
+      dimensions (a.dimensions)
     {
     }
 
   // No type conversion case.
   Array (const Array<T>& a)
-    : rep (a.rep), dimensions (a.dimensions), idx (0), idx_count (0)
+    : rep (a.rep), dimensions (a.dimensions)
     {
       rep->count++;
     }
@@ -217,11 +203,11 @@
 
   Array (const dim_vector& dv)
     : rep (new typename Array<T>::ArrayRep (get_size (dv))),
-      dimensions (dv), idx (0), idx_count (0) { }
+      dimensions (dv) { }
 
   Array (const dim_vector& dv, const T& val)
     : rep (new typename Array<T>::ArrayRep (get_size (dv))),
-      dimensions (dv), idx (0), idx_count (0)
+      dimensions (dv)
     {
       fill (val);
     }
@@ -419,43 +405,6 @@
   Array<T> ipermute (const Array<octave_idx_type>& vec) const
     { return permute (vec, true); }
 
-  void resize_no_fill (octave_idx_type n);
-  void resize_and_fill (octave_idx_type n, const T& val);
-
-  // !!! WARNING !!! -- the following resize_no_fill and
-  // resize_and_fill functions are public because template friends
-  // don't work properly with versions of gcc earlier than 3.3.  You
-  // should use these functions only in classes that are derived
-  // from Array<T>.
-
-  // protected:
-
-  void resize_no_fill (octave_idx_type r, octave_idx_type c);
-  void resize_and_fill (octave_idx_type r, octave_idx_type c, const T& val);
-
-  void resize_no_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p);
-  void resize_and_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val);
-
-  void resize_no_fill (const dim_vector& dv);
-  void resize_and_fill (const dim_vector& dv, const T& val);
-
-public:
-
-  void resize (octave_idx_type n) { resize_no_fill (n); }
-
-  void resize (octave_idx_type n, const T& val) { resize_and_fill (n, val); }
-
-  void resize (const dim_vector& dv) { resize_no_fill (dv); }
-
-  void resize (const dim_vector& dv, const T& val)
-    { resize_and_fill (dv, val); }
-
-  Array<T>& insert (const Array<T>& a, octave_idx_type r, octave_idx_type c);
-  Array<T>& insert2 (const Array<T>& a, octave_idx_type r, octave_idx_type c);
-  Array<T>& insertN (const Array<T>& a, octave_idx_type r, octave_idx_type c);
-
-  Array<T>& insert (const Array<T>& a, const Array<octave_idx_type>& idx);
-
   bool is_square (void) const { return (dim1 () == dim2 ()); }
 
   bool is_empty (void) const { return numel () == 0; }
@@ -482,47 +431,98 @@
 
   void maybe_delete_dims (void);
 
-  void clear_index (void) const;
+  // Indexing without resizing.
+
+  Array<T> index (const idx_vector& i) const;
 
-  void set_index (const idx_vector& i) const;
+  Array<T> index (const idx_vector& i, const idx_vector& j) const;
+
+  Array<T> index (const Array<idx_vector>& ia) const;
 
-  int index_count (void) const { return idx_count; }
+  static T resize_fill_value (); 
+
+  // Resizing (with fill).
 
-  idx_vector *get_idx (void) const { return idx; }
+  void resize_fill (octave_idx_type n, const T& rfv);
 
-  void maybe_delete_elements (idx_vector& i);
+  void resize_fill (octave_idx_type nr, octave_idx_type nc, const T& rfv);
+
+  void resize_fill (const dim_vector& dv, const T& rfv);
 
-  void maybe_delete_elements_1 (idx_vector& i);
+  // Resizing with default fill.
+  // Rationale: 
+  // These use the default fill value rather than leaving memory uninitialized.
+  // Resizing without fill leaves the resulting array in a rather weird state,
+  // where part of the data is initialized an part isn't.
 
-  void maybe_delete_elements_2 (idx_vector& i);
+  void resize (octave_idx_type n)
+    { resize_fill (n, resize_fill_value ()); }
 
-  void maybe_delete_elements (idx_vector& i, idx_vector& j);
+  // FIXME: This method cannot be defined here because it would clash with
+  //   void resize (octave_idx_type, const T&)
+  // (these become undistinguishable when T = octave_idx_type). 
+  // In the future, I think the resize (.., const T& rfv) overloads should go away
+  // in favor of using resize_fill.
 
-  void maybe_delete_elements (idx_vector& i, idx_vector& j, idx_vector& k);
+  // void resize (octave_idx_type nr, octave_idx_type nc)
+  //  { resize_fill (nr, nc, resize_fill_value ()); }
 
-  void maybe_delete_elements (Array<idx_vector>& ra_idx);
+  void resize (dim_vector dv)
+    { resize_fill (dv, resize_fill_value ()); }
+
+  // FIXME: these are here for backward compatibility. They should go away in favor
+  // of using resize_fill directly.
+  void resize (octave_idx_type n, const T& rfv)
+    { resize_fill (n, static_cast<T> (rfv)); }
 
-  Array<T> value (void) const;
+  void resize (octave_idx_type nr, octave_idx_type nc, const T& rfv)
+    { resize_fill (nr, nc, rfv); }
 
-  Array<T> index (idx_vector& i, int resize_ok = 0,
-		  const T& rfv = resize_fill_value (T ())) const;
+  void resize (dim_vector dv, const T& rfv)
+    { resize_fill (dv, rfv); }
+
+  // Indexing with possible resizing and fill
+  // FIXME: This is really a corner case, that should better be handled
+  // directly in liboctinterp.
 
-  Array<T> index1 (idx_vector& i, int resize_ok = 0,
-		   const T& rfv = resize_fill_value (T ())) const;
+  Array<T> index (const idx_vector& i, bool resize_ok,
+                  const T& rfv = resize_fill_value ()) const;
+
+  Array<T> index (const idx_vector& i, const idx_vector& j, 
+                  bool resize_ok, const T& rfv = resize_fill_value ()) const;
 
-  Array<T> index2 (idx_vector& i, int resize_ok = 0,
-		   const T& rfv = resize_fill_value (T ())) const;
+  Array<T> index (const Array<idx_vector>& ia,
+                  bool resize_ok, const T& rfv = resize_fill_value ()) const;
+
+  // Indexed assignment (always with resize & fill).
+
+  void assign (const idx_vector& i, const Array<T>& rhs, 
+               const T& rfv = resize_fill_value ());
 
-  Array<T> indexN (idx_vector& i, int resize_ok = 0,
-		   const T& rfv = resize_fill_value (T ())) const;
+  void assign (const idx_vector& i, const idx_vector& j, const Array<T>& rhs,
+               const T& rfv = resize_fill_value ());
+
+  void assign (const Array<idx_vector>& ia, const Array<T>& rhs,
+               const T& rfv = resize_fill_value ());
+
+  // Deleting elements.
+
+  // A(I) = [] (with a single subscript)
+  void delete_elements (const idx_vector& i);
 
-  Array<T> index (idx_vector& i, idx_vector& j, int resize_ok = 0,
-		  const T& rfv = resize_fill_value (T ())) const;
+  // A(:,...,I,...,:) = [] (>= 2 subscripts, one of them is non-colon)
+  void delete_elements (int dim, const idx_vector& i);
+
+  // Dispatcher to the above two.
+  void delete_elements (const Array<idx_vector>& ia);
 
-  Array<T> index (Array<idx_vector>& ra_idx, int resize_ok = 0,
-		  const T& rfv = resize_fill_value (T ())) const;
+  // FIXME: are these required? What exactly are they supposed to do?.
 
-  //  static T resize_fill_value (void) { return T (); }
+  Array<T>& insert (const Array<T>& a, octave_idx_type r, octave_idx_type c);
+  Array<T>& insert2 (const Array<T>& a, octave_idx_type r, octave_idx_type c);
+  Array<T>& insertN (const Array<T>& a, octave_idx_type r, octave_idx_type c);
+
+  Array<T>& insert (const Array<T>& a, const Array<octave_idx_type>& idx);
 
   void print_info (std::ostream& os, const std::string& prefix) const;
 
@@ -558,48 +558,18 @@
   }
 };
 
-// NOTE: these functions should be friends of the Array<T> class and
-// Array<T>::dimensions should be protected, not public, but we can't
-// do that because of bugs in gcc prior to 3.3.
-
-template <class LT, class RT>
-/* friend */ int
-assign (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv);
-
-template <class LT, class RT>
-/* friend */ int
-assign1 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv);
-
-template <class LT, class RT>
-/* friend */ int
-assign2 (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv);
-
-template <class LT, class RT>
-/* friend */ int
-assignN (Array<LT>& lhs, const Array<RT>& rhs, const LT& rfv);
+#define INSTANTIATE_ARRAY(T, API) \
+  template class API Array<T>
 
-template <class LT, class RT>
-int
-assign (Array<LT>& lhs, const Array<RT>& rhs)
-{
-  return assign (lhs, rhs, resize_fill_value (LT ()));
-}
+// FIXME: These are here for compatibility. In current implementation, only
+// homogeneous array assignments are actually instantiated. I think heterogeneous
+// indexed assignments are rare enough to be implemented via conversion first.
+// This decision may still be revised, that's why these macros stay here.
+#define INSTANTIATE_ARRAY_AND_ASSIGN(T, API) \
+  INSTANTIATE_ARRAY(T, API)
 
-#define INSTANTIATE_ARRAY_ASSIGN(LT, RT, API) \
-  template API int assign (Array<LT>&, const Array<RT>&, const LT&); \
-  template API int assign1 (Array<LT>&, const Array<RT>&, const LT&); \
-  template API int assign2 (Array<LT>&, const Array<RT>&, const LT&); \
-  template API int assignN (Array<LT>&, const Array<RT>&, const LT&); \
-  template API int assign (Array<LT>&, const Array<RT>&)
-
-
-#define INSTANTIATE_ARRAY(T, API) \
-  template class API Array<T>; \
-  template API T resize_fill_value (const T&); \
-
-#define INSTANTIATE_ARRAY_AND_ASSIGN(T, API) \
-  INSTANTIATE_ARRAY (T, API); \
-  INSTANTIATE_ARRAY_ASSIGN (T, T, API)
+#define INSTANTIATE_ARRAY_ASSIGN(LT, RT, API)
+  // do nothing
 
 #define INSTANTIATE_ARRAY_SORT(T) \
   template class octave_sort<T>; \