octave-nkf: liboctave/Array.cc comparison

comparison liboctave/Array.cc @ 8290:7cbe01c21986

improve dense array indexing

author	Jaroslav Hajek <highegg@gmail.com>
date	Mon, 20 Oct 2008 16:54:28 +0200
parents	6c08e3921d3e
children	f2e050b62199

comparison

equal deleted inserted replaced

-:ac7f334d9652
+:7cbe01c21986
 // Template array classes
 /*
 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.
 Octave is free software; you can redistribute it and/or modify it
 under the terms of the GNU General Public License as published by the
 #include <climits>
 #include <iostream>
 #include <sstream>
 #include <vector>
+#include <algorithm>
 #include <new>
 #include "Array.h"
 #include "Array-util.h"
 #include "idx-vector.h"
 // One dimensional array class.  Handles the reference counting for
 // all the derived classes.
 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++;
 if (a.numel () < dv.numel ())
 (*current_liboctave_error_handler)
 template <class T>
 Array<T>::~Array (void)
 {
 if (--rep->count <= 0)
 delete rep;
-delete [] idx;
 }
 template <class T>
 Array<T>&
 Array<T>::operator = (const Array<T>& a)
 rep = a.rep;
 rep->count++;
 dimensions = a.dimensions;
-delete [] idx;
-idx_count = 0;
-idx = 0;
 }
 return *this;
 }
 retval.chop_trailing_singletons ();
 return retval;
 }
+// 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
+{
+octave_idx_type *dim, *cdim;
+idx_vector *idx;
+int top;
+public:
+rec_index_helper (const dim_vector& dv, const Array<idx_vector>& ia)
+{
+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;
+}
+// Recursive N-d indexed assignment
+template <class T>
+void do_fill (const T& val, T *dest, int lev) const
+{
+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);
+}
+}
+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>
+Array<T>
+Array<T>::index (const idx_vector& i) const
+{
+octave_idx_type n = numel ();
+Array<T> retval;
+if (i.is_colon ())
+{
+// A(:) produces a shallow copy as a column vector.
+retval.dimensions = dim_vector (n, 1);
+rep->count++;
+retval.rep = rep;
+}
+else if (i.extent (n) != n)
+{
+gripe_index_out_of_range ();
+}
+else
+{
+// 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);
+}
+}
+return retval;
+}
+template <class T>
+Array<T>
+Array<T>::index (const Array<idx_vector>& ia) const
+{
+int ial = ia.length ();
+Array<T> retval;
+// 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)
+{
+// Get dimensions, allowing Fortran indexing in the last dim.
+dim_vector dv = dimensions.redim (ial);
+// 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 ();
+}
+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 ();
+// Don't use resize here to avoid useless initialization for POD types.
+retval = Array<T> (rdv);
+// Prepare for recursive indexing
+rec_index_helper rh (dv, ia);
+// Do it.
+rh.index (data (), retval.fortran_vec ());
+}
+}
+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 n)
+Array<T>::resize_fill (octave_idx_type n, const T& rfv)
 {
-if (n < 0)
+if (n >= 0 && ndims () == 2)
+{
+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 ();
+octave_idx_type n0 = std::min (n, nx), n1 = n - n0;
+dest = std::copy (data (), data () + n0, dest);
+std::fill (dest, dest + n1, rfv);
+*this = tmp;
+}
+}
+}
+else
+gripe_invalid_resize ();
+}
+template <class T>
+void
+Array<T>::resize_fill (octave_idx_type r, octave_idx_type c, const T& rfv)
+{
+if (r >= 0 && c >= 0 && ndims () == 2)
+{
+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> ();
+}
+return tmp.index (i);
+}
+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)
+{
+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 (tmp.rows () != rx || tmp.columns () != cx)
+return Array<T> ();
+}
+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>::assign (const idx_vector& i, const Array<T>& rhs, const T& rfv)
+{
+octave_idx_type n = numel (), rhl = rhs.numel ();
+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 ();
+}
+// 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 ());
+}
+}
+}
+else
+gripe_invalid_assignment_size ();
+}
+template <class T>
+void
+Array<T>::assign (const idx_vector& i, const idx_vector& j,
+const Array<T>& rhs, const T& rfv)
+{
+// 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
+{
+rdv(0) = i.extent (dv(0));
+rdv(1) = j.extent (dv(1));
+}
+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 (match)
+{
+// Resize if requested.
+if (rdv != dv)
+{
+resize (rdv, rfv);
+dv = dimensions;
+}
+// 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);
+const T* src = rhs.data ();
+T *dest = fortran_vec ();
+// 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));
+}
+}
+}
+}
+}
+else
+gripe_assignment_dimension_mismatch ();
+}
+template <class T>
+void
+Array<T>::assign (const Array<idx_vector>& ia,
+const Array<T>& rhs, const T& rfv)
+{
+int ial = ia.length ();
+// 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)
+{
+// 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 ();
+}
+// 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);
+// Do it.
+if (isfill)
+rh.fill (rhs(0), fortran_vec ());
+else
+rh.assign (rhs.data (), fortran_vec ());
+}
+}
+}
+else
+gripe_assignment_dimension_mismatch ();
+}
+}
+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)
+{
+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)
-	("can't resize to negative dimension");
+("invalid dimension in delete_elements");
 return;
 }
-if (n == length ())
+octave_idx_type n = dimensions (dim);
-return;
+if (i.is_colon ())
+{
-typename Array<T>::ArrayRep *old_rep = rep;
+*this = Array<T> ();
-const T *old_data = data ();
+}
-octave_idx_type old_len = length ();
+else if (i.extent (n) != n)
+{
-rep = new typename Array<T>::ArrayRep (n);
+gripe_index_out_of_range ();
+}
-dimensions = dim_vector (n);
+else if (i.length (n) != 0)
+{
-if (n > 0 && old_data && old_len > 0)
+octave_idx_type l, u;
-{
-octave_idx_type min_len = old_len < n ? old_len : n;
+if (i.is_cont_range (n, l, u))
+{
-for (octave_idx_type i = 0; i < min_len; i++)
+// Special case deleting a contiguous range.
-	xelem (i) = old_data[i];
+octave_idx_type nd = n + l - u, dl = 1, du = 1;
-}
+dim_vector rdv = dimensions;
+rdv(dim) = nd;
-if (--old_rep->count <= 0)
+for (int k = 0; k < dim; k++) dl *= dimensions(k);
-delete old_rep;
+for (int k = dim + 1; k < ndims (); k++) du *= dimensions(k);
-}
+// Special case deleting a contiguous range.
-template <class T>
+Array<T> tmp = Array<T> (rdv);
-void
+const T *src = data ();
-Array<T>::resize_no_fill (const dim_vector& dv)
+T *dest = tmp.fortran_vec ();
-{
+l *= dl; u *= dl; n *= dl;
-octave_idx_type n = dv.length ();
+for (octave_idx_type k = 0; k < du; k++)
+{
-for (octave_idx_type i = 0; i < n; i++)
+dest = std::copy (src, src + l, dest);
-{
+dest = std::copy (src + u, src + n, dest);
-if (dv(i) < 0)
+src += n;
-	{
+}
-	  (*current_liboctave_error_handler)
-	    ("can't resize to negative dimension");
+*this = tmp;
-	  return;
+}
-	}
+else
-}
+{
+// Use index.
-bool same_size = true;
+Array<idx_vector> ia (ndims (), idx_vector::colon);
+ia (dim) = i.complement (n);
-if (dimensions.length () != n)
+*this = index (ia);
-{
+}
-same_size = false;
+}
 }
-else
-{
+template <class T>
-for (octave_idx_type i = 0; i < n; i++)
+void
-	{
+Array<T>::delete_elements (const Array<idx_vector>& ia)
-	  if (dv(i) != dimensions(i))
+{
-	    {
+if (ia.length () == 1)
-	      same_size = false;
+delete_elements (ia(0));
-	      break;
+else
-	    }
+{
-	}
+int len = ia.length (), k, dim = -1;
-}
+for (k = 0; k < len; k++)
+{
-if (same_size)
+if (! ia(k).is_colon ())
-return;
+{
+if (dim < 0)
-typename Array<T>::ArrayRep *old_rep = rep;
+dim = k;
-const T *old_data = data ();
+else
+break;
-octave_idx_type ts = get_size (dv);
+}
+}
-rep = new typename Array<T>::ArrayRep (ts);
+if (dim < 0)
+{
-dim_vector dv_old = dimensions;
+dim_vector dv = dimensions;
-octave_idx_type  dv_old_orig_len = dv_old.length ();
+dv(0) = 0;
-dimensions = dv;
+*this = Array<T> (dv);
-octave_idx_type ts_old = get_size (dv_old);
+}
+else if (k == len)
-if (ts > 0 && ts_old > 0 && dv_old_orig_len > 0)
+{
-{
+delete_elements (dim, ia(dim));
-Array<octave_idx_type> ra_idx (dimensions.length (), 0);
+}
+else
-if (n > dv_old_orig_len)
+{
-	{
+(*current_liboctave_error_handler)
-	  dv_old.resize (n);
+("A null assignment can only have one non-colon index.");
+}
-	  for (octave_idx_type i = dv_old_orig_len; i < n; i++)
+}
-	    dv_old.elem (i) = 1;
-	}
+}
-for (octave_idx_type i = 0; i < ts; i++)
+// FIXME: Remove these methods or implement them using assign.
-	{
-	  if (index_in_bounds (ra_idx, dv_old))
-	    rep->elem (i) = old_data[get_scalar_idx (ra_idx, dv_old)];
-	  increment_index (ra_idx, dimensions);
-	}
-}
-if (--old_rep->count <= 0)
-delete old_rep;
-}
-template <class T>
-void
-Array<T>::resize_no_fill (octave_idx_type r, octave_idx_type c)
-{
-if (r < 0 || c < 0)
-{
-(*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;
-}
-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 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];
-}
-if (--old_rep->count <= 0)
-delete old_rep;
-}
-template <class T>
-void
-Array<T>::resize_and_fill (octave_idx_type n, const T& val)
-{
-if (n < 0)
-{
-(*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-return;
-}
-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)
-{
-octave_idx_type min_len = old_len < n ? old_len : n;
-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 (--old_rep->count <= 0)
-delete old_rep;
-}
-template <class T>
-void
-Array<T>::resize_and_fill (octave_idx_type r, octave_idx_type c, const T& val)
-{
-if (r < 0 || c < 0)
-{
-(*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-return;
-}
-if (ndims () == 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)
-{
-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 (--old_rep->count <= 0)
-delete old_rep;
-}
-template <class T>
-void
-Array<T>::resize_and_fill (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val)
-{
-if (r < 0 || c < 0 || p < 0)
-{
-(*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-return;
-}
-if (ndims () == 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 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 (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.
-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;
-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;
-}
-template <class T>
-void
-Array<T>::resize_and_fill (const dim_vector& dv, const T& val)
-{
-octave_idx_type n = dv.length ();
-for (octave_idx_type i = 0; i < n; i++)
-{
-if (dv(i) < 0)
-	{
-	  (*current_liboctave_error_handler)
-	    ("can't resize to negative dimension");
-	  return;
-	}
-}
-bool same_size = true;
-if (dimensions.length () != n)
-{
-same_size = false;
-}
-else
-{
-for (octave_idx_type i = 0; i < n; i++)
-	{
-	  if (dv(i) != dimensions(i))
-	    {
-	      same_size = false;
-	      break;
-	    }
-	}
-}
-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);
-dim_vector dv_old = dimensions;
-octave_idx_type dv_old_orig_len = dv_old.length ();
-dimensions = dv;
-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);
-	  for (octave_idx_type i = dv_old_orig_len; i < n; i++)
-	    dv_old.elem (i) = 1;
-	}
-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);
-	}
-}
-else
-for (octave_idx_type i = 0; i < len; i++)
-rep->elem (i) = val;
-if (--old_rep->count <= 0)
-delete old_rep;
-}
 template <class T>
 Array<T>&
 Array<T>::insert (const Array<T>& a, octave_idx_type r, octave_idx_type c)
 {
 (*current_liboctave_error_handler)
 ("Array<T>::insert: invalid indexing operation");
 return *this;
 }
 template <class T>
 Array<T>
 Array<T>::transpose (void) const
 {
 	new_dims(i) = dimensions(i);
 }
 if (nd != new_dims.length ())
 dimensions = new_dims;
-}
-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)
 	    }
 	  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);
 	    }
 	  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);
 	    }
 	}
 }
 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

Mercurial > hg > octave-nkf

comparison liboctave/Array.cc @ 8290:7cbe01c21986