Mercurial > hg > octave-nkf

changeset 8375:e3c9102431a9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix design problems of diag & perm matrix classes
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 05 Dec 2008 10:20:18 +0100
parents 0f8e810c6950
children c43481a19bfe
files liboctave/CDiagMatrix.cc liboctave/CDiagMatrix.h liboctave/ChangeLog liboctave/CmplxQRP.cc liboctave/DiagArray2.cc liboctave/DiagArray2.h liboctave/MDiagArray2.h liboctave/PermMatrix.cc liboctave/PermMatrix.h liboctave/dDiagMatrix.cc liboctave/dDiagMatrix.h liboctave/dMatrix.cc liboctave/dbleQRP.cc liboctave/fCDiagMatrix.cc liboctave/fCDiagMatrix.h liboctave/fCmplxQRP.cc liboctave/fDiagMatrix.cc liboctave/fDiagMatrix.h liboctave/fMatrix.cc liboctave/floatQRP.cc liboctave/idx-vector.cc liboctave/mx-op-defs.h
diffstat 22 files changed, 155 insertions(+), 204 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CDiagMatrix.cc
+++ b/liboctave/CDiagMatrix.cc
@@ -529,47 +529,6 @@
 
 // other operations
 
-ComplexColumnVector
-ComplexDiagMatrix::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;
-}
-
 ComplexDET
 ComplexDiagMatrix::determinant (void) const
 {
--- a/liboctave/CDiagMatrix.h
+++ b/liboctave/CDiagMatrix.h
@@ -122,7 +122,8 @@
 
   // other operations
 
-  ComplexColumnVector diag (octave_idx_type k = 0) const;
+  ComplexColumnVector diag (octave_idx_type k = 0) const
+    { return MDiagArray2<Complex>::diag (k); }
 
   ComplexDET determinant (void) const;
   double rcond (void) const;
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,21 @@
+2008-12-04  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DiagArray2.h (DiagArray2<T>): Inherit Array<T> privately.
+	(DiagArray2<T>::dim1, dim2, rows, columns, cols, length,
+	nelem, numel, byte_size, dims): New methods.
+	(DiagArray2<T>::diag): New method decl.
+	* DiagArray2.cc (DiagArray2<T>::diag): New method.
+	* MDiagArray2.h (MDiagArray2<T>::diag): New method.
+	* dDiagMatrix.cc (DiagMatrix::diag): Remove.
+	* fDiagMatrix.cc (FloatDiagMatrix::diag): Remove.
+	* CDiagMatrix.cc (ComplexDiagMatrix::diag): Remove.
+	* fCDiagMatrix.cc (FloatComplexDiagMatrix::diag): Remove.
+
+	* PermMatrix.h (PermMatrix): Inherit Array<octave_idx_type> privately.
+	(PermMatrix::dim1, dim2, rows, columns, cols, length,
+	nelem, numel, byte_size, dims): New methods.
+
+	
 2008-12-04  Jaroslav Hajek  <highegg@gmail.com>
 	
 	* dDiagMatrix.cc (DiagMatrix::determinant, DiagMatrix::rcond): New
--- a/liboctave/CmplxQRP.cc
+++ b/liboctave/CmplxQRP.cc
@@ -129,7 +129,7 @@
 ColumnVector
 ComplexQRP::Pvec (void) const
 {
-  Array<double> pa (p);
+  Array<double> pa (p.pvec ());
   ColumnVector pv (MArray<double> (pa) + 1.0);
   return pv;
 }
--- a/liboctave/DiagArray2.cc
+++ b/liboctave/DiagArray2.cc
@@ -30,16 +30,41 @@
 
 #include <iostream>
 
+#include <algorithm>
+
 #include "DiagArray2.h"
 
 #include "lo-error.h"
 
 template <class T>
+Array<T>
+DiagArray2<T>::diag (octave_idx_type k) const
+{
+  Array<T> d;
+
+  if (k == 0)
+    {
+      // The main diagonal is shallow-copied.
+      d = *this;
+      d.dimensions = dim_vector (length ());
+    }
+  else if (k > 0 && k < cols ())
+    d = Array<T> (std::min (cols () - k, rows ()), T ());
+  else if (k < 0 && -k < rows ())
+    d = Array<T> (std::min (rows () + k, cols ()), T ());
+  else
+    (*current_liboctave_error_handler)
+      ("diag: requested diagonal out of range");
+
+  return d;
+}
+
+template <class T>
 DiagArray2<T>
 DiagArray2<T>::transpose (void) const
 {
   DiagArray2<T> retval (*this);
-  retval.dimensions = dim_vector (this->dim2 (), this->dim1 ());
+  retval.dimensions = dim_vector (dim2 (), dim1 ());
   return retval;
 }
 
@@ -47,7 +72,7 @@
 DiagArray2<T>
 DiagArray2<T>::hermitian (T (* fcn) (const T&)) const
 {
-  DiagArray2<T> retval (this->dim2 (), this->dim1 ());
+  DiagArray2<T> retval (dim2 (), dim1 ());
   const T *p = this->data ();
   T *q = retval.fortran_vec ();
   for (octave_idx_type i = 0; i < this->length (); i++)
@@ -61,7 +86,7 @@
 T
 DiagArray2<T>::checkelem (octave_idx_type r, octave_idx_type c) const
 {
-  if (r < 0 || c < 0 || r >= this->dim1 () || c >= this->dim2 ())
+  if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
     {
       (*current_liboctave_error_handler) ("range error in DiagArray2");
       return T ();
@@ -79,7 +104,7 @@
       return;
     }
 
-  if (r == this->dim1 () && c == this->dim2 ())
+  if (r == dim1 () && c == dim2 ())
     return;
 
   typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
@@ -114,7 +139,7 @@
       return;
     }
 
-  if (r == this->dim1 () && c == this->dim2 ())
+  if (r == dim1 () && c == dim2 ())
     return;
 
   typename Array<T>::ArrayRep *old_rep = Array<T>::rep;
--- a/liboctave/DiagArray2.h
+++ b/liboctave/DiagArray2.h
@@ -3,6 +3,7 @@
 
 Copyright (C) 1996, 1997, 2000, 2002, 2003, 2004, 2005, 2006, 2007
               John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -31,8 +32,6 @@
 #include "Array.h"
 #include "lo-error.h"
 
-class idx_vector;
-
 // A two-dimensional array with diagonal elements only.
 //
 // Idea and example code for Proxy class and functions from:
@@ -46,9 +45,13 @@
 // James Kanze                             email: kanze@us-es.sel.de
 // GABI Software, Sarl., 8 rue du Faisan, F-67000 Strasbourg, France
 
+// Array<T> is inherited privately because we abuse the dimensions variable
+// for true dimensions. Therefore, the inherited Array<T> object is not a valid
+// Array<T> object, and should not be publicly accessible.
+
 template <class T>
 class
-DiagArray2 : public Array<T>
+DiagArray2 : private Array<T>
 {
 private:
 
@@ -112,6 +115,8 @@
 
 public:
 
+  typedef T element_type;
+
   DiagArray2 (void) : Array<T> (dim_vector (0, 0)) { }
 
   DiagArray2 (octave_idx_type r, octave_idx_type c) : Array<T> (r < c ? r : c)
@@ -131,7 +136,7 @@
     { this->dimensions = a.dims (); }
 
   template <class U>
-  DiagArray2 (const DiagArray2<U>& a) : Array<T> (a)
+  DiagArray2 (const DiagArray2<U>& a) : Array<T> (a.diag ())
     { this->dimensions = a.dims (); }
 
   ~DiagArray2 (void) { }
@@ -144,6 +149,24 @@
       return *this;
     }
 
+
+  octave_idx_type dim1 (void) const { return Array<T>::dimensions(0); }
+  octave_idx_type dim2 (void) const { return Array<T>::dimensions(1); }
+
+  octave_idx_type rows (void) const { return dim1 (); }
+  octave_idx_type cols (void) const { return dim2 (); }
+  octave_idx_type columns (void) const { return dim2 (); }
+
+  octave_idx_type length (void) const { return Array<T>::length (); }
+  octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
+  octave_idx_type numel (void) const { return nelem (); }
+
+  size_t byte_size (void) const { return length () * sizeof (T); }
+
+  dim_vector dims (void) const { return Array<T>::dimensions; }
+
+  Array<T> diag (octave_idx_type k = 0) const;
+
   Proxy elem (octave_idx_type r, octave_idx_type c)
     {
       return Proxy (this, r, c);
@@ -151,7 +174,7 @@
 
   Proxy checkelem (octave_idx_type r, octave_idx_type c)
     {
-      if (r < 0 || c < 0 || r >= this->dim1 () || c >= this->dim2 ())
+      if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
 	{
 	  (*current_liboctave_error_handler) ("range error in DiagArray2");
 	  return Proxy (0, r, c);
@@ -162,7 +185,7 @@
 
   Proxy operator () (octave_idx_type r, octave_idx_type c)
     {
-      if (r < 0 || c < 0 || r >= this->dim1 () || c >= this->dim2 ())
+      if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
 	{
 	  (*current_liboctave_error_handler) ("range error in DiagArray2");
 	  return Proxy (0, r, c);
@@ -204,6 +227,15 @@
 
   DiagArray2<T> transpose (void) const;
   DiagArray2<T> hermitian (T (*fcn) (const T&) = 0) const;
+
+  const T *data (void) const { return Array<T>::data (); }
+
+  const T *fortran_vec (void) const { return Array<T>::fortran_vec (); }
+
+  T *fortran_vec (void) { return Array<T>::fortran_vec (); }
+
+  void print_info (std::ostream& os, const std::string& prefix) const
+    { Array<T>::print_info (os, prefix); }
 };
 
 #endif
--- a/liboctave/MDiagArray2.h
+++ b/liboctave/MDiagArray2.h
@@ -27,6 +27,7 @@
 
 #include "DiagArray2.h"
 #include "MArray2.h"
+#include "MArray.h"
 
 // Two dimensional diagonal array with math ops.
 
@@ -88,9 +89,9 @@
     {
       octave_idx_type retval = 0;
 
-      const T *d = this->Array<T>::data ();
+      const T *d = this->data ();
 
-      octave_idx_type nel = this->Array<T>::numel ();
+      octave_idx_type nel = this->length ();
 
       for (octave_idx_type i = 0; i < nel; i++)
 	{
@@ -101,6 +102,9 @@
       return retval;
     }
 
+  MArray<T> diag (octave_idx_type k = 0) const
+    { return DiagArray2<T>::diag (k); }
+
   MDiagArray2<T> transpose (void) const { return DiagArray2<T>::transpose (); }
   MDiagArray2<T> hermitian (T (*fcn) (const T&) = 0) const { return DiagArray2<T>::hermitian (fcn); }
 
--- a/liboctave/PermMatrix.cc
+++ b/liboctave/PermMatrix.cc
@@ -123,7 +123,7 @@
 PermMatrix 
 operator *(const PermMatrix& a, const PermMatrix& b)
 {
-  const Array<octave_idx_type>& ia = a, ib = b;
+  const Array<octave_idx_type> ia = a.pvec (), ib = b.pvec ();
   PermMatrix r;
   octave_idx_type n = a.columns ();
   if (n != b.rows ())
--- a/liboctave/PermMatrix.h
+++ b/liboctave/PermMatrix.h
@@ -26,7 +26,11 @@
 #include "Array.h"
 #include "mx-defs.h"
 
-class PermMatrix : public Array<octave_idx_type>
+// Array<T> is inherited privately because we abuse the dimensions variable
+// for true dimensions. Therefore, the inherited Array<T> object is not a valid
+// Array<T> object, and should not be publicly accessible.
+
+class PermMatrix : private Array<octave_idx_type>
 {
 private:
 
@@ -42,11 +46,35 @@
               bool check = true);
 
   PermMatrix (const PermMatrix& m)
-    : Array<octave_idx_type> (m), _colp(m._colp) 
-    { this->dimensions = m.dims (); }
+    : Array<octave_idx_type> (m), _colp(m._colp) { }
   
   PermMatrix (const idx_vector& idx, bool colp = false, octave_idx_type n = 0); 
 
+  octave_idx_type dim1 (void) const 
+    { return Array<octave_idx_type>::dimensions(0); }
+  octave_idx_type dim2 (void) const 
+    { return Array<octave_idx_type>::dimensions(1); }
+
+  octave_idx_type rows (void) const { return dim1 (); }
+  octave_idx_type cols (void) const { return dim2 (); }
+  octave_idx_type columns (void) const { return dim2 (); }
+
+  octave_idx_type length (void) const 
+    { return Array<octave_idx_type>::length (); }
+  octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
+  octave_idx_type numel (void) const { return nelem (); }
+
+  size_t byte_size (void) const { return length () * sizeof (octave_idx_type); }
+
+  dim_vector dims (void) const { return Array<octave_idx_type>::dimensions; }
+
+  Array<octave_idx_type> pvec (void) const
+    {
+      Array<octave_idx_type> retval (*this);
+      retval.dimensions = dim_vector (length ());
+      return retval;
+    }
+
   octave_idx_type 
   elem (octave_idx_type i, octave_idx_type j) const
     {
@@ -80,6 +108,18 @@
 
   friend PermMatrix operator *(const PermMatrix& a, const PermMatrix& b);
 
+  const octave_idx_type *data (void) const 
+    { return Array<octave_idx_type>::data (); }
+
+  const octave_idx_type *fortran_vec (void) const 
+    { return Array<octave_idx_type>::fortran_vec (); }
+
+  octave_idx_type *fortran_vec (void) 
+    { return Array<octave_idx_type>::fortran_vec (); }
+
+  void print_info (std::ostream& os, const std::string& prefix) const
+    { Array<octave_idx_type>::print_info (os, prefix); }
+
 private:
   bool _colp;
 };
--- a/liboctave/dDiagMatrix.cc
+++ b/liboctave/dDiagMatrix.cc
@@ -342,51 +342,6 @@
 
 // other operations
 
-ColumnVector
-DiagMatrix::diag (octave_idx_type k) const
-{
-  ColumnVector d;
-
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-
-  if (nnr == 0 || nnc == 0)
-    return d;
-    
-  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 (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;
-}
-
 DET
 DiagMatrix::determinant (void) const
 {
--- a/liboctave/dDiagMatrix.h
+++ b/liboctave/dDiagMatrix.h
@@ -97,7 +97,8 @@
 
   // other operations
 
-  ColumnVector diag (octave_idx_type k = 0) const;
+  ColumnVector diag (octave_idx_type k = 0) const
+    { return MDiagArray2<double>::diag (k); }
 
   DET determinant (void) const;
   double rcond (void) const;
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -251,7 +251,7 @@
 Matrix::Matrix (const PermMatrix& a)
   : MArray2<double> (a.rows (), a.cols (), 0.0)
 {
-  const Array<octave_idx_type> ia (a);
+  const Array<octave_idx_type> ia (a.pvec ());
   octave_idx_type len = a.rows ();
   if (a.is_col_perm ())
     for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/dbleQRP.cc
+++ b/liboctave/dbleQRP.cc
@@ -123,7 +123,7 @@
 ColumnVector
 QRP::Pvec (void) const
 {
-  Array<double> pa (p);
+  Array<double> pa (p.pvec ());
   ColumnVector pv (MArray<double> (pa) + 1.0);
   return pv;
 }
--- a/liboctave/fCDiagMatrix.cc
+++ b/liboctave/fCDiagMatrix.cc
@@ -550,47 +550,6 @@
 
 // other operations
 
-FloatComplexColumnVector
-FloatComplexDiagMatrix::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;
-
-  FloatComplexColumnVector 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;
-}
-
 FloatComplexDET
 FloatComplexDiagMatrix::determinant (void) const
 {
--- a/liboctave/fCDiagMatrix.h
+++ b/liboctave/fCDiagMatrix.h
@@ -122,7 +122,8 @@
 
   // other operations
 
-  FloatComplexColumnVector diag (octave_idx_type k = 0) const;
+  FloatComplexColumnVector diag (octave_idx_type k = 0) const
+    { return MDiagArray2<FloatComplex>::diag (k); }
 
   FloatComplexDET determinant (void) const;
   float rcond (void) const;
--- a/liboctave/fCmplxQRP.cc
+++ b/liboctave/fCmplxQRP.cc
@@ -129,7 +129,7 @@
 FloatColumnVector
 FloatComplexQRP::Pvec (void) const
 {
-  Array<float> pa (p);
+  Array<float> pa (p.pvec ());
   FloatColumnVector pv (MArray<float> (pa) + 1.0f);
   return pv;
 }
--- a/liboctave/fDiagMatrix.cc
+++ b/liboctave/fDiagMatrix.cc
@@ -349,51 +349,6 @@
 
 // other operations
 
-FloatColumnVector
-FloatDiagMatrix::diag (octave_idx_type k) const
-{
-  FloatColumnVector d;
-
-  octave_idx_type nnr = rows ();
-  octave_idx_type nnc = cols ();
-
-  if (nnr == 0  || nnc == 0)
-    return d;
-    
-  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 (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;
-}
-
 FloatDET
 FloatDiagMatrix::determinant (void) const
 {
--- a/liboctave/fDiagMatrix.h
+++ b/liboctave/fDiagMatrix.h
@@ -97,7 +97,8 @@
 
   // other operations
 
-  FloatColumnVector diag (octave_idx_type k = 0) const;
+  FloatColumnVector diag (octave_idx_type k = 0) const
+    { return MDiagArray2<float>::diag (k); }
 
   FloatDET determinant (void) const;
   float rcond (void) const;
--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -250,7 +250,7 @@
 FloatMatrix::FloatMatrix (const PermMatrix& a)
   : MArray2<float> (a.rows (), a.cols (), 0.0)
 {
-  const Array<octave_idx_type> ia (a);
+  const Array<octave_idx_type> ia (a.pvec ());
   octave_idx_type len = a.rows ();
   if (a.is_col_perm ())
     for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/floatQRP.cc
+++ b/liboctave/floatQRP.cc
@@ -123,7 +123,7 @@
 FloatColumnVector
 FloatQRP::Pvec (void) const
 {
-  Array<float> pa (p);
+  Array<float> pa (p.pvec ());
   FloatColumnVector pv (MArray<float> (pa) + 1.0f);
   return pv;
 }
--- a/liboctave/idx-vector.cc
+++ b/liboctave/idx-vector.cc
@@ -557,7 +557,7 @@
             }
         }
 
-      delete left;
+      delete [] left;
     }
 
   return retval;
--- a/liboctave/mx-op-defs.h
+++ b/liboctave/mx-op-defs.h
@@ -1301,10 +1301,10 @@
       if (p.is_col_perm ()) \
         { \
           result = M (nr, nc); \
-          result.assign (idx_vector (p), idx_vector::colon, x); \
+          result.assign (p.pvec (), idx_vector::colon, x); \
         } \
       else \
-        result = x.index (idx_vector (p), idx_vector::colon); \
+        result = x.index (p.pvec (), idx_vector::colon); \
     } \
   \
   return result; \
@@ -1320,11 +1320,11 @@
   else \
     { \
       if (p.is_col_perm ()) \
-        result = x.index (idx_vector::colon, idx_vector (p)); \
+        result = x.index (idx_vector::colon, p.pvec ()); \
       else \
         { \
           result = M (nr, nc); \
-          result.assign (idx_vector::colon, idx_vector (p), x); \
+          result.assign (idx_vector::colon, p.pvec (), x); \
         } \
     } \
   \