Mercurial > hg > octave-lyh

changeset 10363:a0728e81ed25

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve diag matrix interface & implementation
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 26 Feb 2010 11:44:38 +0100
parents b47ab50a6aa8
children 96ed7c629bbd
files liboctave/CColVector.cc liboctave/CColVector.h liboctave/CDiagMatrix.cc liboctave/CDiagMatrix.h liboctave/CMatrix.cc liboctave/CRowVector.cc liboctave/ChangeLog liboctave/DiagArray2.cc liboctave/DiagArray2.h liboctave/MDiagArray2.h liboctave/dColVector.cc liboctave/dColVector.h liboctave/dDiagMatrix.cc liboctave/dDiagMatrix.h liboctave/dMatrix.cc liboctave/dRowVector.cc liboctave/fCColVector.cc liboctave/fCColVector.h liboctave/fCDiagMatrix.cc liboctave/fCDiagMatrix.h liboctave/fCMatrix.cc liboctave/fCRowVector.cc liboctave/fColVector.cc liboctave/fColVector.h liboctave/fDiagMatrix.cc liboctave/fDiagMatrix.h liboctave/fMatrix.cc liboctave/fRowVector.cc src/ChangeLog src/ov-base-diag.cc src/xpow.cc
diffstat 31 files changed, 146 insertions(+), 287 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CColVector.cc
+++ b/liboctave/CColVector.cc
@@ -231,14 +231,16 @@
   return MArray<Complex>::transpose ();
 }
 
+ColumnVector
+ComplexColumnVector::abs (void) const
+{
+  return do_mx_unary_map<double, Complex, std::abs> (*this);
+}
+
 ComplexColumnVector
 conj (const ComplexColumnVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  ComplexColumnVector retval;
-  if (a_len > 0)
-    retval = ComplexColumnVector (mx_inline_conj_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_map<Complex, Complex, std::conj> (a);
 }
 
 // resize is the destructive equivalent for this one
--- a/liboctave/CColVector.h
+++ b/liboctave/CColVector.h
@@ -121,6 +121,8 @@
   Complex min (void) const;
   Complex max (void) const;
 
+  ColumnVector abs (void) const;
+
   // i/o
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ComplexColumnVector& a);
--- a/liboctave/CDiagMatrix.cc
+++ b/liboctave/CDiagMatrix.cc
@@ -237,21 +237,13 @@
 DiagMatrix
 ComplexDiagMatrix::abs (void) const
 {
-  DiagMatrix retval (rows (), cols ());
-  for (octave_idx_type i = 0; i < rows (); i++)
-    retval(i, i) = std::abs (elem (i, i));
-  return retval;
+  return DiagMatrix (diag ().abs (), rows (), columns ());
 }
 
 ComplexDiagMatrix
 conj (const ComplexDiagMatrix& a)
 {
-  ComplexDiagMatrix retval;
-  octave_idx_type a_len = a.length ();
-  if (a_len > 0)
-    retval = ComplexDiagMatrix (mx_inline_conj_dup (a.data (), a_len),
-                                a.rows (), a.cols ());
-  return retval;
+  return ComplexDiagMatrix (conj (a.diag ()), a.rows (), a.columns ());
 }
 
 // resize is the destructive analog for this one
--- a/liboctave/CDiagMatrix.h
+++ b/liboctave/CDiagMatrix.h
@@ -47,17 +47,14 @@
   ComplexDiagMatrix (octave_idx_type r, octave_idx_type c, const Complex& val)
     : MDiagArray2<Complex> (r, c, val) { }
 
-  explicit ComplexDiagMatrix (const RowVector& a)
-    : MDiagArray2<Complex> (ComplexRowVector (a)) { }
-
-  explicit ComplexDiagMatrix (const ComplexRowVector& a)
+  explicit ComplexDiagMatrix (const Array<Complex>& a)
     : MDiagArray2<Complex> (a) { }
 
-  explicit ComplexDiagMatrix (const ColumnVector& a)
-    : MDiagArray2<Complex> (ComplexColumnVector (a)) { }
+  explicit ComplexDiagMatrix (const Array<double>& a)
+    : MDiagArray2<Complex> (Array<Complex> (a)) { }
 
-  explicit ComplexDiagMatrix (const ComplexColumnVector& a)
-    : MDiagArray2<Complex> (a) { }
+  ComplexDiagMatrix (const Array<Complex>& a, octave_idx_type r, octave_idx_type c) 
+    : MDiagArray2<Complex> (a, r, c) { }
 
   explicit ComplexDiagMatrix (const DiagMatrix& a);
 
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -892,8 +892,7 @@
 ComplexMatrix
 conj (const ComplexMatrix& a)
 {
-  return ComplexMatrix (mx_inline_conj_dup (a.data (), a.length ()),
-                        a.rows (), a.cols ());
+  return do_mx_unary_map<Complex, Complex, std::conj> (a);
 }
 
 // resize is the destructive equivalent for this one
@@ -3243,8 +3242,7 @@
 
 Matrix ComplexMatrix::abs (void) const
 {
-  return Matrix (mx_inline_cabs_dup (data (), length ()),
-                 rows (), cols ());
+  return do_mx_unary_map<double, Complex, std::abs> (*this);
 }
 
 ComplexMatrix
--- a/liboctave/CRowVector.cc
+++ b/liboctave/CRowVector.cc
@@ -232,11 +232,7 @@
 ComplexRowVector
 conj (const ComplexRowVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  ComplexRowVector retval;
-  if (a_len > 0)
-    retval = ComplexRowVector (mx_inline_conj_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_map<Complex, Complex, std::conj> (a);
 }
 
 // resize is the destructive equivalent for this one
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,43 @@
+2010-02-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DiagArray2.h (DiagArray2): Clean up ctors. Remove Proxy class.
+	Remove non-const two-index elem methods.
+	* DiagArray2.cc: Update.
+	* MDiagArray2.h: Update.
+
+	* dDiagArray2.h: Update.
+	* dDiagArray2.cc: Update.
+	* fDiagArray2.h: Update.
+	* fDiagArray2.cc: Update.
+	* CDiagArray2.h: Update.
+	* CDiagArray2.cc: Update.
+	* fCDiagArray2.h: Update.
+	* fCDiagArray2.cc: Update.
+
+	* dMatrix.cc: Update.
+	* fMatrix.cc: Update.
+	* CMatrix.cc: Update.
+	* fCMatrix.cc: Update.
+
+	* dRowVector.cc: Update.
+	* fRowVector.cc: Update.
+	* CRowVector.cc: Update.
+	* fCRowVector.cc: Update.
+
+	* dColVector.cc: Update.
+	* fColVector.cc: Update.
+	* CColVector.cc: Update.
+	* fCColVector.cc: Update.
+
+	* dColVector.cc (ColumnVector::abs): New method.
+	* dColVector.h: Declare it.
+	* fColVector.cc (FloatColumnVector::abs): New method.
+	* fColVector.h: Declare it.
+	* CColVector.cc (ComplexColumnVector::abs): New method.
+	* CColVector.h: Declare it.
+	* fCColVector.cc (FloatComplexColumnVector::abs): New method.
+	* fCColVector.h: Declare it.
+
 2010-02-26  Jaroslav Hajek  <highegg@gmail.com>
 
 	* mx-inlines.cc: Parameterize all appliers by value types rather than
--- a/liboctave/DiagArray2.cc
+++ b/liboctave/DiagArray2.cc
@@ -37,31 +37,12 @@
 #include "lo-error.h"
 
 template <class T>
-const typename DiagArray2<T>::Proxy& 
-DiagArray2<T>::Proxy::operator = (const T& val) const
+DiagArray2<T>::DiagArray2 (const Array<T>& a, octave_idx_type r, octave_idx_type c)
+  : Array<T> (a.as_column ()), d1 (r), d2 (c)
 {
-  if (i == j)
-    {
-      if (object)
-        object->set (val, i);
-    }
-  else
-    (*current_liboctave_error_handler)
-      ("invalid assignment to off-diagonal in diagonal array");
-
-  return *this;
-}
-
-template <class T>
-DiagArray2<T>::Proxy::operator T () const
-{
-  if (object && i == j)
-    return object->get (i);
-  else
-    {
-      static T foo;
-      return foo;
-    }
+  octave_idx_type rcmin = std::min (r, c);
+  if (rcmin != a.length ())
+      Array<T>::resize (rcmin, 1);
 }
 
 template <class T>
@@ -88,22 +69,14 @@
 DiagArray2<T>
 DiagArray2<T>::transpose (void) const
 {
-  DiagArray2<T> retval (*this);
-  retval.d1 = d2;
-  retval.d2 = d1;
-  return retval;
+  return DiagArray2<T> (*this, d2, d1);
 }
 
 template <class T>
 DiagArray2<T>
 DiagArray2<T>::hermitian (T (* fcn) (const T&)) const
 {
-  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++)
-    q [i] = fcn (p [i]);
-  return retval;
+  return DiagArray2<T> (Array<T>::template map<T> (fcn), d2, d1);
 }
 
 // A two-dimensional array with diagonal elements only.
@@ -113,24 +86,9 @@
 DiagArray2<T>::checkelem (octave_idx_type r, octave_idx_type c) const
 {
   if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
-    {
-      (*current_liboctave_error_handler) ("range error in DiagArray2");
-      return T ();
-    }
-  return elem (r, c);
-}
+    (*current_liboctave_error_handler) ("range error in DiagArray2");
 
-template <class T>
-typename DiagArray2<T>::Proxy
-DiagArray2<T>::checkelem (octave_idx_type r, octave_idx_type c) 
-{
-  if (r < 0 || c < 0 || r >= dim1 () || c >= dim2 ())
-    {
-      (*current_liboctave_error_handler) ("range error in DiagArray2");
-      return Proxy (0, r, c);
-    }
-  else
-    return Proxy (this, r, c);
+  return elem (r, c);
 }
 
 template <class T>
@@ -152,7 +110,7 @@
 }
 
 template <class T>
-DiagArray2<T>::operator Array<T> (void) const
+Array<T> DiagArray2<T>::array_value (void) const
 {
   Array<T> result (dim1 (), dim2 ());
   for (octave_idx_type i = 0, len = length (); i < len; i++)
--- a/liboctave/DiagArray2.h
+++ b/liboctave/DiagArray2.h
@@ -30,19 +30,6 @@
 #include <cstdlib>
 
 #include "Array.h"
-#include "lo-error.h"
-
-// A two-dimensional array with diagonal elements only.
-// Idea and example code for Proxy class and functions from:
-//
-// From: kanze@us-es.sel.de (James Kanze)
-// Subject: Re: How to overload [] to do READ/WRITE differently ?
-// Message-ID: <KANZE.93Nov29151407@slsvhdt.us-es.sel.de>
-// Sender: news@us-es.sel.de
-// Date: 29 Nov 1993 14:14:07 GMT
-// --
-// James Kanze                             email: kanze@us-es.sel.de
-// GABI Software, Sarl., 8 rue du Faisan, F-67000 Strasbourg, France
 
 // Array<T> is inherited privately so that some methods, like index, don't
 // produce unexpected results.
@@ -51,40 +38,6 @@
 class
 DiagArray2 : protected Array<T>
 {
-private:
-
-  T get (octave_idx_type i) { return Array<T>::xelem (i); }
-
-  void set (const T& val, octave_idx_type i) { Array<T>::xelem (i) = val; }
-
-  class Proxy
-  {
-  public:
-
-    Proxy (DiagArray2<T> *ref, octave_idx_type r, octave_idx_type c)
-      : i (r), j (c), object (ref) { } 
-
-    const Proxy& operator = (const T& val) const;
-
-    operator T () const;
-
-  private:
-
-    // FIXME -- this is declared private to keep the user from
-    // taking the address of a Proxy.  Maybe it should be implemented
-    // by means of a companion function in the DiagArray2 class.
-
-    T *operator& () const { assert (0); return 0; }
-
-    octave_idx_type i;
-    octave_idx_type j;
-
-    DiagArray2<T> *object;
-
-  };
-
-  friend class Proxy;
-
 protected:
   octave_idx_type d1, d2;
 
@@ -104,18 +57,10 @@
   DiagArray2 (octave_idx_type r, octave_idx_type c, const T& val) 
     : Array<T> (std::min (r, c), 1, val), d1 (r), d2 (c) { }
 
-  DiagArray2 (const dim_vector& dv)
-    : Array<T> (std::min (dv(0), dv(1)), 1), d1 (dv(0)), d2 (dv(0))
-    {
-      if (dv.length () != 2)
-        (*current_liboctave_error_handler) ("too many dimensions");
-    }
-
-  DiagArray2 (const Array<T>& a) 
+  explicit DiagArray2 (const Array<T>& a) 
     : Array<T> (a.as_column ()), d1 (a.numel ()), d2 (a.numel ()) { }
 
-  DiagArray2 (const Array<T>& a, octave_idx_type r, octave_idx_type c) 
-    : Array<T> (a.as_column ()), d1 (r), d2 (c) { }
+  DiagArray2 (const Array<T>& a, octave_idx_type r, octave_idx_type c);
 
   DiagArray2 (const DiagArray2<T>& a) 
     : Array<T> (a), d1 (a.d1), d2 (a.d2) { }
@@ -178,7 +123,6 @@
     { return Array<T>::elem (i); }
 
   T checkelem (octave_idx_type r, octave_idx_type c) const;
-  Proxy checkelem (octave_idx_type r, octave_idx_type c);
 
   T operator () (octave_idx_type r, octave_idx_type c) const
     {
@@ -189,19 +133,6 @@
 #endif
     }
 
-  // FIXME: can this cause problems?
-#if defined (BOUNDS_CHECKING)
-  Proxy operator () (octave_idx_type r, octave_idx_type c)
-    {
-      return checkelem (r, c);
-    }
-#else
-  T& operator () (octave_idx_type r, octave_idx_type c) 
-    {
-      return elem (r, c);
-    }
-#endif
-
   // No checking.
 
   T xelem (octave_idx_type r, octave_idx_type c) const
@@ -221,7 +152,7 @@
   DiagArray2<T> transpose (void) const;
   DiagArray2<T> hermitian (T (*fcn) (const T&) = 0) const;
 
-  operator Array<T> (void) const;
+  Array<T> array_value (void) const;
 
   const T *data (void) const { return Array<T>::data (); }
 
--- a/liboctave/MDiagArray2.h
+++ b/liboctave/MDiagArray2.h
@@ -52,8 +52,6 @@
 
   MDiagArray2 (octave_idx_type r, octave_idx_type c, const T& val) : DiagArray2<T> (r, c, val) { }
 
-  MDiagArray2 (const dim_vector& dv) : DiagArray2<T> (dv) { }
-
   MDiagArray2 (const MDiagArray2<T>& a) : DiagArray2<T> (a) { }
 
   MDiagArray2 (const DiagArray2<T>& a) : DiagArray2<T> (a) { }
@@ -74,9 +72,9 @@
       return *this;
     }
 
-  operator MArray<T> () const
+  MArray<T> array_value () const
     {
-      return DiagArray2<T>::operator Array<T> ();
+      return DiagArray2<T>::array_value ();
     }
 
   octave_idx_type nnz (void) const
--- a/liboctave/dColVector.cc
+++ b/liboctave/dColVector.cc
@@ -146,23 +146,21 @@
 }
 
 ColumnVector
+ColumnVector::abs (void) const
+{
+  return do_mx_unary_map<double, double, std::abs> (*this);
+}
+
+ColumnVector
 real (const ComplexColumnVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  ColumnVector retval;
-  if (a_len > 0)
-    retval = ColumnVector (mx_inline_real_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<double, Complex> (a, mx_inline_real);
 }
 
 ColumnVector
 imag (const ComplexColumnVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  ColumnVector retval;
-  if (a_len > 0)
-    retval = ColumnVector (mx_inline_imag_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
 }
 
 // resize is the destructive equivalent for this one
--- a/liboctave/dColVector.h
+++ b/liboctave/dColVector.h
@@ -90,6 +90,8 @@
   double min (void) const;
   double max (void) const;
 
+  ColumnVector abs (void) const;
+
   // i/o
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const ColumnVector& a);
--- a/liboctave/dDiagMatrix.cc
+++ b/liboctave/dDiagMatrix.cc
@@ -142,32 +142,19 @@
 DiagMatrix
 DiagMatrix::abs (void) const
 {
-  DiagMatrix retval (rows (), cols ());
-  for (octave_idx_type i = 0; i < rows (); i++)
-    retval(i, i) = std::abs (elem (i, i));
-  return retval;
+  return DiagMatrix (diag ().abs (), rows (), columns ());
 }
 
 DiagMatrix
 real (const ComplexDiagMatrix& a)
 {
-  DiagMatrix retval;
-  octave_idx_type a_len = a.length ();
-  if (a_len > 0)
-    retval = DiagMatrix (mx_inline_real_dup (a.data (), a_len), a.rows (),
-                         a.cols ());
-  return retval;
+  return DiagMatrix (real (a.diag ()), a.rows (), a.cols ());
 }
 
 DiagMatrix
 imag (const ComplexDiagMatrix& a)
 {
-  DiagMatrix retval;
-  octave_idx_type a_len = a.length ();
-  if (a_len > 0)
-    retval = DiagMatrix (mx_inline_imag_dup (a.data (), a_len), a.rows (),
-                         a.cols ());
-  return retval;
+  return DiagMatrix (imag (a.diag ()), a.rows (), a.cols ());
 }
 
 Matrix
--- a/liboctave/dDiagMatrix.h
+++ b/liboctave/dDiagMatrix.h
@@ -54,9 +54,10 @@
   template <class U>
   DiagMatrix (const DiagArray2<U>& a) : MDiagArray2<double> (a) { }
 
-  explicit DiagMatrix (const RowVector& a) : MDiagArray2<double> (a) { }
+  explicit DiagMatrix (const Array<double>& a) : MDiagArray2<double> (a) { }
 
-  explicit DiagMatrix (const ColumnVector& a) : MDiagArray2<double> (a) { }
+  DiagMatrix (const Array<double>& a, octave_idx_type r, octave_idx_type c) 
+    : MDiagArray2<double> (a, r, c) { }
 
   DiagMatrix& operator = (const DiagMatrix& a)
     {
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -565,15 +565,13 @@
 Matrix
 real (const ComplexMatrix& a)
 {
-  return Matrix (mx_inline_real_dup (a.data (), a.length ()),
-                 a.rows (), a.cols ());
+  return do_mx_unary_op<double, Complex> (a, mx_inline_real);
 }
 
 Matrix
 imag (const ComplexMatrix& a)
 {
-  return Matrix (mx_inline_imag_dup (a.data (), a.length ()),
-                 a.rows (), a.cols ());
+  return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
 }
 
 Matrix
@@ -2822,8 +2820,7 @@
 Matrix
 Matrix::abs (void) const
 {
-  return Matrix (mx_inline_fabs_dup (data (), length ()),
-                 rows (), cols ());
+  return do_mx_unary_map<double, double, std::abs> (*this);
 }
 
 Matrix
--- a/liboctave/dRowVector.cc
+++ b/liboctave/dRowVector.cc
@@ -150,21 +150,13 @@
 RowVector
 real (const ComplexRowVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  RowVector retval;
-  if (a_len > 0)
-    retval = RowVector (mx_inline_real_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<double, Complex> (a, mx_inline_real);
 }
 
 RowVector
 imag (const ComplexRowVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  RowVector retval;
-  if (a_len > 0)
-    retval = RowVector (mx_inline_imag_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
 }
 
 RowVector
--- a/liboctave/fCColVector.cc
+++ b/liboctave/fCColVector.cc
@@ -231,14 +231,16 @@
   return MArray<FloatComplex>::transpose ();
 }
 
+FloatColumnVector
+FloatComplexColumnVector::abs (void) const
+{
+  return do_mx_unary_map<float, FloatComplex, std::abs> (*this);
+}
+
 FloatComplexColumnVector
 conj (const FloatComplexColumnVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  FloatComplexColumnVector retval;
-  if (a_len > 0)
-    retval = FloatComplexColumnVector (mx_inline_conj_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_map<FloatComplex, FloatComplex, std::conj> (a);
 }
 
 // resize is the destructive equivalent for this one
--- a/liboctave/fCColVector.h
+++ b/liboctave/fCColVector.h
@@ -121,6 +121,8 @@
   FloatComplex min (void) const;
   FloatComplex max (void) const;
 
+  FloatColumnVector abs (void) const;
+
   // i/o
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatComplexColumnVector& a);
--- a/liboctave/fCDiagMatrix.cc
+++ b/liboctave/fCDiagMatrix.cc
@@ -237,21 +237,13 @@
 FloatDiagMatrix
 FloatComplexDiagMatrix::abs (void) const
 {
-  FloatDiagMatrix retval (rows (), cols ());
-  for (octave_idx_type i = 0; i < rows (); i++)
-    retval(i, i) = std::abs (elem (i, i));
-  return retval;
+  return FloatDiagMatrix (diag ().abs (), rows (), columns ());
 }
 
 FloatComplexDiagMatrix
 conj (const FloatComplexDiagMatrix& a)
 {
-  FloatComplexDiagMatrix retval;
-  octave_idx_type a_len = a.length ();
-  if (a_len > 0)
-    retval = FloatComplexDiagMatrix (mx_inline_conj_dup (a.data (), a_len),
-                                a.rows (), a.cols ());
-  return retval;
+  return FloatComplexDiagMatrix (conj (a.diag ()), a.rows (), a.columns ());
 }
 
 // resize is the destructive analog for this one
--- a/liboctave/fCDiagMatrix.h
+++ b/liboctave/fCDiagMatrix.h
@@ -47,17 +47,14 @@
   FloatComplexDiagMatrix (octave_idx_type r, octave_idx_type c, const FloatComplex& val)
     : MDiagArray2<FloatComplex> (r, c, val) { }
 
-  explicit FloatComplexDiagMatrix (const FloatRowVector& a)
-    : MDiagArray2<FloatComplex> (FloatComplexRowVector (a)) { }
-
-  explicit FloatComplexDiagMatrix (const FloatComplexRowVector& a)
+  explicit FloatComplexDiagMatrix (const Array<FloatComplex>& a)
     : MDiagArray2<FloatComplex> (a) { }
 
-  explicit FloatComplexDiagMatrix (const FloatColumnVector& a)
-    : MDiagArray2<FloatComplex> (FloatComplexColumnVector (a)) { }
+  FloatComplexDiagMatrix (const Array<FloatComplex>& a, octave_idx_type r, octave_idx_type c) 
+    : MDiagArray2<FloatComplex> (a, r, c) { }
 
-  explicit FloatComplexDiagMatrix (const FloatComplexColumnVector& a)
-    : MDiagArray2<FloatComplex> (a) { }
+  explicit FloatComplexDiagMatrix (const Array<float>& a)
+    : MDiagArray2<FloatComplex> (Array<FloatComplex> (a)) { }
 
   explicit FloatComplexDiagMatrix (const FloatDiagMatrix& a);
 
--- a/liboctave/fCMatrix.cc
+++ b/liboctave/fCMatrix.cc
@@ -891,8 +891,7 @@
 FloatComplexMatrix
 conj (const FloatComplexMatrix& a)
 {
-  return FloatComplexMatrix (mx_inline_conj_dup (a.data (), a.length ()),
-                             a.rows (), a.cols ());
+  return do_mx_unary_map<FloatComplex, FloatComplex, std::conj> (a);
 }
 
 // resize is the destructive equivalent for this one
@@ -3236,8 +3235,7 @@
 
 FloatMatrix FloatComplexMatrix::abs (void) const
 {
-  return FloatMatrix (mx_inline_cabs_dup (data (), length ()),
-                      rows (), cols ());
+  return do_mx_unary_map<float, FloatComplex, std::abs> (*this);
 }
 
 FloatComplexMatrix
--- a/liboctave/fCRowVector.cc
+++ b/liboctave/fCRowVector.cc
@@ -232,11 +232,7 @@
 FloatComplexRowVector
 conj (const FloatComplexRowVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  FloatComplexRowVector retval;
-  if (a_len > 0)
-    retval = FloatComplexRowVector (mx_inline_conj_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_map<FloatComplex, FloatComplex, std::conj> (a);
 }
 
 // resize is the destructive equivalent for this one
--- a/liboctave/fColVector.cc
+++ b/liboctave/fColVector.cc
@@ -146,23 +146,21 @@
 }
 
 FloatColumnVector
+FloatColumnVector::abs (void) const
+{
+  return do_mx_unary_map<float, float, std::abs> (*this);
+}
+
+FloatColumnVector
 real (const FloatComplexColumnVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  FloatColumnVector retval;
-  if (a_len > 0)
-    retval = FloatColumnVector (mx_inline_real_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_real);
 }
 
 FloatColumnVector
 imag (const FloatComplexColumnVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  FloatColumnVector retval;
-  if (a_len > 0)
-    retval = FloatColumnVector (mx_inline_imag_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_imag);
 }
 
 // resize is the destructive equivalent for this one
--- a/liboctave/fColVector.h
+++ b/liboctave/fColVector.h
@@ -90,6 +90,8 @@
   float min (void) const;
   float max (void) const;
 
+  FloatColumnVector abs (void) const;
+
   // i/o
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatColumnVector& a);
--- a/liboctave/fDiagMatrix.cc
+++ b/liboctave/fDiagMatrix.cc
@@ -142,32 +142,19 @@
 FloatDiagMatrix
 FloatDiagMatrix::abs (void) const
 {
-  FloatDiagMatrix retval (rows (), cols ());
-  for (octave_idx_type i = 0; i < rows (); i++)
-    retval(i, i) = std::abs (elem (i, i));
-  return retval;
+  return FloatDiagMatrix (diag ().abs (), rows (), columns ());
 }
 
 FloatDiagMatrix
 real (const FloatComplexDiagMatrix& a)
 {
-  FloatDiagMatrix retval;
-  octave_idx_type a_len = a.length ();
-  if (a_len > 0)
-    retval = FloatDiagMatrix (mx_inline_real_dup (a.data (), a_len), a.rows (),
-                         a.cols ());
-  return retval;
+  return FloatDiagMatrix (real (a.diag ()), a.rows (), a.columns ());
 }
 
 FloatDiagMatrix
 imag (const FloatComplexDiagMatrix& a)
 {
-  FloatDiagMatrix retval;
-  octave_idx_type a_len = a.length ();
-  if (a_len > 0)
-    retval = FloatDiagMatrix (mx_inline_imag_dup (a.data (), a_len), a.rows (),
-                         a.cols ());
-  return retval;
+  return FloatDiagMatrix (imag (a.diag ()), a.rows (), a.columns ());
 }
 
 FloatMatrix
--- a/liboctave/fDiagMatrix.h
+++ b/liboctave/fDiagMatrix.h
@@ -54,9 +54,10 @@
   template <class U>
   FloatDiagMatrix (const DiagArray2<U>& a) : MDiagArray2<float> (a) { }
 
-  explicit FloatDiagMatrix (const FloatRowVector& a) : MDiagArray2<float> (a) { }
+  explicit FloatDiagMatrix (const Array<double>& a) : MDiagArray2<float> (a) { }
 
-  explicit FloatDiagMatrix (const FloatColumnVector& a) : MDiagArray2<float> (a) { }
+  FloatDiagMatrix (const Array<float>& a, octave_idx_type r, octave_idx_type c) 
+    : MDiagArray2<float> (a, r, c) { }
 
   FloatDiagMatrix& operator = (const FloatDiagMatrix& a)
     {
--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -564,15 +564,13 @@
 FloatMatrix
 real (const FloatComplexMatrix& a)
 {
-  return FloatMatrix (mx_inline_real_dup (a.data (), a.length ()),
-                      a.rows (), a.cols ());
+  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_real);
 }
 
 FloatMatrix
 imag (const FloatComplexMatrix& a)
 {
-  return FloatMatrix (mx_inline_imag_dup (a.data (), a.length ()),
-                      a.rows (), a.cols ());
+  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_imag);
 }
 
 FloatMatrix
@@ -2821,8 +2819,7 @@
 FloatMatrix
 FloatMatrix::abs (void) const
 {
-  return FloatMatrix (mx_inline_fabs_dup (data (), length ()),
-                      rows (), cols ());
+  return do_mx_unary_map<float, float, std::abs> (*this);
 }
 
 FloatMatrix
--- a/liboctave/fRowVector.cc
+++ b/liboctave/fRowVector.cc
@@ -150,21 +150,13 @@
 FloatRowVector
 real (const FloatComplexRowVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  FloatRowVector retval;
-  if (a_len > 0)
-    retval = FloatRowVector (mx_inline_real_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_real);
 }
 
 FloatRowVector
 imag (const FloatComplexRowVector& a)
 {
-  octave_idx_type a_len = a.length ();
-  FloatRowVector retval;
-  if (a_len > 0)
-    retval = FloatRowVector (mx_inline_imag_dup (a.data (), a_len), a_len);
-  return retval;
+  return do_mx_unary_op<float, FloatComplex> (a, mx_inline_imag);
 }
 
 FloatRowVector
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2010-02-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* src/xpow.cc: Update.
+	* src/ov-base-diag.cc: Update.
+
 2010-02-25  Jaroslav Hajek  <highegg@gmail.com>
 
 	* DLD-FUNCTIONS/str2double.cc (single_number): New helper func.
--- a/src/ov-base-diag.cc
+++ b/src/ov-base-diag.cc
@@ -82,8 +82,7 @@
 
       if (idx0.is_scalar () && idx1.is_scalar ())
         {
-          // FIXME: the proxy mechanism of DiagArray2 causes problems here.
-          retval = el_type (matrix.checkelem (idx0(0), idx1(0)));
+          retval = matrix.checkelem (idx0(0), idx1(0));
         }
       else
         {
@@ -131,7 +130,7 @@
                     && i0(0) < matrix.rows () && i1(0) < matrix.cols ()
                     && chk_valid_scalar (rhs, val))
                   {
-                    matrix (i0(0), i1(0)) = val;                    
+                    matrix.dgelem (i0(0)) = val;                    
                     retval = this;
                     this->count++;
                     // invalidate cache
--- a/src/xpow.cc
+++ b/src/xpow.cc
@@ -281,14 +281,14 @@
         {
           DiagMatrix r (nr, nc);
           for (octave_idx_type i = 0; i < nc; i++)
-            r(i, i) = std::pow (a(i, i), b);
+            r.dgelem (i) = std::pow (a.dgelem (i), b);
           retval = r;
         }
       else
         {
           ComplexDiagMatrix r (nr, nc);
           for (octave_idx_type i = 0; i < nc; i++)
-            r(i, i) = std::pow (static_cast<Complex> (a(i, i)), b);
+            r.dgelem (i) = std::pow (static_cast<Complex> (a.dgelem (i)), b);
           retval = r;
         }
     }
@@ -1694,14 +1694,14 @@
         {
           FloatDiagMatrix r (nr, nc);
           for (octave_idx_type i = 0; i < nc; i++)
-            r(i, i) = std::pow (a(i, i), b);
+            r.dgelem (i) = std::pow (a.dgelem (i), b);
           retval = r;
         }
       else
         {
           FloatComplexDiagMatrix r (nr, nc);
           for (octave_idx_type i = 0; i < nc; i++)
-            r(i, i) = std::pow (static_cast<FloatComplex> (a(i, i)), b);
+            r.dgelem (i) = std::pow (static_cast<FloatComplex> (a.dgelem (i)), b);
           retval = r;
         }
     }