--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -46,6 +46,26 @@
 // all the derived classes.
 
 template <class T>
+void
+Array<T>::make_unique (void)
+{
+  if (rep->count > 1)
+    {
+      --rep->count;
+      rep = new ArrayRep (slice_data, slice_len, true);
+      slice_data = rep->data;
+    }
+  else if (slice_len != rep->len)
+    {
+      // Possibly economize here.
+      ArrayRep *new_rep = new ArrayRep (slice_data, slice_len, true);
+      delete rep;
+      rep = new_rep;
+      slice_data = rep->data;
+    }
+}
+
+template <class T>
 Array<T>::Array (const Array<T>& a, const dim_vector& dv)
   : rep (a.rep), dimensions (dv), 
     slice_data (a.slice_data), slice_len (a.slice_len)
@@ -85,6 +105,20 @@
 }
 
 template <class T>
+void
+Array<T>::fill (const T& val)
+{
+  if (rep->count > 1)
+    {
+      --rep->count;
+      rep = new ArrayRep (length (), val);
+      slice_data = rep->data;
+    }
+  else
+    std::fill (slice_data, slice_data + slice_len, val);
+}
+
+template <class T>
 Array<T>
 Array<T>::squeeze (void) const
 {
@@ -131,13 +165,7 @@
 	    }
 	}
 
-      // FIXME -- it would be better if we did not have to do
-      // this, so we could share the data while still having different
-      // dimension vectors.
-
-      retval.make_unique ();
-
-      retval.dimensions = new_dimensions;
+      retval = Array<T> (*this, new_dimensions);
     }
 
   return retval;

--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -97,50 +97,16 @@
 
 public:
 
-  // !!! WARNING !!! -- these should be protected, not public.  You
-  // should not access these methods directly!
-
-  void make_unique (void)
-    {
-      if (rep->count > 1)
-	{
-	  --rep->count;
-	  rep = new ArrayRep (slice_data, slice_len, true);
-          slice_data = rep->data;
-	}
-      else if (slice_len != rep->len)
-        {
-          // Possibly economize here.
-          ArrayRep *new_rep = new ArrayRep (slice_data, slice_len, true);
-          delete rep;
-          rep = new_rep;
-          slice_data = rep->data;
-        }
-      }
-
-  void make_unique (const T& val)
-    {
-      if (rep->count > 1)
-	{
-	  --rep->count;
-	  rep = new ArrayRep (slice_len, val);
-          slice_data = rep->data;
-	}
-      else
-        std::fill (slice_data, slice_data + slice_len, val);
-    }
+  void make_unique (void);
 
   typedef T element_type;
 
-  // !!! WARNING !!! -- these should be protected, not public.  You
-  // should not access these data members directly!
+protected:
 
   typename Array<T>::ArrayRep *rep;
 
   dim_vector dimensions;
 
-protected:
-
   T* slice_data;
   octave_idx_type slice_len;
 
@@ -220,7 +186,7 @@
   template <class U>
   Array (const Array<U>& a)
     : rep (new typename Array<T>::ArrayRep (coerce (a.data (), a.length ()), a.length ())),
-      dimensions (a.dimensions)
+      dimensions (a.dims ())
     {
       slice_data = rep->data;
       slice_len = rep->len;
@@ -260,7 +226,7 @@
 
   Array<T>& operator = (const Array<T>& a);
 
-  void fill (const T& val) { make_unique (val); }
+  void fill (const T& val); 
 
   octave_idx_type capacity (void) const { return slice_len; }
   octave_idx_type length (void) const { return capacity (); }

--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,17 @@
+2009-01-14  Jaroslav Hajek  <highegg@gmail.com>
+	
+	* Array.h (Array<T>::rep, Array<T>::dimensions): Make protected.
+	* Array.cc (Array<T>::make_unique): Move implementation here.
+	(Array<T>::fill): Dtto.
+	* DiagArray2.h (DiagArray2<T>): Reimplement without abusing
+	Array<T> internals.
+	(DiagArray2<T>::operator Array2<T>): New method.
+	* DiagArray2.cc (DiagArray2<T>): Update methods.
+	* MDiagArray2.h (MDiagArray2<T>::operator Array2<T>): Simplify.
+	* PermMatrix.h (PermMatrix): Reimplement without abusing
+	Array<T> internals.
+	* PermMatrix.cc (PermMatrix): Update methods.
+
 2009-01-14  Jaroslav Hajek  <highegg@gmail.com>
 
 	* Array.cc, Array.h (all Array<T> constructors): Handle slice_data and

--- a/liboctave/DiagArray2.cc
+++ b/liboctave/DiagArray2.cc
@@ -37,17 +37,42 @@
 #include "lo-error.h"
 
 template <class T>
+const typename DiagArray2<T>::Proxy& 
+DiagArray2<T>::Proxy::operator = (const T& val) const
+{
+  if (i == j)
+    {
+      if (object)
+        object->set (val, i);
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("invalid assignment to off-diagonal in diagonal array");
+
+  return *this;
+}
+
+template <class T>
+DiagArray2<T>::Proxy::operator T () const
+{
+  if (object && i == j)
+    return object->get (i);
+  else
+    {
+      static T foo;
+      return foo;
+    }
+}
+
+template <class T>
 Array<T>
 DiagArray2<T>::diag (octave_idx_type k) const
 {
   Array<T> d;
 
   if (k == 0)
-    {
-      // The main diagonal is shallow-copied.
-      d = *this;
-      d.dimensions = dim_vector (length ());
-    }
+    // The main diagonal is shallow-copied.
+    d = *this;
   else if (k > 0 && k < cols ())
     d = Array<T> (std::min (cols () - k, rows ()), T ());
   else if (k < 0 && -k < rows ())
@@ -64,7 +89,8 @@
 DiagArray2<T>::transpose (void) const
 {
   DiagArray2<T> retval (*this);
-  retval.dimensions = dim_vector (dim2 (), dim1 ());
+  retval.d1 = d2;
+  retval.d2 = d1;
   return retval;
 }
 
@@ -91,7 +117,20 @@
       (*current_liboctave_error_handler) ("range error in DiagArray2");
       return T ();
     }
-  return (r == c) ? Array<T>::xelem (r) : T (0);
+  return elem (r, c);
+}
+
+template <class T>
+typename DiagArray2<T>::Proxy
+DiagArray2<T>::checkelem (octave_idx_type r, octave_idx_type c) 
+{
+  if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
+    {
+      (*current_liboctave_error_handler) ("range error in DiagArray2");
+      return Proxy (0, r, c);
+    }
+  else
+    return Proxy (this, r, c);
 }
 
 template <class T>
@@ -104,37 +143,16 @@
       return;
     }
 
-  if (r == dim1 () && c == dim2 ())
-    return;
-
-  // FIXME: this is a mess. DiagArray2 really needs a rewrite.
-  typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
-  const T *old_data = this->data ();
-  octave_idx_type old_len = this->length ();
-
-  octave_idx_type new_len = r < c ? r : c;
-
-  Array<T>::rep = new typename Array<T>::ArrayRep (new_len);
-  Array<T>::slice_data = Array<T>::rep->data;
-  Array<T>::slice_len = Array<T>::rep->len;
-
-  this->dimensions = dim_vector (r, c);
-
-  if (old_data && old_len > 0)
+  if (r != dim1 () || c != dim2 ())
     {
-      octave_idx_type min_len = old_len < new_len ? old_len : new_len;
-
-      for (octave_idx_type i = 0; i < min_len; i++)
-	xelem (i, i) = old_data[i];
+      Array<T>::resize (std::min (r, c));
+      d1 = r; d2 = c;
     }
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
 }
 
 template <class T>
 void
-DiagArray2<T>::resize (octave_idx_type r, octave_idx_type c, const T& val)
+DiagArray2<T>::resize_fill (octave_idx_type r, octave_idx_type c, const T& val)
 {
   if (r < 0 || c < 0)
     {
@@ -142,32 +160,21 @@
       return;
     }
 
-  if (r == dim1 () && c == dim2 ())
-    return;
-
-  typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
-  const T *old_data = this->data ();
-  octave_idx_type old_len = this->length ();
-
-  octave_idx_type new_len = r < c ? r : c;
-
-  Array<T>::rep = new typename Array<T>::ArrayRep (new_len);
-
-  this->dimensions = dim_vector (r, c);
+  if (r != dim1 () || c != dim2 ())
+    {
+      Array<T>::resize_fill (std::min (r, c), val);
+      d1 = r; d2 = c;
+    }
+}
 
-  octave_idx_type min_len = old_len < new_len ? old_len : new_len;
+template <class T>
+DiagArray2<T>::operator Array2<T> (void) const
+{
+  Array2<T> result (dim1 (), dim2 ());
+  for (octave_idx_type i = 0, len = length (); i < len; i++)
+    result.xelem (i, i) = dgelem (i);
 
-  if (old_data && old_len > 0)
-    {
-      for (octave_idx_type i = 0; i < min_len; i++)
-	xelem (i, i) = old_data[i];
-    }
-
-  for (octave_idx_type i = min_len; i < new_len; i++)
-    xelem (i, i) = val;
-
-  if (--old_rep->count <= 0)
-    delete old_rep;
+  return result;
 }
 
 /*

--- a/liboctave/DiagArray2.h
+++ b/liboctave/DiagArray2.h
@@ -3,7 +3,7 @@
 
 Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007
               John W. Eaton
-Copyright (C) 2008 Jaroslav Hajek
+Copyright (C) 2008, 2009 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -30,10 +30,10 @@
 #include <cstdlib>
 
 #include "Array.h"
+#include "Array2.h"
 #include "lo-error.h"
 
 // A two-dimensional array with diagonal elements only.
-//
 // Idea and example code for Proxy class and functions from:
 //
 // From: kanze@us-es.sel.de (James Kanze)
@@ -45,13 +45,12 @@
 // James Kanze                             email: kanze@us-es.sel.de
 // GABI Software, Sarl., 8 rue du Faisan, F-67000 Strasbourg, France
 
-// Array<T> is inherited privately because we abuse the dimensions variable
-// for true dimensions. Therefore, the inherited Array<T> object is not a valid
-// Array<T> object, and should not be publicly accessible.
+// Array<T> is inherited privately so that some methods, like index, don't
+// produce unexpected results.
 
 template <class T>
 class
-DiagArray2 : private Array<T>
+DiagArray2 : protected Array<T>
 {
 private:
 
@@ -66,30 +65,9 @@
     Proxy (DiagArray2<T> *ref, octave_idx_type r, octave_idx_type c)
       : i (r), j (c), object (ref) { } 
 
-    const Proxy& operator = (const T& val) const
-      {
-	if (i == j)
-	  {
-	    if (object)
-	      object->set (val, i);
-	  }
-	else
-	  (*current_liboctave_error_handler)
-	    ("invalid assignment to off-diagonal in diagonal array");
+    const Proxy& operator = (const T& val) const;
 
-	return *this;
-      }
-
-    operator T () const
-      {
-	if (object && i == j)
-	  return object->get (i);
-	else
-	  {
-	    static T foo;
-	    return foo;
-	  }
-      }
+    operator T () const;
 
   private:
 
@@ -106,100 +84,92 @@
 
   };
 
-friend class Proxy;
+  friend class Proxy;
 
 protected:
+  octave_idx_type d1, d2;
 
-  DiagArray2 (T *d, octave_idx_type r, octave_idx_type c) : Array<T> (d, r < c ? r : c)
-    { Array<T>::dimensions = dim_vector (r, c); }
+  DiagArray2 (T *d, octave_idx_type r, octave_idx_type c) 
+    : Array<T> (d, std::min (r, c)), d1 (r), d2 (c) { }
 
 public:
 
   typedef T element_type;
 
-  DiagArray2 (void) : Array<T> (dim_vector (0, 0)) { }
+  DiagArray2 (void) 
+    : Array<T> (), d1 (0), d2 (0) { }
 
-  DiagArray2 (octave_idx_type r, octave_idx_type c) : Array<T> (r < c ? r : c)
-    { this->dimensions = dim_vector (r, c); }
-
-  DiagArray2 (octave_idx_type r, octave_idx_type c, const T& val) : Array<T> (r < c ? r : c)
-    {
-      this->dimensions = dim_vector (r, c);
+  DiagArray2 (octave_idx_type r, octave_idx_type c) 
+    : Array<T> (std::min (r, c)), d1 (r), d2 (c) { }
 
-      Array<T>::fill (val);
-    }
+  DiagArray2 (octave_idx_type r, octave_idx_type c, const T& val) 
+    : Array<T> (std::min (r, c), val), d1 (r), d2 (c) { }
 
-  DiagArray2 (const Array<T>& a) : Array<T> (a)
-    { this->dimensions = dim_vector (a.length (), a.length ()); }
+  DiagArray2 (const Array<T>& a) 
+    : Array<T> (a), d1 (a.numel ()), d2 (a.numel ()) { }
 
-  DiagArray2 (const DiagArray2<T>& a) : Array<T> (a)
-    { this->dimensions = a.dims (); }
+  DiagArray2 (const DiagArray2<T>& a) 
+    : Array<T> (a), d1 (a.d1), d2 (a.d2) { }
 
   template <class U>
-  DiagArray2 (const DiagArray2<U>& a) : Array<T> (a.diag ())
-    { this->dimensions = a.dims (); }
+  DiagArray2 (const DiagArray2<U>& a) 
+  : Array<T> (a.diag ()), d1 (a.dim1 ()), d2 (a.dim2 ()) { }
 
   ~DiagArray2 (void) { }
 
   DiagArray2<T>& operator = (const DiagArray2<T>& a)
     {
       if (this != &a)
-	Array<T>::operator = (a);
+        {
+          Array<T>::operator = (a);
+          d1 = a.d1;
+          d2 = a.d2;
+        }
 
       return *this;
     }
 
-
-  octave_idx_type dim1 (void) const { return Array<T>::dimensions(0); }
-  octave_idx_type dim2 (void) const { return Array<T>::dimensions(1); }
+  octave_idx_type dim1 (void) const { return d1; }
+  octave_idx_type dim2 (void) const { return d2; }
 
   octave_idx_type rows (void) const { return dim1 (); }
   octave_idx_type cols (void) const { return dim2 (); }
   octave_idx_type columns (void) const { return dim2 (); }
 
+  // FIXME: a dangerous ambiguity?
   octave_idx_type length (void) const { return Array<T>::length (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
   octave_idx_type numel (void) const { return nelem (); }
 
   size_t byte_size (void) const { return length () * sizeof (T); }
 
-  dim_vector dims (void) const { return Array<T>::dimensions; }
+  dim_vector dims (void) const { return dim_vector (d1, d2); }
 
   Array<T> diag (octave_idx_type k = 0) const;
 
-  Proxy elem (octave_idx_type r, octave_idx_type c)
-    {
-      return Proxy (this, r, c);
-    }
-
-  Proxy checkelem (octave_idx_type r, octave_idx_type c)
-    {
-      if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
-	{
-	  (*current_liboctave_error_handler) ("range error in DiagArray2");
-	  return Proxy (0, r, c);
-	}
-      else
-	return Proxy (this, r, c);
-    }
-
-  Proxy operator () (octave_idx_type r, octave_idx_type c)
-    {
-      if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
-	{
-	  (*current_liboctave_error_handler) ("range error in DiagArray2");
-	  return Proxy (0, r, c);
-	}
-      else
-	return Proxy (this, r, c);
-  }
+  // Warning: the non-const two-index versions will silently ignore assignments
+  // to off-diagonal elements. 
 
   T elem (octave_idx_type r, octave_idx_type c) const
     {
-      return (r == c) ? Array<T>::xelem (r) : T (0);
+      return (r == c) ? Array<T>::elem (r) : T (0);
+    }
+
+  T& elem (octave_idx_type r, octave_idx_type c)
+    {
+      static T zero (0);
+      return (r == c) ? Array<T>::elem (r) : zero;
     }
 
+  T dgelem (octave_idx_type i) const
+    { return Array<T>::elem (i); }
+
+  T& dgelem (octave_idx_type i) 
+    { return Array<T>::elem (i); }
+
   T checkelem (octave_idx_type r, octave_idx_type c) const;
+  Proxy checkelem (octave_idx_type r, octave_idx_type c);
+
   T operator () (octave_idx_type r, octave_idx_type c) const
     {
 #if defined (BOUNDS_CHECKING)
@@ -209,25 +179,40 @@
 #endif
     }
 
-  // No checking.
-
-  T& xelem (octave_idx_type r, octave_idx_type c)
+  // FIXME: can this cause problems?
+#if defined (BOUNDS_CHECKING)
+  Proxy operator () (octave_idx_type r, octave_idx_type c)
     {
-      static T foo (0);
-      return (r == c) ? Array<T>::xelem (r) : foo;
+      return checkelem (r, c);
     }
+#else
+  T& operator () (octave_idx_type r, octave_idx_type c) 
+    {
+      return elem (r, c);
+    }
+#endif
+
+  // No checking.
 
   T xelem (octave_idx_type r, octave_idx_type c) const
     {
       return (r == c) ? Array<T>::xelem (r) : T (0);
     }
 
+  T& dgxelem (octave_idx_type i)
+    { return Array<T>::xelem (i); }
+
+  T dgxelem (octave_idx_type i) const
+    { return Array<T>::xelem (i); }
+
   void resize (octave_idx_type n, octave_idx_type m);
-  void resize (octave_idx_type n, octave_idx_type m, const T& val);
+  void resize_fill (octave_idx_type n, octave_idx_type m, const T& val);
 
   DiagArray2<T> transpose (void) const;
   DiagArray2<T> hermitian (T (*fcn) (const T&) = 0) const;
 
+  operator Array2<T> (void) const;
+
   const T *data (void) const { return Array<T>::data (); }
 
   const T *fortran_vec (void) const { return Array<T>::fortran_vec (); }

--- a/liboctave/MDiagArray2.h
+++ b/liboctave/MDiagArray2.h
@@ -72,17 +72,7 @@
 
   operator MArray2<T> () const
     {
-      octave_idx_type nr = DiagArray2<T>::dim1 ();
-      octave_idx_type nc = DiagArray2<T>::dim2 ();
-
-      MArray2<T> retval (nr, nc,  T (0));
-
-      octave_idx_type len = nr < nc ? nr : nc;
-
-      for (octave_idx_type i = 0; i < len; i++)
-	retval.xelem (i, i) = this->xelem (i, i);
-
-      return retval;
+      return DiagArray2<T>::operator Array2<T> ();
     }
 
   octave_idx_type nnz (void) const

--- a/liboctave/PermMatrix.cc
+++ b/liboctave/PermMatrix.cc
@@ -39,7 +39,6 @@
 PermMatrix::PermMatrix (const Array<octave_idx_type>& p, bool colp, bool check)
   : Array<octave_idx_type> (p), _colp(colp)
 {
-  this->dimensions = dim_vector (p.length (), p.length ());
   if (check)
     {
       if (! idx_vector (p).is_permutation (p.length ()))
@@ -61,14 +60,12 @@
       Array<octave_idx_type> idxa (len);
       for (octave_idx_type i = 0; i < len; i++) idxa(i) = idx(i);
       Array<octave_idx_type>::operator = (idxa);
-      this->dimensions = dim_vector (len, len);
     }
 }
 
 PermMatrix::PermMatrix (octave_idx_type n)
   : Array<octave_idx_type> (n), _colp (false)
 {
-  this->dimensions = dim_vector (n, n);
   for (octave_idx_type i = 0; i < n; i++) xelem (i) = i;
 }
 

--- a/liboctave/PermMatrix.h
+++ b/liboctave/PermMatrix.h
@@ -26,15 +26,11 @@
 #include "Array.h"
 #include "mx-defs.h"
 
-// Array<T> is inherited privately because we abuse the dimensions variable
-// for true dimensions. Therefore, the inherited Array<T> object is not a valid
-// Array<T> object, and should not be publicly accessible.
+// Array<T> is inherited privately so that some methods, like index, don't
+// produce unexpected results.
 
-class PermMatrix : private Array<octave_idx_type>
+class PermMatrix : protected Array<octave_idx_type>
 {
-private:
-
-  octave_idx_type get (octave_idx_type i) const { return Array<octave_idx_type>::xelem (i); }
 
 public:
 
@@ -51,36 +47,34 @@
   PermMatrix (const idx_vector& idx, bool colp = false, octave_idx_type n = 0); 
 
   octave_idx_type dim1 (void) const 
-    { return Array<octave_idx_type>::dimensions(0); }
+    { return Array<octave_idx_type>::length (); }
   octave_idx_type dim2 (void) const 
-    { return Array<octave_idx_type>::dimensions(1); }
+    { return Array<octave_idx_type>::length (); }
 
   octave_idx_type rows (void) const { return dim1 (); }
   octave_idx_type cols (void) const { return dim2 (); }
   octave_idx_type columns (void) const { return dim2 (); }
 
+  octave_idx_type perm_length (void) const 
+    { return Array<octave_idx_type>::length (); }
   octave_idx_type length (void) const 
-    { return Array<octave_idx_type>::length (); }
+    { return dim1 () * dim2 (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
   octave_idx_type numel (void) const { return nelem (); }
 
-  size_t byte_size (void) const { return length () * sizeof (octave_idx_type); }
+  size_t byte_size (void) const { return perm_length () * sizeof (octave_idx_type); }
 
-  dim_vector dims (void) const { return Array<octave_idx_type>::dimensions; }
+  dim_vector dims (void) const { return dim_vector (dim1 (), dim2 ()); }
 
   Array<octave_idx_type> pvec (void) const
-    {
-      Array<octave_idx_type> retval (*this);
-      retval.dimensions = dim_vector (length ());
-      return retval;
-    }
+    { return *this; }
 
   octave_idx_type 
   elem (octave_idx_type i, octave_idx_type j) const
     {
       return (_colp 
-              ? ((get(j) != i) ? 1 : 0)
-              : ((get(i) != j) ? 1 : 0));
+              ? ((Array<octave_idx_type>::elem (j) != i) ? 1 : 0)
+              : ((Array<octave_idx_type>::elem (i) != j) ? 1 : 0));
     }
 
   octave_idx_type 

--- a/src/strfns.cc
+++ b/src/strfns.cc
@@ -395,7 +395,7 @@
 	    {
 	      // Broadcast the string.
 
-	      boolNDArray output (cell.dimensions, false);
+	      boolNDArray output (cell.dims (), false);
 
 	      std::string s = r == 0 ? std::string () : str[0];
 
@@ -430,7 +430,7 @@
 		{
 		  // Must match in all dimensions.
 
-		  boolNDArray output (cell.dimensions, false);
+		  boolNDArray output (cell.dims (), false);
 
 		  if (cell.length () == r)
 		    {
@@ -471,8 +471,8 @@
 	      r2 = cell2.length ();
 	    }
 
-	  const dim_vector size1 = cell1.dimensions;
-	  const dim_vector size2 = cell2.dimensions;
+	  const dim_vector size1 = cell1.dims ();
+	  const dim_vector size2 = cell2.dims ();
 
 	  boolNDArray output (size1, false);
 
@@ -671,7 +671,7 @@
 	    {
 	      // Broadcast the string.
 
-	      boolNDArray output (cell.dimensions, false);
+	      boolNDArray output (cell.dims (), false);
 
 	      if (c < n)
 		{
@@ -724,7 +724,7 @@
 		{
 		  // Must match in all dimensions.
 
-		  boolNDArray output (cell.dimensions, false);
+		  boolNDArray output (cell.dims (), false);
 
 		  if (cell.numel () == r)
 		    {
@@ -774,8 +774,8 @@
 	      r2 = cell2.length ();
 	    }
 
-	  const dim_vector size1 = cell1.dimensions;
-	  const dim_vector size2 = cell2.dimensions;
+	  const dim_vector size1 = cell1.dims ();
+	  const dim_vector size2 = cell2.dims ();
 
 	  boolNDArray output (size1, false);