octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: goals for 3.1


From: David Bateman
Subject: Re: goals for 3.1
Date: Fri, 21 Mar 2008 00:13:08 +0100
User-agent: Thunderbird 2.0.0.12 (X11/20080306)

John W. Eaton wrote:
> On 20-Mar-2008, David Bateman wrote:
> 
> | David Bateman wrote:
> | 
> | > That is what I just implemented.. Will send a patch after I feed the
> | > kids ...
> | > 
> | > D.
> | >
> | 
> | Changeset attached. I implemented the cell diag function in the same
> | manner as the existing diag function. I will refactor this code later to
> | get rid of the nest of if/else statements..
> 
> I applied it.
> 
> Thanks,
> 
> jwe
> 

Ok and here is the patch that moves the diag function into the
octave_value classes. It also has template versions of the diag function
in the Array<T> and Sparse<T> classes rather than in the derived classes
 to simplify the code.

I believe point 7 on your list for 3.1 is then fully addressed with this
and the previous patches..

D.

# HG changeset patch
# User David Bateman <address@hidden>
# Date 1206054504 -3600
# Node ID ec3993370ae71e19d35b3f0a14788b318dc52b3b
# Parent  489c110b13859e0ef7008b8b28e365a11ffb7bfd
Move diag function into the octave_value class

diff --git a/liboctave/Array.cc b/liboctave/Array.cc
--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -2623,6 +2623,93 @@ Array<double> Array<double>::sort (Array
                                   octave_idx_type dim, sortmode mode) const;
 
 #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.
 
diff --git a/liboctave/Array.h b/liboctave/Array.h
--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -549,6 +549,8 @@ public:
   Array<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const;
   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>
diff --git a/liboctave/Array2.h b/liboctave/Array2.h
--- a/liboctave/Array2.h
+++ b/liboctave/Array2.h
@@ -136,6 +136,11 @@ public:
       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
   {
diff --git a/liboctave/ArrayN.h b/liboctave/ArrayN.h
--- a/liboctave/ArrayN.h
+++ b/liboctave/ArrayN.h
@@ -149,6 +149,11 @@ public:
       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
   {
diff --git a/liboctave/CDiagMatrix.cc b/liboctave/CDiagMatrix.cc
--- a/liboctave/CDiagMatrix.cc
+++ b/liboctave/CDiagMatrix.cc
@@ -499,14 +499,6 @@ operator * (const DiagMatrix& a, const C
 }
 
 // other operations
-
-ComplexColumnVector
-ComplexDiagMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-// Could be optimized...
 
 ComplexColumnVector
 ComplexDiagMatrix::diag (octave_idx_type k) const
diff --git a/liboctave/CDiagMatrix.h b/liboctave/CDiagMatrix.h
--- a/liboctave/CDiagMatrix.h
+++ b/liboctave/CDiagMatrix.h
@@ -114,8 +114,7 @@ public:
 
   // other operations
 
-  ComplexColumnVector diag (void) const;
-  ComplexColumnVector diag (octave_idx_type k) const;
+  ComplexColumnVector diag (octave_idx_type k = 0) const;
 
   // i/o
 
diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -3316,51 +3316,10 @@ Matrix ComplexMatrix::abs (void) const
   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
diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -334,8 +334,7 @@ public:
   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;
diff --git a/liboctave/CNDArray.cc b/liboctave/CNDArray.cc
--- a/liboctave/CNDArray.cc
+++ b/liboctave/CNDArray.cc
@@ -989,6 +989,12 @@ ComplexNDArray::compute_index (Array<oct
                               const dim_vector& dimensions)
 {
   return ::compute_index (ra_idx, dimensions);
+}
+
+ComplexNDArray
+ComplexNDArray::diag (octave_idx_type k) const
+{
+  return MArrayN<Complex>::diag (k);
 }
 
 NDArray
diff --git a/liboctave/CNDArray.h b/liboctave/CNDArray.h
--- a/liboctave/CNDArray.h
+++ b/liboctave/CNDArray.h
@@ -116,6 +116,8 @@ public:
   //  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&);
diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -7359,88 +7359,7 @@ SparseComplexMatrix
 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
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,74 @@ 2008-03-19  John W. Eaton  <address@hidden
+2008-03-20  David Bateman  <address@hidden>
+
+       * 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  <address@hidden>
 
        * oct-env.cc (octave_env::do_base_pathname): Also handle rooted
diff --git a/liboctave/MArray2.h b/liboctave/MArray2.h
--- a/liboctave/MArray2.h
+++ b/liboctave/MArray2.h
@@ -81,6 +81,11 @@ public:
 
   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
   {
diff --git a/liboctave/MArrayN.h b/liboctave/MArrayN.h
--- a/liboctave/MArrayN.h
+++ b/liboctave/MArrayN.h
@@ -98,6 +98,11 @@ public:
 
   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
   {
diff --git a/liboctave/MSparse.h b/liboctave/MSparse.h
--- a/liboctave/MSparse.h
+++ b/liboctave/MSparse.h
@@ -112,7 +112,12 @@ public:
     { 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);
diff --git a/liboctave/Range.cc b/liboctave/Range.cc
--- a/liboctave/Range.cc
+++ b/liboctave/Range.cc
@@ -176,6 +176,40 @@ Range::sort_internal (Array<octave_idx_t
   for (octave_idx_type i = 0; i < nel; i++, tmp += stp)
     psidx[i] = tmp;
 
+}
+
+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
diff --git a/liboctave/Range.h b/liboctave/Range.h
--- a/liboctave/Range.h
+++ b/liboctave/Range.h
@@ -65,6 +65,8 @@ Range
   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,
diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc
--- a/liboctave/Sparse.cc
+++ b/liboctave/Sparse.cc
@@ -2251,6 +2251,156 @@ Sparse<T>::sort (Array<octave_idx_type> 
     }
 
   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
diff --git a/liboctave/Sparse.h b/liboctave/Sparse.h
--- a/liboctave/Sparse.h
+++ b/liboctave/Sparse.h
@@ -522,6 +522,8 @@ public:
   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
diff --git a/liboctave/boolMatrix.cc b/liboctave/boolMatrix.cc
--- a/liboctave/boolMatrix.cc
+++ b/liboctave/boolMatrix.cc
@@ -78,50 +78,9 @@ boolMatrix::operator ! (void) const
 // 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
diff --git a/liboctave/boolMatrix.h b/liboctave/boolMatrix.h
--- a/liboctave/boolMatrix.h
+++ b/liboctave/boolMatrix.h
@@ -64,8 +64,7 @@ public:
 
   // 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;
diff --git a/liboctave/boolNDArray.cc b/liboctave/boolNDArray.cc
--- a/liboctave/boolNDArray.cc
+++ b/liboctave/boolNDArray.cc
@@ -129,6 +129,12 @@ boolNDArray::compute_index (Array<octave
   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, )
 
diff --git a/liboctave/boolNDArray.h b/liboctave/boolNDArray.h
--- a/liboctave/boolNDArray.h
+++ b/liboctave/boolNDArray.h
@@ -109,6 +109,8 @@ public:
       return retval;
     }
 
+  boolNDArray diag (octave_idx_type k = 0) const;
+
 private:
 
   boolNDArray (bool *d, dim_vector& dv) : ArrayN<bool> (d, dv) { }
diff --git a/liboctave/boolSparse.cc b/liboctave/boolSparse.cc
--- a/liboctave/boolSparse.cc
+++ b/liboctave/boolSparse.cc
@@ -143,88 +143,7 @@ SparseBoolMatrix
 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
diff --git a/liboctave/chMatrix.cc b/liboctave/chMatrix.cc
--- a/liboctave/chMatrix.cc
+++ b/liboctave/chMatrix.cc
@@ -192,50 +192,9 @@ charMatrix::extract (octave_idx_type r1,
 }
 
 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
diff --git a/liboctave/chMatrix.h b/liboctave/chMatrix.h
--- a/liboctave/chMatrix.h
+++ b/liboctave/chMatrix.h
@@ -73,8 +73,7 @@ public:
 
   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;
diff --git a/liboctave/chNDArray.cc b/liboctave/chNDArray.cc
--- a/liboctave/chNDArray.cc
+++ b/liboctave/chNDArray.cc
@@ -145,6 +145,12 @@ charNDArray::compute_index (Array<octave
   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
 {
diff --git a/liboctave/chNDArray.h b/liboctave/chNDArray.h
--- a/liboctave/chNDArray.h
+++ b/liboctave/chNDArray.h
@@ -89,6 +89,8 @@ public:
 
   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;
diff --git a/liboctave/dDiagMatrix.cc b/liboctave/dDiagMatrix.cc
--- a/liboctave/dDiagMatrix.cc
+++ b/liboctave/dDiagMatrix.cc
@@ -347,18 +347,13 @@ operator * (const DiagMatrix& a, const D
 // 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)
diff --git a/liboctave/dDiagMatrix.h b/liboctave/dDiagMatrix.h
--- a/liboctave/dDiagMatrix.h
+++ b/liboctave/dDiagMatrix.h
@@ -92,8 +92,7 @@ public:
 
   // other operations
 
-  ColumnVector diag (void) const;
-  ColumnVector diag (octave_idx_type k) const;
+  ColumnVector diag (octave_idx_type k = 0) const;
 
   // i/o
 
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -2843,51 +2843,10 @@ Matrix::abs (void) const
   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
diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -291,8 +291,7 @@ public:
   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;
diff --git a/liboctave/dNDArray.cc b/liboctave/dNDArray.cc
--- a/liboctave/dNDArray.cc
+++ b/liboctave/dNDArray.cc
@@ -965,6 +965,12 @@ NDArray::compute_index (Array<octave_idx
                        const dim_vector& dimensions)
 {
   return ::compute_index (ra_idx, dimensions);
+}
+
+NDArray
+NDArray::diag (octave_idx_type k) const
+{
+  return MArrayN<double>::diag (k);
 }
 
 NDArray
diff --git a/liboctave/dNDArray.h b/liboctave/dNDArray.h
--- a/liboctave/dNDArray.h
+++ b/liboctave/dNDArray.h
@@ -124,6 +124,8 @@ public:
 
   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);
diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -7437,88 +7437,7 @@ SparseMatrix
 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
diff --git a/liboctave/intNDArray.cc b/liboctave/intNDArray.cc
--- a/liboctave/intNDArray.cc
+++ b/liboctave/intNDArray.cc
@@ -60,66 +60,11 @@ intNDArray<T>::any_element_not_one_or_ze
   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.
diff --git a/liboctave/intNDArray.h b/liboctave/intNDArray.h
--- a/liboctave/intNDArray.h
+++ b/liboctave/intNDArray.h
@@ -65,8 +65,7 @@ public:
 
   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.
 
diff --git a/src/Cell.cc b/src/Cell.cc
--- a/src/Cell.cc
+++ b/src/Cell.cc
@@ -239,49 +239,9 @@ Cell::map (ctype_mapper fcn) const
 }
 
 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);
 }
 
 /*
diff --git a/src/Cell.h b/src/Cell.h
--- a/src/Cell.h
+++ b/src/Cell.h
@@ -115,8 +115,7 @@ public:
 
   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); }
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,5 +1,28 @@ 2008-03-20  David Bateman  <address@hidden
 2008-03-20  David Bateman  <address@hidden>
 
+       * 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&,
diff --git a/src/data.cc b/src/data.cc
--- a/src/data.cc
+++ b/src/data.cc
@@ -883,305 +883,6 @@ same orientation as @var{x}.\n\
   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 @@ Given a matrix argument, instead of a ve
   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 ();
 
diff --git a/src/ov-base-mat.h b/src/ov-base-mat.h
--- a/src/ov-base-mat.h
+++ b/src/ov-base-mat.h
@@ -110,6 +110,9 @@ public:
   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,
diff --git a/src/ov-base-scalar.h b/src/ov-base-scalar.h
--- a/src/ov-base-scalar.h
+++ b/src/ov-base-scalar.h
@@ -93,6 +93,9 @@ public:
 
   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,
diff --git a/src/ov-base-sparse.h b/src/ov-base-sparse.h
--- a/src/ov-base-sparse.h
+++ b/src/ov-base-sparse.h
@@ -116,6 +116,9 @@ octave_base_sparse : public octave_base_
   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,
diff --git a/src/ov-base.cc b/src/ov-base.cc
--- a/src/ov-base.cc
+++ b/src/ov-base.cc
@@ -876,6 +876,14 @@ octave_base_value::as_mxArray (void) con
   gripe_wrong_type_arg ("octave_base_value::as_mxArray ()", type_name ());
 
   return 0;
+}
+
+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
diff --git a/src/ov-base.h b/src/ov-base.h
--- a/src/ov-base.h
+++ b/src/ov-base.h
@@ -457,6 +457,8 @@ public:
   virtual octave_idx_type *mex_get_jc (void) const { return 0; }
 
   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;
diff --git a/src/ov-range.h b/src/ov-range.h
--- a/src/ov-range.h
+++ b/src/ov-range.h
@@ -131,6 +131,9 @@ public:
   octave_value all (int dim = 0) const;
 
   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); }
diff --git a/src/ov.h b/src/ov.h
--- a/src/ov.h
+++ b/src/ov.h
@@ -868,6 +868,9 @@ public:
 
   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,

reply via email to

[Prev in Thread] Current Thread [Next in Thread]