--- a/liboctave/CColVector.h
+++ b/liboctave/CColVector.h
@@ -41,6 +41,8 @@
 
   explicit ComplexColumnVector (octave_idx_type n) : MArray<Complex> (n) { }
 
+  explicit ComplexColumnVector (const dim_vector& dv) : MArray<Complex> (dv) { }
+
   ComplexColumnVector (octave_idx_type n, const Complex& val)
     : MArray<Complex> (n, val) { }
 

--- a/liboctave/CNDArray.cc
+++ b/liboctave/CNDArray.cc
@@ -1083,8 +1083,9 @@
 ComplexNDArray& operator *= (ComplexNDArray& a, double s)
 {
   if (a.is_shared ())
-    return a = a * s;
-  DO_VS_OP2 (Complex, a, *=, s)
+    a = a * s;
+  else
+    do_ms_inplace_op<ComplexNDArray, double> (a, s, mx_inline_mul2);
   return a;
 }
 
@@ -1092,7 +1093,8 @@
 {
   if (a.is_shared ())
     return a = a / s;
-  DO_VS_OP2 (Complex, a, /=, s)
+  else
+    do_ms_inplace_op<ComplexNDArray, double> (a, s, mx_inline_div2);
   return a;
 }
 

--- a/liboctave/CRowVector.h
+++ b/liboctave/CRowVector.h
@@ -40,6 +40,8 @@
 
   explicit ComplexRowVector (octave_idx_type n) : MArray<Complex> (n) { }
 
+  explicit ComplexRowVector (const dim_vector& dv) : MArray<Complex> (dv) { }
+
   ComplexRowVector (octave_idx_type n, const Complex& val) : MArray<Complex> (n, val) { }
 
   ComplexRowVector (const ComplexRowVector& a) : MArray<Complex> (a) { }

--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,88 @@
+2009-09-04  Jaroslav Hajek  <highegg@gmail.com>
+
+	* mx-inlines.cc (DEFMXBOOLOPEQ): New macro.
+	(mx_inline_and2, mx_inline_or2): New loops.
+
+	* MArray.h (MArray<T>::MArray (const dim_vector&)): New constructor.
+	(MArray<T>::MArray (const Array<T>&)): Ensure column vector.
+	(MArray<T>::resize): New method.
+	* DiagArray2.h (DiagArray2<T>::DiagArray2 (const dim_vector&)): New
+	constructor.
+	* MDiagArray2.h (MDiagArray2<T>::MDiagArray2 (const dim_vector&)): New
+	constructor.
+	* dColVector.h (ColumnVector::ColumnVector (const dim_vector&)): New
+	constructor.
+	* fColVector.h (FloatColumnVector::FloatColumnVector (const
+	dim_vector&)): New constructor.
+	* CColVector.h (ComplexColumnVector::ComplexColumnVector (const
+	dim_vector&)): New constructor.
+	* fCColVector.h (FloatComplexColumnVector::FloatComplexColumnVector
+	(const dim_vector&)): New constructor.
+	* dRowVector.h (RowVector::RowVector (const dim_vector&)): New
+	constructor.
+	* fRowVector.h (FloatRowVector::FloatRowVector (const dim_vector&)):
+	New constructor.
+	* CRowVector.h (ComplexRowVector::ComplexRowVector (const
+	dim_vector&)): New constructor.
+	* fCRowVector.h (FloatComplexRowVector::FloatComplexRowVector (const
+	dim_vector&)): New constructor.
+
+	* CNDArray.cc (operator *= (const ComplexNDArray, double)): Simplify.
+	* fCNDArray.cc (operator *= (const FloatComplexNDArray, float)): Simplify.
+
+	* MArray.cc
+	(operator += (MArray<T>&, const T&)): Simplify.
+	(operator -= (MArray<T>&, const T&)): Simplify.
+	(operator *= (MArray<T>&, const T&)): Simplify.
+	(operator /= (MArray<T>&, const T&)): Simplify.
+	(operator += (MArray<T>&, const MArray<T>&)): Simplify.
+	(operator -= (MArray<T>&, const MArray<T>&)): Simplify.
+	(product (MArray<T>&, const MArray<T>&)): Simplify.
+	(quotient (MArray<T>&, const MArray<T>&)): Simplify.
+	(MARRAY_AS_OP): Simplify.
+	(MARRAY_SA_OP): Simplify.
+	(MARRAY_AA_OP): Simplify.
+	(operator - (const MArray<T>&)): Simplify.
+
+	* MArray2.cc
+	(operator += (MArray2<T>&, const T&)): Simplify.
+	(operator -= (MArray2<T>&, const T&)): Simplify.
+	(operator *= (MArray2<T>&, const T&)): Simplify.
+	(operator /= (MArray2<T>&, const T&)): Simplify.
+	(operator += (MArray2<T>&, const MArray<T>&)): Simplify.
+	(operator -= (MArray2<T>&, const MArray<T>&)): Simplify.
+	(product (MArray2<T>&, const MArray<T>&)): Simplify.
+	(quotient (MArray2<T>&, const MArray<T>&)): Simplify.
+	(MARRAY_A2S_OP): Simplify.
+	(MARRAY_SA2_OP): Simplify.
+	(MARRAY_A2A2_OP): Simplify.
+	(operator - (const MArray2<T>&)): Simplify.
+
+	* MArrayN.cc
+	(operator += (MArrayN<T>&, const T&)): Simplify.
+	(operator -= (MArrayN<T>&, const T&)): Simplify.
+	(operator *= (MArrayN<T>&, const T&)): Simplify.
+	(operator /= (MArrayN<T>&, const T&)): Simplify.
+	(operator += (MArrayN<T>&, const MArray<T>&)): Simplify.
+	(operator -= (MArrayN<T>&, const MArray<T>&)): Simplify.
+	(product (MArrayN<T>&, const MArray<T>&)): Simplify.
+	(quotient (MArrayN<T>&, const MArray<T>&)): Simplify.
+	(MARRAY_NDS_OP): Simplify.
+	(MARRAY_SND_OP): Simplify.
+	(MARRAY_NDND_OP): Simplify.
+	(operator - (const MArrayN<T>&)): Simplify.
+
+	* MDiagArray2.cc
+	(operator += (MDiagArray2<T>&, const MDiagArray2<T>&)): Simplify.
+	(operator -= (MDiagArray2<T>&, const MDiagArray2<T>&)): Simplify.
+	(operator *= (MDiagArray2<T>&, T)): Simplify.
+	(operator /= (MDiagArray2<T>&, T)): Simplify.
+	(MARRAY_DAS_OP, MARRAY_SDA_OP, MARRAY_DADA_OP): Simplify.
+	* boolNDArray.cc (mx_el_and_assign, mx_el_or_assign): Simplify.
+
+	* MArray-defs.h (DO_VS_OP, DO_SV_OP, DO_VV_OP,
+	DO_VS_OP2, DO_VV_OP2, NEGV): Remove.
+
 2009-09-03  Jaroslav Hajek  <highegg@gmail.com>
 
 	* mx-inlines.cc (DEFMXUNOPEQ): New macro.

--- a/liboctave/DiagArray2.h
+++ b/liboctave/DiagArray2.h
@@ -105,6 +105,13 @@
   DiagArray2 (octave_idx_type r, octave_idx_type c, const T& val) 
     : Array<T> (std::min (r, c), val), d1 (r), d2 (c) { }
 
+  DiagArray2 (const dim_vector& dv)
+    : Array<T> (std::min (dv(0), dv(1))), d1 (dv(0)), d2 (dv(0))
+    {
+      if (dv.length () != 2)
+	(*current_liboctave_error_handler) ("too many dimensions");
+    }
+
   DiagArray2 (const Array<T>& a) 
     : Array<T> (a), d1 (a.numel ()), d2 (a.numel ()) { }
 

--- a/liboctave/MArray-defs.h
+++ b/liboctave/MArray-defs.h
@@ -25,54 +25,7 @@
 #if !defined (octave_MArray_defs_h)
 #define octave_MArray_defs_h 1
 
-// Nothing like a little CPP abuse to brighten everyone's day.
-
-#define DO_VS_OP(r, l, v, OP, s) \
-  if (l > 0) \
-    { \
-      for (octave_idx_type i = 0; i < l; i++) \
-	r[i] = v[i] OP s; \
-    }
-
-#define DO_SV_OP(r, l, s, OP, v) \
-  if (l > 0) \
-    { \
-      for (octave_idx_type i = 0; i < l; i++) \
-	r[i] = s OP v[i]; \
-    }
-
-#define DO_VV_OP(r, l, x, OP, y) \
-  if (l > 0) \
-    { \
-      for (octave_idx_type i = 0; i < l; i++) \
-	r[i] = x[i] OP y[i]; \
-    }
-
-#define NEG_V(r, l, x) \
-  if (l > 0) \
-    { \
-      for (octave_idx_type i = 0; i < l; i++) \
-	r[i] = -x[i]; \
-    }
-
-#define DO_VS_OP2(T, a, OP, s) \
-  octave_idx_type l = a.length (); \
-  if (l > 0) \
-    { \
-      T *tmp = a.fortran_vec (); \
-      for (octave_idx_type i = 0; i < l; i++) \
-	tmp[i] OP s; \
-    }
-
-#define DO_VV_OP2(T, a, OP, b) \
-  do \
-    { \
-      T *a_tmp = a.fortran_vec (); \
-      const T *b_tmp = b.data (); \
-      for (octave_idx_type i = 0; i < l; i++) \
-	a_tmp[i] OP b_tmp[i]; \
-    } \
-  while (0)
+#include "mx-inlines.cc"
 
 // Instantiate the OP= operators.
 #define MARRAY_OP_ASSIGN_DEFS(A_T, E_T, RHS_T, API) \

--- a/liboctave/MArray.cc
+++ b/liboctave/MArray.cc
@@ -117,8 +117,9 @@
 operator += (MArray<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a + s;
-  DO_VS_OP2 (T, a, +=, s)
+    a = a + s;
+  else
+    do_ms_inplace_op<MArray<T>, T> (a, s, mx_inline_add2);
   return a;
 }
 
@@ -127,8 +128,9 @@
 operator -= (MArray<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a - s;
-  DO_VS_OP2 (T, a, -=, s)
+    a = a - s;
+  else
+    do_ms_inplace_op<MArray<T>, T> (a, s, mx_inline_sub2);
   return a;
 }
 
@@ -137,8 +139,9 @@
 operator *= (MArray<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a * s;
-  DO_VS_OP2 (T, a, *=, s)
+    a = a * s;
+  else
+    do_ms_inplace_op<MArray<T>, T> (a, s, mx_inline_mul2);
   return a;
 }
 
@@ -147,8 +150,9 @@
 operator /= (MArray<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a / s;
-  DO_VS_OP2 (T, a, /=, s)
+    a = a / s;
+  else
+    do_ms_inplace_op<MArray<T>, T> (a, s, mx_inline_div2);
   return a;
 }
 
@@ -159,16 +163,9 @@
 operator += (MArray<T>& a, const MArray<T>& b)
 {
   if (a.is_shared ())
-    return a = a + b;
-  octave_idx_type l = a.length ();
-  if (l > 0)
-    {
-      octave_idx_type bl = b.length ();
-      if (l != bl)
-	gripe_nonconformant ("operator +=", l, bl);
-      else
-	DO_VV_OP2 (T, a, +=, b);
-    }
+    a = a + b;
+  else
+    do_mm_inplace_op<MArray<T>, MArray<T> > (a, b, mx_inline_add2, "+=");
   return a;
 }
 
@@ -177,16 +174,9 @@
 operator -= (MArray<T>& a, const MArray<T>& b)
 {
   if (a.is_shared ())
-    return a = a - b;
-  octave_idx_type l = a.length ();
-  if (l > 0)
-    {
-      octave_idx_type bl = b.length ();
-      if (l != bl)
-	gripe_nonconformant ("operator -=", l, bl);
-      else
-	DO_VV_OP2 (T, a, -=, b);
-    }
+    a = a - b;
+  else
+    do_mm_inplace_op<MArray<T>, MArray<T> > (a, b, mx_inline_sub2, "-=");
   return a;
 }
 
@@ -196,15 +186,8 @@
 {
   if (a.is_shared ())
     return a = product (a, b);
-  octave_idx_type l = a.length ();
-  if (l > 0)
-    {
-      octave_idx_type bl = b.length ();
-      if (l != bl)
-	gripe_nonconformant ("operator .*=", l, bl);
-      else
-	DO_VV_OP2 (T, a, *=, b);
-    }
+  else
+    do_mm_inplace_op<MArray<T>, MArray<T> > (a, b, mx_inline_mul2, ".*=");
   return a;
 }
 
@@ -214,86 +197,55 @@
 {
   if (a.is_shared ())
     return a = quotient (a, b);
-  octave_idx_type l = a.length ();
-  if (l > 0)
-    {
-      octave_idx_type bl = b.length ();
-      if (l != bl)
-	gripe_nonconformant ("operator ./=", l, bl);
-      else
-	DO_VV_OP2 (T, a, /=, b);
-    }
+  else
+    do_mm_inplace_op<MArray<T>, MArray<T> > (a, b, mx_inline_div2, "./=");
   return a;
 }
 
 // Element by element MArray by scalar ops.
 
-#define MARRAY_AS_OP(OP) \
+#define MARRAY_AS_OP(OP, FN) \
   template <class T> \
   MArray<T> \
   operator OP (const MArray<T>& a, const T& s) \
   { \
-    MArray<T> result (a.length ()); \
-    T *r = result.fortran_vec (); \
-    octave_idx_type l = a.length (); \
-    const T *v = a.data (); \
-    DO_VS_OP (r, l, v, OP, s); \
-    return result; \
+    return do_ms_binary_op<MArray<T>, MArray<T>, T> (a, s, FN); \
   }
 
-MARRAY_AS_OP (+)
-MARRAY_AS_OP (-)
-MARRAY_AS_OP (*)
-MARRAY_AS_OP (/)
+MARRAY_AS_OP (+, mx_inline_add)
+MARRAY_AS_OP (-, mx_inline_sub)
+MARRAY_AS_OP (*, mx_inline_mul)
+MARRAY_AS_OP (/, mx_inline_div)
 
 // Element by element scalar by MArray ops.
 
-#define MARRAY_SA_OP(OP) \
+#define MARRAY_SA_OP(OP, FN) \
   template <class T> \
   MArray<T> \
   operator OP (const T& s, const MArray<T>& a) \
   { \
-    MArray<T> result (a.length ()); \
-    T *r = result.fortran_vec (); \
-    octave_idx_type l = a.length (); \
-    const T *v = a.data (); \
-    DO_SV_OP (r, l, s, OP, v); \
-    return result; \
+    return do_sm_binary_op<MArray<T>, T, MArray<T> > (s, a, FN); \
   }
 
-MARRAY_SA_OP(+)
-MARRAY_SA_OP(-)
-MARRAY_SA_OP(*)
-MARRAY_SA_OP(/)
+MARRAY_SA_OP(+, mx_inline_add)
+MARRAY_SA_OP(-, mx_inline_sub)
+MARRAY_SA_OP(*, mx_inline_mul)
+MARRAY_SA_OP(/, mx_inline_div)
 
 // Element by element MArray by MArray ops.
 
-#define MARRAY_AA_OP(FCN, OP) \
+#define MARRAY_AA_OP(FCN, OP, FN) \
   template <class T> \
   MArray<T> \
   FCN (const MArray<T>& a, const MArray<T>& b) \
   { \
-    octave_idx_type l = a.length (); \
-    octave_idx_type bl = b.length (); \
-    if (l != bl) \
-      { \
-	gripe_nonconformant (#FCN, l, bl); \
-	return MArray<T> (); \
-      } \
-    if (l == 0) \
-      return MArray<T> (); \
-    MArray<T> result (l); \
-    T *r = result.fortran_vec (); \
-    const T *x = a.data (); \
-    const T *y = b.data (); \
-    DO_VV_OP (r, l, x, OP, y); \
-    return result; \
+    return do_mm_binary_op<MArray<T>, MArray<T>, MArray<T> > (a, b, FN, #FCN); \
   }
 
-MARRAY_AA_OP (operator +, +)
-MARRAY_AA_OP (operator -, -)
-MARRAY_AA_OP (product,    *)
-MARRAY_AA_OP (quotient,   /)
+MARRAY_AA_OP (operator +, +, mx_inline_add)
+MARRAY_AA_OP (operator -, -, mx_inline_sub)
+MARRAY_AA_OP (product,    *, mx_inline_mul)
+MARRAY_AA_OP (quotient,   /, mx_inline_div)
 
 // Unary MArray ops.
 
@@ -308,12 +260,7 @@
 MArray<T>
 operator - (const MArray<T>& a)
 {
-  octave_idx_type l = a.length ();
-  MArray<T> result (l);
-  T *r = result.fortran_vec ();
-  const T *x = a.data ();
-  NEG_V (r, l, x);
-  return result;
+  return do_mx_unary_op<MArray<T>, MArray<T> > (a, mx_inline_uminus); 
 }
 
 /*

--- a/liboctave/MArray.h
+++ b/liboctave/MArray.h
@@ -53,7 +53,10 @@
 
   MArray (const MArray<T>& a) : Array<T> (a) { }
 
-  MArray (const Array<T>& a) : Array<T> (a) { }
+  // FIXME: kluge.
+  MArray (const Array<T>& a) : Array<T> (a, dim_vector (a.length ())) { }
+
+  MArray (const dim_vector& dv) : Array<T> (dv) { }
 
   ~MArray (void) { }
 
@@ -63,6 +66,10 @@
       return *this;
     }
 
+  // FIXME: kluge again. This design really sucks.
+  void resize (octave_idx_type n, const T& val = Array<T>::resize_fill_value ())
+    { Array<T>::resize_fill (n, 1, val); }
+
   MArray<T> transpose (void) const { return Array<T>::transpose (); }
   MArray<T> hermitian (T (*fcn) (const T&) = 0) const { return Array<T>::hermitian (fcn); }
 

--- a/liboctave/MArray2.cc
+++ b/liboctave/MArray2.cc
@@ -40,8 +40,9 @@
 operator += (MArray2<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a + s;
-  DO_VS_OP2 (T, a, +=, s)
+    a = a + s;
+  else
+    do_ms_inplace_op<MArray2<T>, T> (a, s, mx_inline_add2);
   return a;
 }
 
@@ -50,8 +51,9 @@
 operator -= (MArray2<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a - s;
-  DO_VS_OP2 (T, a, -=, s)
+    a = a - s;
+  else
+    do_ms_inplace_op<MArray2<T>, T> (a, s, mx_inline_sub2);
   return a;
 }
 
@@ -60,8 +62,9 @@
 operator *= (MArray2<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a * s;
-  DO_VS_OP2 (T, a, *=, s)
+    a = a * s;
+  else
+    do_ms_inplace_op<MArray2<T>, T> (a, s, mx_inline_mul2);
   return a;
 }
 
@@ -70,8 +73,9 @@
 operator /= (MArray2<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a / s;
-  DO_VS_OP2 (T, a, /=, s)
+    a = a / s;
+  else
+    do_ms_inplace_op<MArray2<T>, T> (a, s, mx_inline_div2);
   return a;
 }
 
@@ -82,21 +86,9 @@
 operator += (MArray2<T>& a, const MArray2<T>& b)
 {
   if (a.is_shared ())
-    return a = a + b;
-  octave_idx_type r = a.rows ();
-  octave_idx_type c = a.cols ();
-  octave_idx_type br = b.rows ();
-  octave_idx_type bc = b.cols ();
-  if (r != br || c != bc)
-    gripe_nonconformant ("operator +=", r, c, br, bc);
+    a = a + b;
   else
-    {
-      if (r > 0 && c > 0)
-	{
-	  octave_idx_type l = a.length ();
-	  DO_VV_OP2 (T, a, +=, b);
-	}
-    }
+    do_mm_inplace_op<MArray2<T>, MArray2<T> > (a, b, mx_inline_add2, "+=");
   return a;
 }
 
@@ -105,21 +97,9 @@
 operator -= (MArray2<T>& a, const MArray2<T>& b)
 {
   if (a.is_shared ())
-    return a = a - b;
-  octave_idx_type r = a.rows ();
-  octave_idx_type c = a.cols ();
-  octave_idx_type br = b.rows ();
-  octave_idx_type bc = b.cols ();
-  if (r != br || c != bc)
-    gripe_nonconformant ("operator -=", r, c, br, bc);
+    a = a - b;
   else
-    {
-      if (r > 0 && c > 0)
-	{
-	  octave_idx_type l = a.length ();
-	  DO_VV_OP2 (T, a, -=, b);
-	}
-    }
+    do_mm_inplace_op<MArray2<T>, MArray2<T> > (a, b, mx_inline_sub2, "-=");
   return a;
 }
 
@@ -130,20 +110,8 @@
 {
   if (a.is_shared ())
     return a = product (a, b);
-  octave_idx_type r = a.rows ();
-  octave_idx_type c = a.cols ();
-  octave_idx_type br = b.rows ();
-  octave_idx_type bc = b.cols ();
-  if (r != br || c != bc)
-    gripe_nonconformant ("operator .*=", r, c, br, bc);
   else
-    {
-      if (r > 0 && c > 0)
-	{
-	  octave_idx_type l = a.length ();
-	  DO_VV_OP2 (T, a, *=, b);
-	}
-    }
+    do_mm_inplace_op<MArray2<T>, MArray2<T> > (a, b, mx_inline_mul2, ".*=");
   return a;
 }
 
@@ -153,94 +121,55 @@
 {
   if (a.is_shared ())
     return a = quotient (a, b);
-  octave_idx_type r = a.rows ();
-  octave_idx_type c = a.cols ();
-  octave_idx_type br = b.rows ();
-  octave_idx_type bc = b.cols ();
-  if (r != br || c != bc)
-    gripe_nonconformant ("operator ./=", r, c, br, bc);
   else
-    {
-      if (r > 0 && c > 0)
-	{
-	  octave_idx_type l = a.length ();
-	  DO_VV_OP2 (T, a, /=, b);
-	}
-    }
+    do_mm_inplace_op<MArray2<T>, MArray2<T> > (a, b, mx_inline_div2, "./=");
   return a;
 }
 
 // Element by element MArray2 by scalar ops.
 
-#define MARRAY_A2S_OP(OP) \
+#define MARRAY_A2S_OP(OP, FN) \
   template <class T> \
   MArray2<T> \
   operator OP (const MArray2<T>& a, const T& s) \
   { \
-    MArray2<T> result (a.rows (), a.cols ()); \
-    T *r = result.fortran_vec (); \
-    octave_idx_type l = a.length (); \
-    const T *v = a.data (); \
-    DO_VS_OP (r, l, v, OP, s); \
-    return result; \
+    return do_ms_binary_op<MArray2<T>, MArray2<T>, T> (a, s, FN); \
   }
 
-MARRAY_A2S_OP (+)
-MARRAY_A2S_OP (-)
-MARRAY_A2S_OP (*)
-MARRAY_A2S_OP (/)
+MARRAY_A2S_OP (+, mx_inline_add)
+MARRAY_A2S_OP (-, mx_inline_sub)
+MARRAY_A2S_OP (*, mx_inline_mul)
+MARRAY_A2S_OP (/, mx_inline_div)
 
 // Element by element scalar by MArray2 ops.
 
-#define MARRAY_SA2_OP(OP) \
+#define MARRAY_SA2_OP(OP, FN) \
   template <class T> \
   MArray2<T> \
   operator OP (const T& s, const MArray2<T>& a) \
   { \
-    MArray2<T> result (a.rows (), a.cols ()); \
-    T *r = result.fortran_vec (); \
-    octave_idx_type l = a.length (); \
-    const T *v = a.data (); \
-    DO_SV_OP (r, l, s, OP, v); \
-    return result; \
+    return do_sm_binary_op<MArray2<T>, T, MArray2<T> > (s, a, FN); \
   }
 
-MARRAY_SA2_OP (+)
-MARRAY_SA2_OP (-)
-MARRAY_SA2_OP (*)
-MARRAY_SA2_OP (/)
+MARRAY_SA2_OP (+, mx_inline_add)
+MARRAY_SA2_OP (-, mx_inline_sub)
+MARRAY_SA2_OP (*, mx_inline_mul)
+MARRAY_SA2_OP (/, mx_inline_div)
 
 // Element by element MArray2 by MArray2 ops.
 
-#define MARRAY_A2A2_OP(FCN, OP) \
+#define MARRAY_A2A2_OP(FCN, OP, FN) \
   template <class T> \
   MArray2<T> \
   FCN (const MArray2<T>& a, const MArray2<T>& b) \
   { \
-    octave_idx_type a_nr = a.rows (); \
-    octave_idx_type a_nc = a.cols (); \
-    octave_idx_type b_nr = b.rows (); \
-    octave_idx_type b_nc = b.cols (); \
-    if (a_nr != b_nr || a_nc != b_nc) \
-      { \
-        gripe_nonconformant (#FCN, a_nr, a_nc, b_nr, b_nc); \
-	return MArray2<T> (); \
-      } \
-    if (a_nr == 0 || a_nc == 0) \
-      return MArray2<T> (a_nr, a_nc); \
-    octave_idx_type l = a.length (); \
-    MArray2<T> result (a_nr, a_nc); \
-    T *r = result.fortran_vec (); \
-    const T *x = a.data (); \
-    const T *y = b.data (); \
-    DO_VV_OP (r, l, x, OP, y); \
-    return result; \
+    return do_mm_binary_op<MArray2<T>, MArray2<T>, MArray2<T> > (a, b, FN, #FCN); \
   }
 
-MARRAY_A2A2_OP (operator +, +)
-MARRAY_A2A2_OP (operator -, -)
-MARRAY_A2A2_OP (product,    *)
-MARRAY_A2A2_OP (quotient,   /)
+MARRAY_A2A2_OP (operator +, +, mx_inline_add)
+MARRAY_A2A2_OP (operator -, -, mx_inline_sub)
+MARRAY_A2A2_OP (product,    *, mx_inline_mul)
+MARRAY_A2A2_OP (quotient,   /, mx_inline_div)
 
 // Unary MArray2 ops.
 
@@ -255,12 +184,7 @@
 MArray2<T>
 operator - (const MArray2<T>& a)
 {
-  octave_idx_type l = a.length ();
-  MArray2<T> result (a.rows (), a.cols ());
-  T *r = result.fortran_vec ();
-  const T *x = a.data ();
-  NEG_V (r, l, x);
-  return result;
+  return do_mx_unary_op<MArray2<T>, MArray2<T> > (a, mx_inline_uminus); 
 }
 
 /*

--- a/liboctave/MArrayN.cc
+++ b/liboctave/MArrayN.cc
@@ -49,8 +49,9 @@
 operator += (MArrayN<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a + s;
-  DO_VS_OP2 (T, a, +=, s)
+    a = a + s;
+  else
+    do_ms_inplace_op<MArrayN<T>, T> (a, s, mx_inline_add2);
   return a;
 }
 
@@ -59,8 +60,9 @@
 operator -= (MArrayN<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a - s;
-  DO_VS_OP2 (T, a, -=, s)
+    a = a - s;
+  else
+    do_ms_inplace_op<MArrayN<T>, T> (a, s, mx_inline_sub2);
   return a;
 }
 
@@ -69,8 +71,9 @@
 operator *= (MArrayN<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a * s;
-  DO_VS_OP2 (T, a, *=, s)
+    a = a * s;
+  else
+    do_ms_inplace_op<MArrayN<T>, T> (a, s, mx_inline_mul2);
   return a;
 }
 
@@ -79,8 +82,9 @@
 operator /= (MArrayN<T>& a, const T& s)
 {
   if (a.is_shared ())
-    return a = a / s;
-  DO_VS_OP2 (T, a, /=, s)
+    a = a / s;
+  else
+    do_ms_inplace_op<MArrayN<T>, T> (a, s, mx_inline_div2);
   return a;
 }
 
@@ -91,21 +95,9 @@
 operator += (MArrayN<T>& a, const MArrayN<T>& b)
 {
   if (a.is_shared ())
-    return a = a + b;
-
-  dim_vector a_dims = a.dims ();
-  dim_vector b_dims = b.dims ();
-
-  if (a_dims != b_dims)
-    gripe_nonconformant ("operator +=", a_dims, b_dims);
-  else 
-    {
-      octave_idx_type l = a.length ();
-
-      if (l > 0)
-        DO_VV_OP2 (T, a, +=, b);
-    }
-
+    a = a + b;
+  else
+    do_mm_inplace_op<MArrayN<T>, MArrayN<T> > (a, b, mx_inline_add2, "+=");
   return a;
 }
 
@@ -114,44 +106,21 @@
 operator -= (MArrayN<T>& a, const MArrayN<T>& b)
 {
   if (a.is_shared ())
-    return a = a - b;
-
-  dim_vector a_dims = a.dims ();
-  dim_vector b_dims = b.dims ();
-
-  if (a_dims != b_dims)
-    gripe_nonconformant ("operator -=", a_dims, b_dims);
-  else 
-    {
-      octave_idx_type l = a.length ();
-
-      if (l > 0)
-        DO_VV_OP2 (T, a, -=, b);
-    }
-
+    a = a - b;
+  else
+    do_mm_inplace_op<MArrayN<T>, MArrayN<T> > (a, b, mx_inline_sub2, "-=");
   return a;
 }
 
+
 template <class T>
 MArrayN<T>&
 product_eq (MArrayN<T>& a, const MArrayN<T>& b)
 {
   if (a.is_shared ())
     return a = product (a, b);
-
-  dim_vector a_dims = a.dims ();
-  dim_vector b_dims = b.dims ();
-
-  if (a_dims != b_dims)
-    gripe_nonconformant ("operator .*=", a_dims, b_dims);
-  else 
-    {
-      octave_idx_type l = a.length ();
-
-      if (l > 0)
-        DO_VV_OP2 (T, a, *=, b);
-    }
-
+  else
+    do_mm_inplace_op<MArrayN<T>, MArrayN<T> > (a, b, mx_inline_mul2, ".*=");
   return a;
 }
 
@@ -161,104 +130,55 @@
 {
   if (a.is_shared ())
     return a = quotient (a, b);
-
-  dim_vector a_dims = a.dims ();
-  dim_vector b_dims = b.dims ();
-
-  if (a_dims != b_dims)
-    gripe_nonconformant ("operator ./=", a_dims, b_dims);
-  else 
-    {
-      octave_idx_type l = a.length ();
-
-      if (l > 0)
-        DO_VV_OP2 (T, a, /=, b);
-    }
-
+  else
+    do_mm_inplace_op<MArrayN<T>, MArrayN<T> > (a, b, mx_inline_div2, "./=");
   return a;
 }
 
 // Element by element MArrayN by scalar ops.
 
-#define MARRAYN_NDS_OP(OP) \
+#define MARRAY_NDS_OP(OP, FN) \
   template <class T> \
   MArrayN<T> \
   operator OP (const MArrayN<T>& a, const T& s) \
-    { \
-      MArrayN<T> result (a.dims ()); \
-      T *r = result.fortran_vec (); \
-      octave_idx_type l = a.length (); \
-      const T *v = a.data (); \
-      DO_VS_OP (r, l, v, OP, s); \
-      return result; \
-    }
+  { \
+    return do_ms_binary_op<MArrayN<T>, MArrayN<T>, T> (a, s, FN); \
+  }
 
-MARRAYN_NDS_OP (+)
-MARRAYN_NDS_OP (-)
-MARRAYN_NDS_OP (*)
-MARRAYN_NDS_OP (/)
+MARRAY_NDS_OP (+, mx_inline_add)
+MARRAY_NDS_OP (-, mx_inline_sub)
+MARRAY_NDS_OP (*, mx_inline_mul)
+MARRAY_NDS_OP (/, mx_inline_div)
 
-// Element by element MArrayN by scalar ops.
+// Element by element scalar by MArrayN ops.
 
-#define MARRAYN_SND_OP(OP) \
+#define MARRAY_SND_OP(OP, FN) \
   template <class T> \
   MArrayN<T> \
   operator OP (const T& s, const MArrayN<T>& a) \
   { \
-    MArrayN<T> result (a.dims ()); \
-    T *r = result.fortran_vec (); \
-    octave_idx_type l = a.length (); \
-    const T *v = a.data (); \
-    DO_SV_OP (r, l, s, OP, v); \
-    return result; \
+    return do_sm_binary_op<MArrayN<T>, T, MArrayN<T> > (s, a, FN); \
   }
 
-MARRAYN_SND_OP (+)
-MARRAYN_SND_OP (-)
-MARRAYN_SND_OP (*)
-MARRAYN_SND_OP (/)
+MARRAY_SND_OP (+, mx_inline_add)
+MARRAY_SND_OP (-, mx_inline_sub)
+MARRAY_SND_OP (*, mx_inline_mul)
+MARRAY_SND_OP (/, mx_inline_div)
+
+// Element by element MArrayN by MArrayN ops.
 
-#define MARRAY_NDND_OP(FCN, OP) \
-template <class T> \
-MArrayN<T> \
-FCN (const MArrayN<T>& a, const MArrayN<T>& b) \
-{ \
-dim_vector a_dims = a.dims (); \
-dim_vector b_dims = b.dims (); \
-int dims_ok = 1; \
-int any_dims_zero = 0; \
-if (a_dims.length () != b_dims.length ()) \
- dims_ok = 0; \
- else \
-   { \
-     for (int i = 0; i < a_dims.length (); i++) \
-       { \
-	 if (a_dims (i) != b_dims (i)) \
-	   { dims_ok = 0; break; } \
-	 if (a_dims (i) == 0) \
-	   any_dims_zero = 1; \
-       } \
-   } \
- if (!dims_ok) \
-   { \
-     gripe_nonconformant (#FCN, a_dims, b_dims); \
-     return MArrayN<T> (); \
-   } \
- if (any_dims_zero) \
-   return MArrayN<T> (a_dims); \
- octave_idx_type l = a.length (); \
- MArrayN<T> result (a_dims); \
- T* r = result.fortran_vec (); \
- const T *x = a.data (); \
- const T *y = b.data (); \
- DO_VV_OP (r, l, x, OP, y); \
- return result; \
-}
+#define MARRAY_NDND_OP(FCN, OP, FN) \
+  template <class T> \
+  MArrayN<T> \
+  FCN (const MArrayN<T>& a, const MArrayN<T>& b) \
+  { \
+    return do_mm_binary_op<MArrayN<T>, MArrayN<T>, MArrayN<T> > (a, b, FN, #FCN); \
+  }
 
-MARRAY_NDND_OP (operator +, +)
-MARRAY_NDND_OP (operator -, -)
-MARRAY_NDND_OP (product,    *)
-MARRAY_NDND_OP (quotient,   /)
+MARRAY_NDND_OP (operator +, +, mx_inline_add)
+MARRAY_NDND_OP (operator -, -, mx_inline_sub)
+MARRAY_NDND_OP (product,    *, mx_inline_mul)
+MARRAY_NDND_OP (quotient,   /, mx_inline_div)
 
 template <class T>
 MArrayN<T>
@@ -271,12 +191,7 @@
 MArrayN<T>
 operator - (const MArrayN<T>& a)
 {
-  octave_idx_type l = a.length ();
-  MArrayN<T> result (a.dims ());
-  T *r = result.fortran_vec ();
-  const T *x = a.data ();
-  NEG_V (r, l, x);
-  return result;
+  return do_mx_unary_op<MArrayN<T>, MArrayN<T> > (a, mx_inline_uminus); 
 }
 
 /*

--- a/liboctave/MDiagArray2.cc
+++ b/liboctave/MDiagArray2.cc
@@ -47,9 +47,6 @@
   return retval;
 }
 
-// Some functions return a reference to this object after a failure.
-template <class T> MDiagArray2<T> MDiagArray2<T>::nil_array;
-
 // Two dimensional diagonal array with math ops.
 
 // Element by element MDiagArray2 by MDiagArray2 ops.
@@ -58,22 +55,10 @@
 MDiagArray2<T>&
 operator += (MDiagArray2<T>& a, const MDiagArray2<T>& b)
 {
-  octave_idx_type r = a.rows ();
-  octave_idx_type c = a.cols ();
-
-  octave_idx_type b_nr = b.rows ();
-  octave_idx_type b_nc = b.cols ();
-
-  if (r != b_nr || c != b_nc)
-    {
-      gripe_nonconformant ("operator +=", r, c, b_nr, b_nc);
-      return MDiagArray2<T>::nil_array;
-    }
+  if (a.is_shared ())
+    a = a + b;
   else
-    {
-      octave_idx_type l = a.length ();
-      DO_VV_OP2 (T, a, +=, b);
-    }
+    do_mm_inplace_op<MDiagArray2<T>, MDiagArray2<T> > (a, b, mx_inline_add2, "+=");
   return a;
 }
 
@@ -81,42 +66,26 @@
 MDiagArray2<T>&
 operator -= (MDiagArray2<T>& a, const MDiagArray2<T>& b)
 {
-  octave_idx_type r = a.rows ();
-  octave_idx_type c = a.cols ();
-
-  octave_idx_type b_nr = b.rows ();
-  octave_idx_type b_nc = b.cols ();
-
-  if (r != b_nr || c != b_nc)
-    {
-      gripe_nonconformant ("operator -=", r, c, b_nr, b_nc);
-      return MDiagArray2<T>::nil_array;
-    }
+  if (a.is_shared ())
+    a = a - b;
   else
-    {
-      octave_idx_type l = a.length ();
-      DO_VV_OP2 (T, a, -=, b);
-    }
+    do_mm_inplace_op<MDiagArray2<T>, MDiagArray2<T> > (a, b, mx_inline_sub2, "-=");
   return a;
 }
 
+
 // Element by element MDiagArray2 by scalar ops.
 
-#define MARRAY_DAS_OP(OP) \
+#define MARRAY_DAS_OP(OP, FN) \
   template <class T> \
   MDiagArray2<T> \
   operator OP (const MDiagArray2<T>& a, const T& s) \
   { \
-    MDiagArray2<T> result (a.rows (), a.cols ()); \
-    T *r = result.fortran_vec (); \
-    octave_idx_type l = a.length (); \
-    const T *v = a.data (); \
-    DO_VS_OP (r, l, v, OP, s); \
-    return result; \
+    return do_ms_binary_op<MDiagArray2<T>, MDiagArray2<T>, T> (a, s, FN); \
   }
 
-MARRAY_DAS_OP (*)
-MARRAY_DAS_OP (/)
+MARRAY_DAS_OP (*, mx_inline_mul)
+MARRAY_DAS_OP (/, mx_inline_div)
 
 // Element by element scalar by MDiagArray2 ops.
 
@@ -124,44 +93,22 @@
 MDiagArray2<T>
 operator * (const T& s, const MDiagArray2<T>& a)
 {
-  MDiagArray2<T> result (a.rows (), a.cols ()); \
-  T *r = result.fortran_vec (); \
-  octave_idx_type l = a.length (); \
-  const T *v = a.data (); \
-  DO_SV_OP (r, l, s, *, v); \
-  return result; \
+  return do_sm_binary_op<MDiagArray2<T>, T, MDiagArray2<T> > (s, a, mx_inline_mul);
 }
 
 // Element by element MDiagArray2 by MDiagArray2 ops.
 
-#define MARRAY_DADA_OP(FCN, OP) \
+#define MARRAY_DADA_OP(FCN, OP, FN) \
   template <class T> \
   MDiagArray2<T> \
   FCN (const MDiagArray2<T>& a, const MDiagArray2<T>& b) \
   { \
-    octave_idx_type a_nr = a.rows (); \
-    octave_idx_type a_nc = a.cols (); \
-    octave_idx_type b_nr = b.rows (); \
-    octave_idx_type b_nc = b.cols (); \
-    if (a_nr != b_nr || a_nc != b_nc) \
-      { \
-        gripe_nonconformant (#FCN, a_nr, a_nc, b_nr, b_nc); \
-	return MDiagArray2<T> (); \
-      } \
-    if (a_nc == 0 || a_nr == 0) \
-      return MDiagArray2<T> (); \
-    octave_idx_type l = a.length (); \
-    MDiagArray2<T> result (a_nr, a_nc); \
-    T *r = result.fortran_vec (); \
-    const T *x = a.data (); \
-    const T *y = b.data (); \
-    DO_VV_OP (r, l, x, OP, y); \
-    return result; \
+    return do_mm_binary_op<MDiagArray2<T>, MDiagArray2<T>, MDiagArray2<T> > (a, b, FN, #FCN); \
   }
 
-MARRAY_DADA_OP (operator +, +)
-MARRAY_DADA_OP (operator -, -)
-MARRAY_DADA_OP (product,    *)
+MARRAY_DADA_OP (operator +, +, mx_inline_add)
+MARRAY_DADA_OP (operator -, -, mx_inline_sub)
+MARRAY_DADA_OP (product,    *, mx_inline_mul)
 
 // Unary MDiagArray2 ops.
 
@@ -176,12 +123,7 @@
 MDiagArray2<T>
 operator - (const MDiagArray2<T>& a)
 {
-  octave_idx_type l = a.length ();
-  MDiagArray2<T> result (a.rows (), a.cols ());
-  T *r = result.fortran_vec ();
-  const T *x = a.data ();
-  NEG_V (r, l, x);
-  return result;
+  return do_mx_unary_op<MDiagArray2<T>, MDiagArray2<T> > (a, mx_inline_uminus); 
 }
 
 /*

--- a/liboctave/MDiagArray2.h
+++ b/liboctave/MDiagArray2.h
@@ -53,6 +53,8 @@
 
   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) { }
@@ -100,8 +102,6 @@
 
   bool is_multiple_of_identity (T val) const;
 
-  static MDiagArray2<T> nil_array;
-
   // Currently, the OPS functions don't need to be friends, but that
   // may change.
 

--- a/liboctave/boolNDArray.cc
+++ b/liboctave/boolNDArray.cc
@@ -166,20 +166,9 @@
 mx_el_and_assign (boolNDArray& a, const boolNDArray& b)
 {
   if (a.is_shared ())
-    return a = mx_el_and (a, b);
-
-  dim_vector a_dims = a.dims ();
-  dim_vector b_dims = b.dims ();
-
-  if (a_dims != b_dims)
-    gripe_nonconformant ("operator &=", a_dims, b_dims);
+    a = mx_el_and (a, b);
   else
-    {
-      octave_idx_type l = a.length ();
-
-      if (l > 0)
-        DO_VV_OP2 (bool, a, &=, b);
-    }
+    do_mm_inplace_op<boolNDArray, boolNDArray> (a, b, mx_inline_and2, "operator &=");
 
   return a;
 }
@@ -188,20 +177,9 @@
 mx_el_or_assign (boolNDArray& a, const boolNDArray& b)
 {
   if (a.is_shared ())
-    return a = mx_el_and (a, b);
-
-  dim_vector a_dims = a.dims ();
-  dim_vector b_dims = b.dims ();
-
-  if (a_dims != b_dims)
-    gripe_nonconformant ("operator |=", a_dims, b_dims);
+    a = mx_el_or (a, b);
   else
-    {
-      octave_idx_type l = a.length ();
-
-      if (l > 0)
-        DO_VV_OP2 (bool, a, |=, b);
-    }
+    do_mm_inplace_op<boolNDArray, boolNDArray> (a, b, mx_inline_or2, "operator |=");
 
   return a;
 }

--- a/liboctave/dColVector.h
+++ b/liboctave/dColVector.h
@@ -38,6 +38,8 @@
 
   explicit ColumnVector (octave_idx_type n) : MArray<double> (n) { }
 
+  explicit ColumnVector (const dim_vector& dv) : MArray<double> (dv) { }
+
   ColumnVector (octave_idx_type n, double val) : MArray<double> (n, val) { }
 
   ColumnVector (const ColumnVector& a) : MArray<double> (a) { }

--- a/liboctave/dRowVector.h
+++ b/liboctave/dRowVector.h
@@ -38,6 +38,8 @@
 
   explicit RowVector (octave_idx_type n) : MArray<double> (n) { }
 
+  explicit RowVector (const dim_vector& dv) : MArray<double> (dv) { }
+
   RowVector (octave_idx_type n, double val) : MArray<double> (n, val) { }
 
   RowVector (const RowVector& a) : MArray<double> (a) { }

--- a/liboctave/fCColVector.h
+++ b/liboctave/fCColVector.h
@@ -41,6 +41,8 @@
 
   explicit FloatComplexColumnVector (octave_idx_type n) : MArray<FloatComplex> (n) { }
 
+  explicit FloatComplexColumnVector (const dim_vector& dv) : MArray<FloatComplex> (dv) { }
+
   FloatComplexColumnVector (octave_idx_type n, const FloatComplex& val)
     : MArray<FloatComplex> (n, val) { }
 

--- a/liboctave/fCNDArray.cc
+++ b/liboctave/fCNDArray.cc
@@ -1078,16 +1078,18 @@
 FloatComplexNDArray& operator *= (FloatComplexNDArray& a, float s)
 {
   if (a.is_shared ())
-    return a = a * s;
-  DO_VS_OP2 (FloatComplex, a, *=, s)
+    a = a * s;
+  else
+    do_ms_inplace_op<FloatComplexNDArray, float> (a, s, mx_inline_mul2);
   return a;
 }
 
 FloatComplexNDArray& operator /= (FloatComplexNDArray& a, float s)
 {
   if (a.is_shared ())
-    return a = a / s;
-  DO_VS_OP2 (FloatComplex, a, /=, s)
+    a = a / s;
+  else
+    do_ms_inplace_op<FloatComplexNDArray, float> (a, s, mx_inline_div2);
   return a;
 }
 

--- a/liboctave/fCRowVector.h
+++ b/liboctave/fCRowVector.h
@@ -40,6 +40,8 @@
 
   explicit FloatComplexRowVector (octave_idx_type n) : MArray<FloatComplex> (n) { }
 
+  explicit FloatComplexRowVector (const dim_vector& dv) : MArray<FloatComplex> (dv) { }
+
   FloatComplexRowVector (octave_idx_type n, const FloatComplex& val) : MArray<FloatComplex> (n, val) { }
 
   FloatComplexRowVector (const FloatComplexRowVector& a) : MArray<FloatComplex> (a) { }

--- a/liboctave/fColVector.h
+++ b/liboctave/fColVector.h
@@ -38,6 +38,8 @@
 
   explicit FloatColumnVector (octave_idx_type n) : MArray<float> (n) { }
 
+  explicit FloatColumnVector (const dim_vector& dv) : MArray<float> (dv) { }
+
   FloatColumnVector (octave_idx_type n, float val) : MArray<float> (n, val) { }
 
   FloatColumnVector (const FloatColumnVector& a) : MArray<float> (a) { }

--- a/liboctave/fRowVector.h
+++ b/liboctave/fRowVector.h
@@ -38,6 +38,8 @@
 
   explicit FloatRowVector (octave_idx_type n) : MArray<float> (n) { }
 
+  explicit FloatRowVector (const dim_vector& dv) : MArray<float> (dv) { }
+
   FloatRowVector (octave_idx_type n, float val) : MArray<float> (n, val) { }
 
   FloatRowVector (const FloatRowVector& a) : MArray<float> (a) { }

--- a/liboctave/mx-inlines.cc
+++ b/liboctave/mx-inlines.cc
@@ -168,6 +168,17 @@
 DEFMXBOOLOP (mx_inline_and_not, , &, !)
 DEFMXBOOLOP (mx_inline_or_not, , |, !)
 
+#define DEFMXBOOLOPEQ(F, OP) \
+template <class X> \
+inline void F (size_t n, bool *r, const X *x) \
+{ \
+  for (size_t i = 0; i < n; i++) \
+    r[i] OP logical_value (x[i]); \
+} \
+
+DEFMXBOOLOPEQ (mx_inline_and2, &=)
+DEFMXBOOLOPEQ (mx_inline_or2, |=)
+
 template <class T> 
 inline bool 
 mx_inline_any_nan (size_t, const T*) { return false; }
@@ -230,7 +241,7 @@
                             const typename XNDA::element_type *))
 {
   RNDA r (x.dims ());
-  op (r.nelem (), r.fortran_vec (), x.data ());
+  op (r.length (), r.fortran_vec (), x.data ());
   return r;
 }
 
@@ -256,7 +267,7 @@
   if (dx == dy)
     {
       RNDA r (dx);
-      op (r.nelem (), r.fortran_vec (), x.data (), y.data ());
+      op (r.length (), r.fortran_vec (), x.data (), y.data ());
       return r;
     }
   else
@@ -273,7 +284,7 @@
                              const typename XNDA::element_type *, YS))
 {
   RNDA r (x.dims ());
-  op (r.nelem (), r.fortran_vec (), x.data (), y);
+  op (r.length (), r.fortran_vec (), x.data (), y);
   return r;
 }
 
@@ -284,7 +295,7 @@
                              const typename YNDA::element_type *))
 {
   RNDA r (y.dims ());
-  op (r.nelem (), r.fortran_vec (), x, y.data ());
+  op (r.length (), r.fortran_vec (), x, y.data ());
   return r;
 }
 
@@ -297,15 +308,10 @@
 {
   dim_vector dr = r.dims (), dx = x.dims ();
   if (dr == dx)
-    {
-      op (r.nelem (), r.fortran_vec (), x.data ());
-      return r;
-    }
+    op (r.length (), r.fortran_vec (), x.data ());
   else
-    {
-      gripe_nonconformant (opname, dr, dx);
-      return RNDA ();
-    }
+    gripe_nonconformant (opname, dr, dx);
+  return r;
 }
 
 template <class RNDA, class XS>
@@ -313,7 +319,7 @@
 do_ms_inplace_op (RNDA& r, const XS& x,
                   void (*op) (size_t, typename RNDA::element_type *, XS))
 {
-  op (r.nelem (), r.fortran_vec (), x);
+  op (r.length (), r.fortran_vec (), x);
   return r;
 }
 

--- a/src/pr-output.cc
+++ b/src/pr-output.cc
@@ -852,7 +852,7 @@
   double r_abs = rp < 0.0 ? -rp : rp;
   double i_abs = ip < 0.0 ? -ip : ip;
 
-  int r_x = (xisinf (rp) || xisnan (rp) || xr_abs == 0.0)
+  int r_x = (xisinf (rp) || xisnan (rp) || r_abs == 0.0)
     ? 0 : static_cast<int> (floor (log10 (r_abs) + 1.0));
 
   int i_x = (xisinf (ip) || xisnan (ip) || i_abs == 0.0)