--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -38,6 +38,7 @@
 #include <sys/types.h>
 #endif
 
+#include "CMatrix.h"
 #include "CmplxAEPBAL.h"
 #include "CmplxDET.h"
 #include "CmplxSCHUR.h"
@@ -48,6 +49,8 @@
 #include "lo-mappers.h"
 #include "lo-utils.h"
 #include "mx-base.h"
+#include "mx-cm-dm.h"
+#include "mx-cm-s.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
 
@@ -159,6 +162,13 @@
 // XXX FIXME XXX -- could we use a templated mixed-type copy function
 // here?
 
+ComplexMatrix::ComplexMatrix (const boolMatrix& a)
+{
+  for (int i = 0; i < a.cols (); i++)
+    for (int j = 0; j < a.rows (); j++)
+      elem (i, j) = a.elem (i, j);
+}
+
 ComplexMatrix::ComplexMatrix (const charMatrix& a)
 {
   for (int i = 0; i < a.cols (); i++)
@@ -1712,94 +1722,6 @@
   return retval;
 }
 
-// diagonal matrix by scalar -> matrix operations
-
-ComplexMatrix
-operator + (const DiagMatrix& a, const Complex& s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator - (const DiagMatrix& a, const Complex& s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), -s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator + (const ComplexDiagMatrix& a, double s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator - (const ComplexDiagMatrix& a, double s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), -s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator + (const ComplexDiagMatrix& a, const Complex& s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator - (const ComplexDiagMatrix& a, const Complex& s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), -s);
-  return a + tmp;
-}
-
-// scalar by diagonal matrix -> matrix operations
-
-ComplexMatrix
-operator + (const Complex& s, const DiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp + a;
-}
-
-ComplexMatrix
-operator - (const Complex& s, const DiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp - a;
-}
-
-ComplexMatrix
-operator + (double s, const ComplexDiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp + a;
-}
-
-ComplexMatrix
-operator - (double s, const ComplexDiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp - a;
-}
-
-ComplexMatrix
-operator + (const Complex& s, const ComplexDiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp + a;
-}
-
-ComplexMatrix
-operator - (const Complex& s, const ComplexDiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp - a;
-}
-
 // matrix by diagonal matrix -> matrix operations
 
 ComplexMatrix&
@@ -1886,411 +1808,6 @@
   return *this;
 }
 
-ComplexMatrix
-operator + (const Matrix& m, const ComplexDiagMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (m);
-  for (int i = 0; i < a.length (); i++)
-    result.elem (i, i) += a.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator - (const Matrix& m, const ComplexDiagMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (m);
-  for (int i = 0; i < a.length (); i++)
-    result.elem (i, i) -= a.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator * (const Matrix& m, const ComplexDiagMatrix& a)
-{
-  ComplexMatrix retval;
-
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-  else
-    {
-      if (nr == 0 || nc == 0 || a_nc == 0)
-	retval.resize (nr, a_nc, 0.0);
-      else
-	{
-	  retval.resize (nr, a_nc);
-	  Complex *c = retval.fortran_vec ();
-
-	  Complex *ctmp = 0;
-
-	  for (int j = 0; j < a.length (); j++)
-	    {
-	      int idx = j * nr;
-	      ctmp = c + idx;
-	      if (a.elem (j, j) == 1.0)
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = m.elem (i, j);
-		}
-	      else if (a.elem (j, j) == 0.0)
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = 0.0;
-		}
-	      else
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
-		}
-	    }
-
-	  if (a_nr < a_nc)
-	    {
-	      for (int i = nr * nc; i < nr * a_nc; i++)
-		ctmp[i] = 0.0;
-	    }
-	}
-    }
-
-  return retval;
-}
-
-// diagonal matrix by matrix -> matrix operations
-
-ComplexMatrix
-operator + (const DiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator - (const DiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (-a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator * (const DiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    {
-      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return ComplexMatrix (nr, nc, 0.0);
-
-  ComplexMatrix c (nr, a_nc);
-
-  for (int i = 0; i < m.length (); i++)
-    {
-      if (m.elem (i, i) == 1.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = a.elem (i, j);
-	}
-      else if (m.elem (i, i) == 0.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = 0.0;
-	}
-      else
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
-	}
-    }
-
-  if (nr > nc)
-    {
-      for (int j = 0; j < a_nc; j++)
-	for (int i = a_nr; i < nr; i++)
-	  c.elem (i, j) = 0.0;
-    }
-
-  return c;
-}
-
-ComplexMatrix
-operator + (const ComplexDiagMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator - (const ComplexDiagMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (-a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator * (const ComplexDiagMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    {
-      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return ComplexMatrix (nr, a_nc, 0.0);
-
-  ComplexMatrix c (nr, a_nc);
-
-  for (int i = 0; i < m.length (); i++)
-    {
-      if (m.elem (i, i) == 1.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = a.elem (i, j);
-	}
-      else if (m.elem (i, i) == 0.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = 0.0;
-	}
-      else
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
-	}
-    }
-
-  if (nr > nc)
-    {
-      for (int j = 0; j < a_nc; j++)
-	for (int i = a_nr; i < nr; i++)
-	  c.elem (i, j) = 0.0;
-    }
-
-  return c;
-}
-
-ComplexMatrix
-operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (-a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    {
-      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return ComplexMatrix (nr, a_nc, 0.0);
-
-  ComplexMatrix c (nr, a_nc);
-
-  for (int i = 0; i < m.length (); i++)
-    {
-      if (m.elem (i, i) == 1.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = a.elem (i, j);
-	}
-      else if (m.elem (i, i) == 0.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = 0.0;
-	}
-      else
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
-	}
-    }
-
-  if (nr > nc)
-    {
-      for (int j = 0; j < a_nc; j++)
-	for (int i = a_nr; i < nr; i++)
-	  c.elem (i, j) = 0.0;
-    }
-
-  return c;
-}
-
 // matrix by matrix -> matrix operations
 
 ComplexMatrix&
@@ -2397,553 +1914,6 @@
   return Matrix (not (data (), length ()), rows (), cols ());
 }
 
-// matrix by scalar -> matrix operations
-
-ComplexMatrix
-operator + (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (add (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator - (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (subtract (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator * (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (multiply (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator / (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (divide (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator + (const ComplexMatrix& a, double s)
-{
-  return ComplexMatrix (add (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator - (const ComplexMatrix& a, double s)
-{
-  return ComplexMatrix (subtract (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator * (const ComplexMatrix& a, double s)
-{
-  return ComplexMatrix (multiply (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator / (const ComplexMatrix& a, double s)
-{
-  return ComplexMatrix (divide (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-// scalar by matrix -> matrix operations
-
-ComplexMatrix
-operator + (double s, const ComplexMatrix& a)
-{
-  return ComplexMatrix (add (a.data (), a.length (), s), a.rows (),
-			a.cols ());
-}
-
-ComplexMatrix
-operator - (double s, const ComplexMatrix& a)
-{
-  return ComplexMatrix (subtract (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator * (double s, const ComplexMatrix& a)
-{
-  return ComplexMatrix (multiply (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator / (double s, const ComplexMatrix& a)
-{
-  return ComplexMatrix (divide (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator + (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (add (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator - (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (subtract (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator * (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (multiply (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator / (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (divide (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-// matrix by diagonal matrix -> matrix operations
-
-ComplexMatrix
-operator + (const ComplexMatrix& m, const DiagMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (m);
-  for (int i = 0; i < a.length (); i++)
-    result.elem (i, i) += a.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator - (const ComplexMatrix& m, const DiagMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (m);
-  for (int i = 0; i < a.length (); i++)
-    result.elem (i, i) -= a.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator * (const ComplexMatrix& m, const DiagMatrix& a)
-{
-  ComplexMatrix retval;
-
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-  else
-    {
-      if (nr == 0 || nc == 0 || a_nc == 0)
-	retval.resize (nr, nc, 0.0);
-      else
-	{
-	  retval.resize (nr, a_nc);
-	  Complex *c = retval.fortran_vec ();
-	  Complex *ctmp = 0;
-
-	  for (int j = 0; j < a.length (); j++)
-	    {
-	      int idx = j * nr;
-	      ctmp = c + idx;
-	      if (a.elem (j, j) == 1.0)
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = m.elem (i, j);
-		}
-	      else if (a.elem (j, j) == 0.0)
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = 0.0;
-		}
-	      else
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
-		}
-	    }
-
-	  if (a.rows () < a_nc)
-	    {
-	      for (int i = nr * nc; i < nr * a_nc; i++)
-		ctmp[i] = 0.0;
-	    }
-	}
-    }
-
-  return retval;
-}
-
-ComplexMatrix
-operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (m);
-  for (int i = 0; i < a.length (); i++)
-    result.elem (i, i) += a.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator - (const ComplexMatrix& m, const ComplexDiagMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  ComplexMatrix result (m);
-  for (int i = 0; i < a.length (); i++)
-    result.elem (i, i) -= a.elem (i, i);
-
-  return result;
-}
-
-ComplexMatrix
-operator * (const ComplexMatrix& m, const ComplexDiagMatrix& a)
-{
-  ComplexMatrix retval;
-
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-  else
-    {
-      if (nr == 0 || nc == 0 || a_nc == 0)
-	retval.resize (nr, nc, 0.0);
-      else
-	{
-	  retval.resize (nr, nc);
-	  Complex *c = retval.fortran_vec ();
-	  Complex *ctmp = 0;
-
-	  for (int j = 0; j < a.length (); j++)
-	    {
-	      int idx = j * nr;
-	      ctmp = c + idx;
-	      if (a.elem (j, j) == 1.0)
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = m.elem (i, j);
-		}
-	      else if (a.elem (j, j) == 0.0)
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = 0.0;
-		}
-	      else
-		{
-		  for (int i = 0; i < nr; i++)
-		    ctmp[i] = a.elem (j, j) * m.elem (i, j);
-		}
-	    }
-
-	  if (a.rows () < a_nc)
-	    {
-	      for (int i = nr * nc; i < nr * a_nc; i++)
-		ctmp[i] = 0.0;
-	    }
-	}
-    }
-
-  return retval;
-}
-
-// matrix by matrix -> matrix operations
-
-ComplexMatrix
-operator + (const ComplexMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-operator - (const ComplexMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-operator + (const Matrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-operator - (const Matrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-operator * (const ComplexMatrix& m, const Matrix& a)
-{
-  ComplexMatrix tmp (a);
-  return m * tmp;
-}
-
-ComplexMatrix
-operator * (const Matrix& m, const ComplexMatrix& a)
-{
-  ComplexMatrix tmp (m);
-  return tmp * a;
-}
-
-ComplexMatrix
-operator * (const ComplexMatrix& m, const ComplexMatrix& a)
-{
-  ComplexMatrix retval;
-
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nc != a_nr)
-    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
-  else
-    {
-      if (nr == 0 || nc == 0 || a_nc == 0)
-	retval.resize (nr, nc, 0.0);
-      else
-	{
-	  int ld  = nr;
-	  int lda = a.rows ();
-
-	  retval.resize (nr, a_nc);
-	  Complex *c = retval.fortran_vec ();
-
-	  F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0,
-				   m.data (), ld, a.data (), lda, 0.0,
-				   c, nr, 1L, 1L));
-
-	  if (f77_exception_encountered)
-	    (*current_liboctave_error_handler)
-	      ("unrecoverable error in zgemm");
-	}
-    }
-
-  return retval;
-}
-
-ComplexMatrix
-product (const ComplexMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("product", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-quotient (const ComplexMatrix& m, const Matrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("quotient", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-product (const Matrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("product", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-quotient (const Matrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-
-  if (nr != a_nr || nc != a_nc)
-    {
-      gripe_nonconformant ("quotient", nr, nc, a_nr, a_nc);
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ComplexMatrix (nr, nc);
-
-  return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc);
-}
-
 // other operations
 
 ComplexMatrix
@@ -3856,6 +2826,58 @@
   return retval;
 }
 
+ComplexMatrix
+operator * (const ComplexMatrix& m, const Matrix& a)
+{
+  ComplexMatrix tmp (a);
+  return m * tmp;
+}
+
+ComplexMatrix
+operator * (const Matrix& m, const ComplexMatrix& a)
+{
+  ComplexMatrix tmp (m);
+  return tmp * a;
+}
+
+ComplexMatrix
+operator * (const ComplexMatrix& m, const ComplexMatrix& a)
+{
+  ComplexMatrix retval;
+
+  int nr = m.rows ();
+  int nc = m.cols ();
+
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+
+  if (nc != a_nr)
+    gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
+  else
+    {
+      if (nr == 0 || nc == 0 || a_nc == 0)
+	retval.resize (nr, nc, 0.0);
+      else
+	{
+	  int ld  = nr;
+	  int lda = a.rows ();
+
+	  retval.resize (nr, a_nc);
+	  Complex *c = retval.fortran_vec ();
+
+	  F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0,
+				   m.data (), ld, a.data (), lda, 0.0,
+				   c, nr, 1L, 1L));
+
+	  if (f77_exception_encountered)
+	    (*current_liboctave_error_handler)
+	      ("unrecoverable error in zgemm");
+	}
+    }
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***