changeset 7620:36594d5bbe13

Move diag function into the octave_value class
author David Bateman <dbateman@free.fr>
date Fri, 21 Mar 2008 00:08:24 +0100
parents 56012914972a
children 4682dda22527
files liboctave/Array.cc liboctave/Array.h liboctave/Array2.h liboctave/ArrayN.h liboctave/CDiagMatrix.cc liboctave/CDiagMatrix.h liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/CNDArray.cc liboctave/CNDArray.h liboctave/CSparse.cc liboctave/ChangeLog liboctave/MArray2.h liboctave/MArrayN.h liboctave/MSparse.h liboctave/Range.cc liboctave/Range.h liboctave/Sparse.cc liboctave/Sparse.h liboctave/boolMatrix.cc liboctave/boolMatrix.h liboctave/boolNDArray.cc liboctave/boolNDArray.h liboctave/boolSparse.cc liboctave/chMatrix.cc liboctave/chMatrix.h liboctave/chNDArray.cc liboctave/chNDArray.h liboctave/dDiagMatrix.cc liboctave/dDiagMatrix.h liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/dNDArray.cc liboctave/dNDArray.h liboctave/dSparse.cc liboctave/intNDArray.cc liboctave/intNDArray.h src/Cell.cc src/Cell.h src/ChangeLog src/data.cc src/ov-base-mat.h src/ov-base-scalar.h src/ov-base-sparse.h src/ov-base.cc src/ov-base.h src/ov-range.h src/ov.h
diffstat 48 files changed, 485 insertions(+), 847 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/Array.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Array.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -2624,6 +2624,93 @@
 
 #endif
 
+template <class T>
+Array<T>
+Array<T>::diag (octave_idx_type k) const
+{
+  dim_vector dv = dims ();
+  octave_idx_type nd = dv.length ();
+  Array<T> d;
+
+  if (nd > 2)
+    (*current_liboctave_error_handler) ("Matrix must be 2-dimensional");    
+  else
+    {
+      octave_idx_type nnr = dv (0);
+      octave_idx_type nnc = dv (1);
+
+      if (nnr == 0 || nnc == 0)
+	; // do nothing
+      else if (nnr != 1 && nnc != 1)
+	{
+	  if (k > 0)
+	    nnc -= k;
+	  else if (k < 0)
+	    nnr += k;
+
+	  if (nnr > 0 && nnc > 0)
+	    {
+	      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
+
+	      d.resize (dim_vector (ndiag, 1));
+
+	      if (k > 0)
+		{
+		  for (octave_idx_type i = 0; i < ndiag; i++)
+		    d.xelem (i) = elem (i, i+k);
+		}
+	      else if (k < 0)
+		{
+		  for (octave_idx_type i = 0; i < ndiag; i++)
+		    d.xelem (i) = elem (i-k, i);
+		}
+	      else
+		{
+		  for (octave_idx_type i = 0; i < ndiag; i++)
+		    d.xelem (i) = elem (i, i);
+		}
+	    }
+	  else
+	    (*current_liboctave_error_handler)
+	      ("diag: requested diagonal out of range");
+	}
+      else if (nnr != 0 && nnc != 0)
+	{
+	  octave_idx_type roff = 0;
+	  octave_idx_type coff = 0;
+	  if (k > 0)
+	    {
+	      roff = 0;
+	      coff = k;
+	    }
+	  else if (k < 0)
+	    {
+	      roff = -k;
+	      coff = 0;
+	    }
+
+	  if (nnr == 1)
+	    {
+	      octave_idx_type n = nnc + std::abs (k);
+	      d = Array<T> (dim_vector (n, n), resize_fill_value (T ()));
+
+	      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 ()));
+
+	      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>
--- a/liboctave/Array.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Array.h	Fri Mar 21 00:08:24 2008 +0100
@@ -550,6 +550,8 @@
   Array<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
 		 sortmode mode = ASCENDING) const;
 
+  Array<T> diag (octave_idx_type k = 0) const;
+
   template <class U, class F>
   Array<U>
   map (F fcn) const
--- a/liboctave/Array2.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Array2.h	Fri Mar 21 00:08:24 2008 +0100
@@ -136,6 +136,11 @@
       return Array2<T> (tmp, tmp.rows (), tmp.columns ());
     }
 
+  Array2<T> diag (octave_idx_type k) const
+  {
+    return Array<T>::diag (k);
+  }
+
   template <class U, class F>
   Array2<U> map (F fcn) const
   {
--- a/liboctave/ArrayN.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/ArrayN.h	Fri Mar 21 00:08:24 2008 +0100
@@ -149,6 +149,11 @@
       return ArrayN<T> (tmp, tmp.dims ());
     }
 
+  ArrayN<T> diag (octave_idx_type k) const
+  {
+    return Array<T>::diag (k);
+  }
+
   template <class U, class F>
   ArrayN<U> map (F fcn) const
   {
--- a/liboctave/CDiagMatrix.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CDiagMatrix.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -501,14 +501,6 @@
 // other operations
 
 ComplexColumnVector
-ComplexDiagMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-// Could be optimized...
-
-ComplexColumnVector
 ComplexDiagMatrix::diag (octave_idx_type k) const
 {
   octave_idx_type nnr = rows ();
--- a/liboctave/CDiagMatrix.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CDiagMatrix.h	Fri Mar 21 00:08:24 2008 +0100
@@ -114,8 +114,7 @@
 
   // other operations
 
-  ComplexColumnVector diag (void) const;
-  ComplexColumnVector diag (octave_idx_type k) const;
+  ComplexColumnVector diag (octave_idx_type k = 0) const;
 
   // i/o
 
--- a/liboctave/CMatrix.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CMatrix.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -3316,51 +3316,10 @@
   return retval;
 }
 
-ComplexColumnVector
-ComplexMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-ComplexColumnVector
+ComplexMatrix
 ComplexMatrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  ComplexColumnVector d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      d.resize (ndiag);
-
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i, i+k);
-	}
-      else if (k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i-k, i);
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i, i);
-	}
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return MArray2<Complex>::diag (k);
 }
 
 bool
--- a/liboctave/CMatrix.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CMatrix.h	Fri Mar 21 00:08:24 2008 +0100
@@ -334,8 +334,7 @@
   ComplexMatrix sumsq (int dim = -1) const;
   Matrix abs (void) const;
 
-  ComplexColumnVector diag (void) const;
-  ComplexColumnVector diag (octave_idx_type k) const;
+  ComplexMatrix diag (octave_idx_type k = 0) const;
 
   bool row_is_real_only (octave_idx_type) const;
   bool column_is_real_only (octave_idx_type) const;
--- a/liboctave/CNDArray.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CNDArray.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -991,6 +991,12 @@
   return ::compute_index (ra_idx, dimensions);
 }
 
+ComplexNDArray
+ComplexNDArray::diag (octave_idx_type k) const
+{
+  return MArrayN<Complex>::diag (k);
+}
+
 NDArray
 ComplexNDArray::map (dmapper fcn) const
 {
--- a/liboctave/CNDArray.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CNDArray.h	Fri Mar 21 00:08:24 2008 +0100
@@ -116,6 +116,8 @@
   //  bool all_elements_are_real (void) const;
   //  bool all_integers (double& max_val, double& min_val) const;
 
+  ComplexNDArray diag (octave_idx_type k = 0) const;
+
   typedef double (*dmapper) (const Complex&);
   typedef Complex (*cmapper) (const Complex&);
   typedef bool (*bmapper) (const Complex&);
--- a/liboctave/CSparse.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/CSparse.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -7359,88 +7359,7 @@
 SparseComplexMatrix
 SparseComplexMatrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  SparseComplexMatrix d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      // Count the number of non-zero elements
-      octave_idx_type nel = 0;
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i, i+k) != 0.)
-	      nel++;
-	}
-      else if ( k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i-k, i) != 0.)
-	      nel++;
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i, i) != 0.)
-	      nel++;
-	}
-      
-      d = SparseComplexMatrix (ndiag, 1, nel);
-      d.xcidx (0) = 0;
-      d.xcidx (1) = nel;
-
-      octave_idx_type ii = 0;
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      Complex tmp = elem (i, i+k);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-      else if ( k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      Complex tmp = elem (i-k, i);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      Complex tmp = elem (i, i);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-    }
-  else
-    (*current_liboctave_error_handler) 
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return MSparse<Complex>::diag (k);
 }
 
 SparseMatrix
--- a/liboctave/ChangeLog	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/ChangeLog	Fri Mar 21 00:08:24 2008 +0100
@@ -2,6 +2,77 @@
 
 	* oct-sparse.h: Add headers for amd.h.
 
+2008-03-20  David Bateman  <dbateman@free.fr>
+
+	* Array.cc (Array<T> Array<T>::diag (octave_idx_type) const): New
+	method for diag function.
+	* Array.h  (Array<T> diag (octave_idx_type) const): Declare it.
+	* Array2.h (Array2<T> diag (octave_idx_type) const): New method.
+	* MArray2.h (MArray2<T> diag (octave_idx_type) const): ditto.
+	* ArrayN.h (ArrayN<T> diag (octave_idx_type) const): ditto.
+	* MArrayN.h (MArrayN<T> diag (octave_idx_type) const): ditto.
+
+	* Sparse.cc (Sparse<T> Sparse<T>::diag (octave_idx_type) const):
+	New method for the diag function.
+	* Sparse.h  (Sparse<T> diag (octave_idx_type) const): Declare it.
+	* MSparse.h (MSparse<T> diag (octave_idx_type) const): New method.
+
+	* Range.cc (Matrix Range::diag (octave_idx_type) const):
+	New method for the diag function.
+	* Range.h  (Matrix diag (octave_idx_type) const): Declare it.
+	
+	* CDiagMatrix.cc (ComplexColumnVector ComplexDiagMatrix::diag
+	(void) const): delete.
+	* dDiagMatrix.cc (ColumnVector DiagMatrix::diag (void) const): delete.
+	* dDiagMatrix.h (ColumnVector diag (void) const): ditto.
+	* CMatrix.cc (ComplexColumnVector ComplexMatrix::diag (void) const):
+	delete.
+	* CMatrix.h (ComplexColumnVector diag (void) const): ditto.
+	* dMatrix.cc (ColumnVector Matrix::diag (void) const): ditto.
+	* dMatrix.h (ColumnVector diag (void) const): ditto.
+	* boolMatrix.cc (boolMatrix boolMatrix::diag (void) const): ditto.
+	* boolMatrix.h (boolMatrix diag (void) const): ditto.
+	* chMatrix.cc (charMatrix charMatrix::diag (void) const): ditto.
+	* chMatrix.h (charMatrix diag (void) const): ditto.
+	* intNDArray.cc (intNDArray<T> intNDArray<T>::diag (void) const): ditto.
+	* intNDArray.h (intNDArray<T> diag (void) const): ditto.
+	
+	* CMatrix.cc (ComplexMatrix ComplexMatrix::diag (octave_idx_type)
+	const): Rewrite in terms of template classes function.
+	* CMatrix.h (ComplexMatrix diag (octave_idx_type)const ): Change
+	return type.
+	* dMatrix.cc (Matrix Matrix::diag (octave_idx_type) const): Rewrite in
+	terms of template classes function.
+	* dMatrix.h (Matrix diag (octave_idx_type) const): Change return type.
+	* boolMatrix.cc (boolMatrix boolMatrix::diag (octave_idx_type) const):
+	Rewrite in terms of template classes function.
+	* boolMatrix.h (boolMatrix diag (octave_idx_type) const): Change
+	return type. 
+	* chMatrix.cc (charMatrix charMatrix::diag (octave_idx_type)
+	const): Rewrite in terms of template classes function.
+	
+	* dSparse.cc (SparseMatrix SparseMatrix::diag (octave_idx_type) const):
+	Rewrite in terms of template classes function.
+	* CSparse.cc (SparseComplexMatrix SparseComplexMatrix::diag
+	(octave_idx_type) const): ditto.
+	* boolSparse.cc (SparseBoolMatrix SparseBoolMatrix::diag
+	(octave_idx_type) const): ditto.
+	* intNDArray.cc (intNDArray<T> intNDArray<T>::diag
+	(octave_idx_type) const): ditto.
+
+	* CNDArray.cc (ComplexNDArray ComplexNDArray::diag
+	(octave_idx_type) const): New method.
+	* CNDArray.h (ComplexNDArray diag (octave_idx_type) const):
+	Declare it.
+	* dNDArray.cc (NDArray NDArray::diag (octave_idx_type) const): New
+	method.
+	* dNDArray.h (NDArray diag (octave_idx_type) const): Declare it.
+	* chNDArray.cc (charNDArray charNDArray::diag
+	(octave_idx_type) const): New method.
+	* chNDArray.h (charNDArray diag (octave_idx_type) const):
+	Declare it.
+
+	
 2008-03-19  John W. Eaton  <jwe@octave.org>
 
 	* oct-env.cc (octave_env::do_base_pathname): Also handle rooted
--- a/liboctave/MArray2.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/MArray2.h	Fri Mar 21 00:08:24 2008 +0100
@@ -81,6 +81,11 @@
 
   MArray2<T> transpose (void) const { return Array2<T>::transpose (); }
 
+  MArray2<T> diag (octave_idx_type k) const
+  {
+    return Array2<T>::diag (k);
+  }
+
   template <class U, class F>
   MArray2<U> map (F fcn) const
   {
--- a/liboctave/MArrayN.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/MArrayN.h	Fri Mar 21 00:08:24 2008 +0100
@@ -98,6 +98,11 @@
 
   MArrayN squeeze (void) const { return ArrayN<T>::squeeze (); }
 
+  MArrayN<T> diag (octave_idx_type k) const
+  {
+    return ArrayN<T>::diag (k);
+  }
+
   template <class U, class F>
   MArrayN<U> map (F fcn) const
   {
--- a/liboctave/MSparse.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/MSparse.h	Fri Mar 21 00:08:24 2008 +0100
@@ -112,7 +112,12 @@
     { return Sparse<T>::ipermute (vec); }
 
 
-  template <class U, class F>
+  MSparse<T> diag (octave_idx_type k = 0) const
+  {
+    return Sparse<T>::diag (k);
+  }
+
+ template <class U, class F>
   MSparse<U> map (F fcn) const
   {
     return Sparse<T>::template map<U> (fcn);
--- a/liboctave/Range.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Range.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -178,6 +178,40 @@
 
 }
 
+Matrix 
+Range::diag (octave_idx_type k) const
+{
+  octave_idx_type nnr = 1;
+  octave_idx_type nnc = nelem ();
+  Matrix d;
+
+  if  (nnr != 0)
+    {
+      octave_idx_type roff = 0;
+      octave_idx_type coff = 0;
+      if (k > 0)
+	{
+	  roff = 0;
+	  coff = k;
+	}
+      else if (k < 0)
+	{
+	  roff = -k;
+	  coff = 0;
+	}
+
+      // Force cached matrix to be created
+      matrix_value ();
+
+      octave_idx_type n = nnc + std::abs (k);
+      d = Matrix (n, n, Matrix::resize_fill_value ());
+      for (octave_idx_type i = 0; i < nnc; i++)
+	d.xelem (i+roff, i+coff) = cache.xelem (i);
+    }
+
+  return d;
+}
+
 Range
 Range::sort (octave_idx_type dim, sortmode mode) const
 {
--- a/liboctave/Range.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Range.h	Fri Mar 21 00:08:24 2008 +0100
@@ -65,6 +65,8 @@
   void sort_internal (bool ascending = true);
   void sort_internal (Array<octave_idx_type>& sidx, bool ascending = true);
 
+  Matrix diag (octave_idx_type k = 0) const;
+
   Range sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const;
 
   Range sort (Array<octave_idx_type>& sidx, octave_idx_type dim = 0,
--- a/liboctave/Sparse.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Sparse.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -2253,6 +2253,156 @@
   return m;
 }
 
+template <class T>
+Sparse<T>
+Sparse<T>::diag (octave_idx_type k) const
+{
+  octave_idx_type nnr = rows ();
+  octave_idx_type nnc = cols ();
+  Sparse<T> d;
+
+  if (nnr == 0 || nnc == 0)
+    ; // do nothing
+  else if (nnr != 1 && nnc != 1)
+    {
+      if (k > 0)
+	nnc -= k;
+      else if (k < 0)
+	nnr += k;
+
+      if (nnr > 0 && nnc > 0)
+	{
+	  octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
+
+	  // Count the number of non-zero elements
+	  octave_idx_type nel = 0;
+	  if (k > 0)
+	    {
+	      for (octave_idx_type i = 0; i < ndiag; i++)
+		if (elem (i, i+k) != 0.)
+		  nel++;
+	    }
+	  else if ( k < 0)
+	    {
+	      for (octave_idx_type i = 0; i < ndiag; i++)
+		if (elem (i-k, i) != 0.)
+		  nel++;
+	    }
+	  else
+	    {
+	      for (octave_idx_type i = 0; i < ndiag; i++)
+		if (elem (i, i) != 0.)
+		  nel++;
+	    }
+      
+	  d = Sparse<T> (ndiag, 1, nel);
+	  d.xcidx (0) = 0;
+	  d.xcidx (1) = nel;
+
+	  octave_idx_type ii = 0;
+	  if (k > 0)
+	    {
+	      for (octave_idx_type i = 0; i < ndiag; i++)
+		{
+		  T tmp = elem (i, i+k);
+		  if (tmp != 0.)
+		    {
+		      d.xdata (ii) = tmp;
+		      d.xridx (ii++) = i;
+		    }
+		}
+	    }
+	  else if ( k < 0)
+	    {
+	      for (octave_idx_type i = 0; i < ndiag; i++)
+		{
+		  T tmp = elem (i-k, i);
+		  if (tmp != 0.)
+		    {
+		      d.xdata (ii) = tmp;
+		      d.xridx (ii++) = i;
+		    }
+		}
+	    }
+	  else
+	    {
+	      for (octave_idx_type i = 0; i < ndiag; i++)
+		{
+		  T tmp = elem (i, i);
+		  if (tmp != 0.)
+		    {
+		      d.xdata (ii) = tmp;
+		      d.xridx (ii++) = i;
+		    }
+		}
+	    }
+	}
+      else
+	(*current_liboctave_error_handler) 
+	  ("diag: requested diagonal out of range");
+    }
+  else if (nnr != 0 && nnc != 0)
+    {
+      octave_idx_type roff = 0;
+      octave_idx_type coff = 0;
+      if (k > 0) 
+	{
+	  roff = 0;
+	  coff = k;
+	} 
+      else if (k < 0) 
+	{
+	  roff = -k;
+	  coff = 0;
+	}
+
+      if (nnr == 1) 
+	{
+	  octave_idx_type n = nnc + std::abs (k);
+	  octave_idx_type nz = nzmax ();
+	  d = Sparse<T> (n, n, nz);
+	  for (octave_idx_type i = 0; i < coff+1; i++)
+	    d.xcidx (i) = 0;
+	  for (octave_idx_type j = 0; j < nnc; j++)
+	    {
+	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
+		  d.xdata (i) = data (i);
+		  d.xridx (i) = j + roff;
+		}
+	      d.xcidx (j + coff + 1) = cidx(j+1);
+	    }
+	  for (octave_idx_type i = nnc + coff + 1; i < n + 1; i++)
+	    d.xcidx (i) = nz;
+	} 
+      else 
+	{
+	  octave_idx_type n = nnr + std::abs (k);
+	  octave_idx_type nz = nzmax ();
+	  octave_idx_type ii = 0;
+	  octave_idx_type ir = ridx(0);
+	  d = Sparse<T> (n, n, nz);
+	  for (octave_idx_type i = 0; i < coff+1; i++)
+	    d.xcidx (i) = 0;
+	  for (octave_idx_type i = 0; i < nnr; i++)
+	    {
+	      if (ir == i)
+		{
+		  d.xdata (ii) = data (ii);
+		  d.xridx (ii++) = ir + roff;
+		  if (ii != nz)
+		    ir = ridx (ii);
+		}
+	      d.xcidx (i + coff + 1) = ii;
+	    }
+	  for (octave_idx_type i = nnr + coff + 1; i < n+1; i++)
+	    d.xcidx (i) = nz;
+	}
+    }
+
+  return d;
+}
+
 // FIXME
 // Unfortunately numel can overflow for very large but very sparse matrices.
 // For now just flag an error when this happens.
--- a/liboctave/Sparse.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/Sparse.h	Fri Mar 21 00:08:24 2008 +0100
@@ -522,6 +522,8 @@
   Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
 		 sortmode mode = ASCENDING) const;
 
+  Sparse<T> diag (octave_idx_type k = 0) const;
+
   template <class U, class F>
   Sparse<U>
   map (F fcn) const
--- a/liboctave/boolMatrix.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/boolMatrix.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -78,50 +78,9 @@
 // other operations
 
 boolMatrix
-boolMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-boolMatrix
 boolMatrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  boolMatrix d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      d.resize (ndiag, 1);
-
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.xelem (i) = elem (i, i+k);
-	}
-      else if (k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.xelem (i) = elem (i-k, i);
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.xelem (i) = elem (i, i);
-	}
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return Array2<bool>::diag (k);
 }
 
 // FIXME Do these really belong here?  Maybe they should be
--- a/liboctave/boolMatrix.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/boolMatrix.h	Fri Mar 21 00:08:24 2008 +0100
@@ -64,8 +64,7 @@
 
   // other operations
 
-  boolMatrix diag (void) const;
-  boolMatrix diag (octave_idx_type k) const;
+  boolMatrix diag (octave_idx_type k = 0) const;
 
   boolMatrix all (int dim = -1) const;
   boolMatrix any (int dim = -1) const;
--- a/liboctave/boolNDArray.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/boolNDArray.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -129,6 +129,12 @@
   return ::compute_index (ra_idx, dimensions);
 }
 
+boolNDArray
+boolNDArray::diag (octave_idx_type k) const
+{
+  return ArrayN<bool>::diag (k);
+}
+
 NDND_BOOL_OPS (boolNDArray, boolNDArray, false)
 NDND_CMP_OPS (boolNDArray, , boolNDArray, )
 
--- a/liboctave/boolNDArray.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/boolNDArray.h	Fri Mar 21 00:08:24 2008 +0100
@@ -109,6 +109,8 @@
       return retval;
     }
 
+  boolNDArray diag (octave_idx_type k = 0) const;
+
 private:
 
   boolNDArray (bool *d, dim_vector& dv) : ArrayN<bool> (d, dv) { }
--- a/liboctave/boolSparse.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/boolSparse.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -143,88 +143,7 @@
 SparseBoolMatrix
 SparseBoolMatrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  SparseBoolMatrix d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      // Count the number of non-zero elements
-      octave_idx_type nel = 0;
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i, i+k) != 0.)
-	      nel++;
-	}
-      else if ( k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i-k, i) != 0.)
-	      nel++;
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i, i) != 0.)
-	      nel++;
-	}
-      
-      d = SparseBoolMatrix (ndiag, 1, nel);
-      d.xcidx (0) = 0;
-      d.xcidx (1) = nel;
-
-      octave_idx_type ii = 0;
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      bool tmp = elem (i, i+k);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-      else if ( k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      bool tmp = elem (i-k, i);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      bool tmp = elem (i, i);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-    }
-  else
-    (*current_liboctave_error_handler) 
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return Sparse<bool>::diag (k);
 }
 
 boolMatrix
--- a/liboctave/chMatrix.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/chMatrix.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -192,50 +192,9 @@
 }
 
 charMatrix
-charMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-charMatrix
 charMatrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  charMatrix d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      d.resize (ndiag, 1);
-
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.xelem (i) = elem (i, i+k);
-	}
-      else if (k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.xelem (i) = elem (i-k, i);
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.xelem (i) = elem (i, i);
-	}
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return MArray2<char>::diag (k);
 }
 
 // FIXME Do these really belong here?  Maybe they should be
--- a/liboctave/chMatrix.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/chMatrix.h	Fri Mar 21 00:08:24 2008 +0100
@@ -73,8 +73,7 @@
 
   charMatrix extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const;
 
-  charMatrix diag (void) const;
-  charMatrix diag (octave_idx_type k) const;
+  charMatrix diag (octave_idx_type k = 0) const;
 
   boolMatrix all (int dim = -1) const;
   boolMatrix any (int dim = -1) const;
--- a/liboctave/chNDArray.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/chNDArray.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -145,6 +145,12 @@
   return ::compute_index (ra_idx, dimensions);
 }
 
+charNDArray
+charNDArray::diag (octave_idx_type k) const
+{
+  return MArrayN<char>::diag (k);
+}
+
 boolNDArray
 charNDArray::bmap (mapper fcn) const
 {
--- a/liboctave/chNDArray.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/chNDArray.h	Fri Mar 21 00:08:24 2008 +0100
@@ -89,6 +89,8 @@
 
   static char resize_fill_value (void) { return '\0'; }
 
+  charNDArray diag (octave_idx_type k = 0) const;
+
   typedef int (*mapper) (int);
   boolNDArray bmap (mapper fcn) const;
   NDArray dmap (mapper fcn) const;
--- a/liboctave/dDiagMatrix.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dDiagMatrix.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -347,18 +347,13 @@
 // other operations
 
 ColumnVector
-DiagMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-// Could be optimized...
-
-ColumnVector
 DiagMatrix::diag (octave_idx_type k) const
 {
   octave_idx_type nnr = rows ();
   octave_idx_type nnc = cols ();
+
+  if (nnr == 0  || nnc == 0)
+    
   if (k > 0)
     nnc -= k;
   else if (k < 0)
--- a/liboctave/dDiagMatrix.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dDiagMatrix.h	Fri Mar 21 00:08:24 2008 +0100
@@ -92,8 +92,7 @@
 
   // other operations
 
-  ColumnVector diag (void) const;
-  ColumnVector diag (octave_idx_type k) const;
+  ColumnVector diag (octave_idx_type k = 0) const;
 
   // i/o
 
--- a/liboctave/dMatrix.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dMatrix.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -2843,51 +2843,10 @@
   return retval;
 }
 
-ColumnVector
-Matrix::diag (void) const
-{
-  return diag (0);
-}
-
-ColumnVector
+Matrix
 Matrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  ColumnVector d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      d.resize (ndiag);
-
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i, i+k);
-	}
-      else if (k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i-k, i);
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i, i);
-	}
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return MArray2<double>::diag (k);
 }
 
 ColumnVector
--- a/liboctave/dMatrix.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dMatrix.h	Fri Mar 21 00:08:24 2008 +0100
@@ -291,8 +291,7 @@
   Matrix sumsq (int dim = -1) const;
   Matrix abs (void) const;
 
-  ColumnVector diag (void) const;
-  ColumnVector diag (octave_idx_type k) const;
+  Matrix diag (octave_idx_type k = 0) const;
 
   ColumnVector row_min (void) const;
   ColumnVector row_max (void) const;
--- a/liboctave/dNDArray.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dNDArray.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -968,6 +968,12 @@
 }
 
 NDArray
+NDArray::diag (octave_idx_type k) const
+{
+  return MArrayN<double>::diag (k);
+}
+
+NDArray
 NDArray::map (dmapper fcn) const
 {
   return MArrayN<double>::map<double> (func_ptr (fcn));
--- a/liboctave/dNDArray.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dNDArray.h	Fri Mar 21 00:08:24 2008 +0100
@@ -124,6 +124,8 @@
 
   static double resize_fill_value (void) { return 0; }
 
+  NDArray diag (octave_idx_type k = 0) const;
+
   typedef double (*dmapper) (double);
   typedef Complex (*cmapper) (const Complex&);
   typedef bool (*bmapper) (double);
--- a/liboctave/dSparse.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/dSparse.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -7437,88 +7437,7 @@
 SparseMatrix
 SparseMatrix::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  SparseMatrix d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      // Count the number of non-zero elements
-      octave_idx_type nel = 0;
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i, i+k) != 0.)
-	      nel++;
-	}
-      else if ( k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i-k, i) != 0.)
-	      nel++;
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    if (elem (i, i) != 0.)
-	      nel++;
-	}
-      
-      d = SparseMatrix (ndiag, 1, nel);
-      d.xcidx (0) = 0;
-      d.xcidx (1) = nel;
-
-      octave_idx_type ii = 0;
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      double tmp = elem (i, i+k);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-      else if ( k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      double tmp = elem (i-k, i);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    {
-	      double tmp = elem (i, i);
-	      if (tmp != 0.)
-		{
-		  d.xdata (ii) = tmp;
-		  d.xridx (ii++) = i;
-		}
-	    }
-	}
-    }
-  else
-    (*current_liboctave_error_handler) 
-      ("diag: requested diagonal out of range");
-
-  return d;
+  return MSparse<double>::diag (k);
 }
 
 Matrix
--- a/liboctave/intNDArray.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/intNDArray.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -60,66 +60,11 @@
   return false;
 }
 
-
-template <class T>
-intNDArray<T>
-intNDArray<T>::diag (void) const
-{
-  return diag (0);
-}
-
 template <class T>
 intNDArray<T>
 intNDArray<T>::diag (octave_idx_type k) const
 {
-  dim_vector dv = this->dims ();
-  octave_idx_type nd = dv.length ();
-
-  if (nd > 2)
-    {
-      (*current_liboctave_error_handler) ("Matrix must be 2-dimensional");    
-      return intNDArray<T>();
-    }
-  else
-    {
-      octave_idx_type nnr = dv (0);
-      octave_idx_type nnc = dv (1);
-
-      if (k > 0)
-	nnc -= k;
-      else if (k < 0)
-	nnr += k;
-
-      intNDArray<T> d;
-
-      if (nnr > 0 && nnc > 0)
-	{
-	  octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-	  d.resize (dim_vector (ndiag, 1));
-
-	  if (k > 0)
-	    {
-	      for (octave_idx_type i = 0; i < ndiag; i++)
-		d.xelem (i) = this->elem (i, i+k);
-	    }
-	  else if (k < 0)
-	    {
-	      for (octave_idx_type i = 0; i < ndiag; i++)
-		d.xelem (i) = this->elem (i-k, i);
-	    }
-	  else
-	    {
-	      for (octave_idx_type i = 0; i < ndiag; i++)
-		d.xelem (i) = this->elem (i, i);
-	    }
-	}
-      else
-	(*current_liboctave_error_handler)
-	  ("diag: requested diagonal out of range");
-
-      return d;
-    }
+  return MArrayN<T>::diag (k);
 }
 
 // FIXME -- this is not quite the right thing.
--- a/liboctave/intNDArray.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/liboctave/intNDArray.h	Fri Mar 21 00:08:24 2008 +0100
@@ -65,8 +65,7 @@
 
   bool any_element_not_one_or_zero (void) const;
 
-  intNDArray diag (void) const;
-  intNDArray diag (octave_idx_type k) const;
+  intNDArray diag (octave_idx_type k = 0) const;
 
   // FIXME -- this is not quite the right thing.
 
--- a/src/Cell.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/Cell.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -239,49 +239,9 @@
 }
 
 Cell
-Cell::diag (void) const
-{
-  return diag (0);
-}
-
-Cell
 Cell::diag (octave_idx_type k) const
 {
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  Cell d;
-
-  if (nnr > 0 && nnc > 0)
-    {
-      octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
-
-      d.resize (dim_vector (ndiag, 1));
-
-      if (k > 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i, i+k);
-	}
-      else if (k < 0)
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i-k, i);
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ndiag; i++)
-	    d.elem (i) = elem (i, i);
-	}
-    }
-  else
-    error ("diag: requested diagonal out of range");
-
-  return d;
+  return ArrayN<octave_value>::diag (k);
 }
 
 /*
--- a/src/Cell.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/Cell.h	Fri Mar 21 00:08:24 2008 +0100
@@ -115,8 +115,7 @@
 
   static octave_value resize_fill_value (void) { return Matrix (); }
 
-  Cell diag (void) const;
-  Cell diag (octave_idx_type k) const;
+  Cell diag (octave_idx_type k = 0) const;
 
   Cell xisalnum (void) const { return map (&octave_value::xisalnum); }
   Cell xisalpha (void) const { return map (&octave_value::xisalpha); }
--- a/src/ChangeLog	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ChangeLog	Fri Mar 21 00:08:24 2008 +0100
@@ -5,6 +5,29 @@
 
 2008-03-20  David Bateman  <dbateman@free.fr>
 
+	* Cell.cc (Cell Cell::diag (void) const): delete.
+	(Cell Cell::diag (octave__idx_type) const):Rewrite in terms of
+	template classes function.
+	* Cell.h (Cell diag (void) const): delete.
+
+	* ov.h (octave_value diag (octave_idx_type) const): New method.
+	* ov-base.h (virtual octave_value diag (octave_idx_type) const):
+	New virtual method.
+	* ov-base.cc (octave_value octave_base_value::diag
+	(octave_idx_type) const): New default method.
+	
+	* ov-base-mat.h (octave_value diag (octave_idx_type) const): New
+	method. 
+	* ov-base-sparse.h (octave_value diag (octave_idx_type) const): New
+	method. 
+	* ov-base-scalar.h (octave_value diag (octave_idx_type) const): New
+	method. 
+	* ov-range.h (octave_value diag (octave_idx_type) const): New
+	method. 
+
+	* data.cc (make_diag, make_spdiag): Delete.
+	(Fdiag): Rewrite in terms of octave_value diag function.
+		
 	* data.cc (static octave_value make_diag (const Cell&,
 	octave_idx_type)): New instantiation of template function.
 	(static octave_value make_diag (const octave_value&,
--- a/src/data.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/data.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -883,305 +883,6 @@
   DATA_REDUCTION (cumsum);
 }
 
-template <class T>
-static octave_value
-make_diag (const T& v, octave_idx_type k)
-{
-  octave_value retval;
-  dim_vector dv = v.dims ();
-  octave_idx_type nd = dv.length ();
-
-  if (nd > 2)
-    error ("diag: expecting 2-dimensional matrix");
-  else
-    {
-      octave_idx_type nr = dv (0);
-      octave_idx_type nc = dv (1);
-
-      if (nr == 0 || nc == 0)
-	retval = T ();
-      else if (nr != 1 && nc != 1)
-	retval = v.diag (k);
-      else
-	{
-	  octave_idx_type roff = 0;
-	  octave_idx_type coff = 0;
-	  if (k > 0)
-	    {
-	      roff = 0;
-	      coff = k;
-	    }
-	  else if (k < 0)
-	    {
-	      roff = -k;
-	      coff = 0;
-	    }
-
-	  if (nr == 1)
-	    {
-	      octave_idx_type n = nc + std::abs (k);
-	      T m (dim_vector (n, n), T::resize_fill_value ());
-
-	      for (octave_idx_type i = 0; i < nc; i++)
-		m (i+roff, i+coff) = v (0, i);
-	      retval = m;
-	    }
-	  else
-	    {
-	      octave_idx_type n = nr + std::abs (k);
-	      T m (dim_vector (n, n), T::resize_fill_value ());
-	      for (octave_idx_type i = 0; i < nr; i++)
-		m (i+roff, i+coff) = v (i, 0);
-	      retval = m;
-	    }
-	}
-    }
-  
-  return retval;
-}
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-static octave_value
-make_diag (const Matrix& v, octave_idx_type k);
-
-static octave_value
-make_diag (const ComplexMatrix& v, octave_idx_type k);
-
-static octave_value
-make_diag (const charMatrix& v, octave_idx_type k);
-
-static octave_value
-make_diag (const boolMatrix& v, octave_idx_type k);
-
-static octave_value
-make_diag (const int8NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const int16NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const int32NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const int64NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const uint8NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const uint16NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const uint32NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const uint64NDArray& v, octave_idx_type k);
-
-static octave_value
-make_diag (const Cell& v, octave_idx_type k);
-#endif
-
-template <class T>
-static octave_value
-make_spdiag (const T& v, octave_idx_type k)
-{
-  octave_value retval;
-  dim_vector dv = v.dims ();
-  octave_idx_type nr = dv (0);
-  octave_idx_type nc = dv (1);
-
-  if (nr == 0 || nc == 0)
-    retval = T ();
-  else if (nr != 1 && nc != 1)
-    retval = v.diag (k);
-  else
-    {
-      octave_idx_type roff = 0;
-      octave_idx_type coff = 0;
-      if (k > 0) 
-	{
-	  roff = 0;
-	  coff = k;
-	} 
-      else if (k < 0) 
-	{
-	  roff = -k;
-	  coff = 0;
-	}
-
-      if (nr == 1) 
-	{
-	  octave_idx_type n = nc + std::abs (k);
-	  octave_idx_type nz = v.nzmax ();
-	  T r (n, n, nz);
-	  for (octave_idx_type i = 0; i < coff+1; i++)
-	    r.xcidx (i) = 0;
-	  for (octave_idx_type j = 0; j < nc; j++)
-	    {
-	      for (octave_idx_type i = v.cidx(j); i < v.cidx(j+1); i++)
-		{
-		  r.xdata (i) = v.data (i);
-		  r.xridx (i) = j + roff;
-		}
-	      r.xcidx (j+coff+1) = v.cidx(j+1);
-	    }
-	  for (octave_idx_type i = nc+coff+1; i < n+1; i++)
-	    r.xcidx (i) = nz;
-	  retval = r;
-	} 
-      else 
-	{
-	  octave_idx_type n = nr + std::abs (k);
-	  octave_idx_type nz = v.nzmax ();
-	  octave_idx_type ii = 0;
-	  octave_idx_type ir = v.ridx(0);
-	  T r (n, n, nz);
-	  for (octave_idx_type i = 0; i < coff+1; i++)
-	    r.xcidx (i) = 0;
-	  for (octave_idx_type i = 0; i < nr; i++)
-	    {
-	      if (ir == i)
-		{
-		  r.xdata (ii) = v.data (ii);
-		  r.xridx (ii++) = ir + roff;
-		  if (ii != nz)
-		    ir = v.ridx (ii);
-		}
-	      r.xcidx (i+coff+1) = ii;
-	    }
-	  for (octave_idx_type i = nr+coff+1; i < n+1; i++)
-	    r.xcidx (i) = nz;
-	  retval = r;
-	}
-    }
-
-  return retval;
-}
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-static octave_value
-make_spdiag (const SparseMatrix& v, octave_idx_type k);
-
-static octave_value
-make_spdiag (const SparseComplexMatrix& v, octave_idx_type k);
-
-static octave_value
-make_spdiag (const SparseBoolMatrix& v, octave_idx_type k);
-#endif
-
-static octave_value
-make_diag (const octave_value& a, octave_idx_type k)
-{
-  octave_value retval;
-  std::string result_type = a.class_name ();
-
-  if (result_type == "double")
-    {
-      if (a.is_sparse_type ())
-	{
-	  if (a.is_real_type ())
-	    {
-	      SparseMatrix m = a.sparse_matrix_value ();
-	      if (!error_state)
-		retval = make_spdiag (m, k);
-	    }
-	  else
-	    {
-	      SparseComplexMatrix m = a.sparse_complex_matrix_value ();
-	      if (!error_state)
-		retval = make_spdiag (m, k);
-	    }
-	}
-      else
-	{
-	  if (a.is_real_type ())
-	    {
-	      Matrix m = a.matrix_value ();
-	      if (!error_state)
-		retval = make_diag (m, k);
-	    }
-	  else
-	    {
-	      ComplexMatrix m = a.complex_matrix_value ();
-	      if (!error_state)
-		retval = make_diag (m, k);
-	    }
-	}
-    }
-#if 0
-  else if (result_type == "single")
-    retval = make_diag (a.single_array_value (), k);
-#endif
-  else if (result_type == "char")
-    {
-      charMatrix m = a.char_matrix_value ();
-      if (!error_state)
-	{
-	  retval = make_diag (m, k);
-	  if (a.is_sq_string ())
-	    retval = octave_value (retval.char_array_value (), true, '\'');
-	}
-    }
-  else if (result_type == "logical")
-    {
-      if (a.is_sparse_type ())
-	{
-	  SparseBoolMatrix m = a.sparse_bool_matrix_value ();
-	  if (!error_state)
-	    retval = make_spdiag (m, k);
-	}
-      else
-	{
-	  boolMatrix m = a.bool_matrix_value ();
-	  if (!error_state)
-	    retval = make_diag (m, k);
-	}
-    }
-  else if (result_type == "int8")
-    retval = make_diag (a.int8_array_value (), k);
-  else if (result_type == "int16")
-    retval = make_diag (a.int16_array_value (), k);
-  else if (result_type == "int32")
-    retval = make_diag (a.int32_array_value (), k);
-  else if (result_type == "int64")
-    retval = make_diag (a.int64_array_value (), k);
-  else if (result_type == "uint8")
-    retval = make_diag (a.uint8_array_value (), k);
-  else if (result_type == "uint16")
-    retval = make_diag (a.uint16_array_value (), k);
-  else if (result_type == "uint32")
-    retval = make_diag (a.uint32_array_value (), k);
-  else if (result_type == "uint64")
-    retval = make_diag (a.uint64_array_value (), k);
-  else if (result_type == "cell")
-    retval = make_diag (a.cell_value (), k);
-  else
-    gripe_wrong_type_arg ("diag", a);
-
-  return retval;
-}
-
-static octave_value
-make_diag (const octave_value& arg)
-{
-  return make_diag (arg, 0);
-}
-
-static octave_value
-make_diag (const octave_value& a, const octave_value& b)
-{
-  octave_value retval;
-
-  octave_idx_type k = b.int_value ();
-
-  if (error_state)
-    error ("diag: invalid second argument");      
-  else
-    retval = make_diag (a, k);
-
-  return retval;
-}
-
 DEFUN (diag, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} diag (@var{v}, @var{k})\n\
@@ -1211,9 +912,16 @@
   int nargin = args.length ();
 
   if (nargin == 1 && args(0).is_defined ())
-    retval = make_diag (args(0));
+    retval = args(0).diag();
   else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
-    retval = make_diag (args(0), args(1));
+    {
+      octave_idx_type k = args(1).int_value ();
+
+      if (error_state)
+	error ("diag: invalid second argument");      
+      else
+	retval = args(0).diag(k);
+    }
   else
     print_usage ();
 
--- a/src/ov-base-mat.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov-base-mat.h	Fri Mar 21 00:08:24 2008 +0100
@@ -110,6 +110,9 @@
   MatrixType matrix_type (const MatrixType& _typ) const
     { MatrixType ret = typ; typ = _typ; return ret; }
 
+  octave_value diag (octave_idx_type k = 0) const
+    { return octave_value (matrix.diag (k)); }
+
   octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const
     { return octave_value (matrix.sort (dim, mode)); }
   octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
--- a/src/ov-base-scalar.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov-base-scalar.h	Fri Mar 21 00:08:24 2008 +0100
@@ -93,6 +93,9 @@
 
   octave_value any (int = 0) const { return (scalar != ST ()); }
 
+  octave_value diag (octave_idx_type k = 0) const 
+    { return octave_value (matrix_value (). diag (k)); }
+
   octave_value sort (octave_idx_type, sortmode) const
     { return octave_value (scalar); }
   octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type,
--- a/src/ov-base-sparse.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov-base-sparse.h	Fri Mar 21 00:08:24 2008 +0100
@@ -116,6 +116,9 @@
   octave_value all (int dim = 0) const { return matrix.all (dim); }
   octave_value any (int dim = 0) const { return matrix.any (dim); }
 
+  octave_value diag (octave_idx_type k = 0) const
+    { return octave_value (matrix.diag (k)); }
+
   octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const
     { return octave_value (matrix.sort (dim, mode)); }
   octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
--- a/src/ov-base.cc	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov-base.cc	Fri Mar 21 00:08:24 2008 +0100
@@ -879,6 +879,14 @@
 }
 
 octave_value
+octave_base_value::diag (octave_idx_type) const
+{
+  gripe_wrong_type_arg ("octave_base_value::diag ()", type_name ());
+
+  return octave_value();
+}
+
+octave_value
 octave_base_value::sort (octave_idx_type, sortmode) const
 {
   gripe_wrong_type_arg ("octave_base_value::sort ()", type_name ());
--- a/src/ov-base.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov-base.h	Fri Mar 21 00:08:24 2008 +0100
@@ -458,6 +458,8 @@
 
   virtual mxArray *as_mxArray (void) const;
 
+  virtual octave_value diag (octave_idx_type k = 0) const;
+
   virtual octave_value sort (octave_idx_type dim = 0, 
 			     sortmode mode = ASCENDING) const;
   virtual octave_value sort (Array<octave_idx_type> &sidx, 
--- a/src/ov-range.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov-range.h	Fri Mar 21 00:08:24 2008 +0100
@@ -132,6 +132,9 @@
 
   octave_value any (int dim = 0) const;
 
+  octave_value diag (octave_idx_type k = 0) const
+    { return octave_value (range.diag (k)); }
+
   octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const
     { return range.sort (dim, mode); }
 
--- a/src/ov.h	Fri Mar 21 13:20:11 2008 -0400
+++ b/src/ov.h	Fri Mar 21 00:08:24 2008 +0100
@@ -868,6 +868,9 @@
 
   mxArray *as_mxArray (void) const { return rep->as_mxArray (); }
 
+  octave_value diag (octave_idx_type k = 0) const
+    { return rep->diag (k); }
+
   octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const
     { return rep->sort (dim, mode); } 
   octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,