new file mode 100644
--- /dev/null
+++ b/liboctave/CColVector.cc
@@ -0,0 +1,659 @@
+// ColumnVector manipulations.                            -*- C++ -*-
+/*
+
+Copyright (C) 1992, 1993, 1994 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#include <iostream.h>
+
+#include <Complex.h>
+
+#include "mx-base.h"
+#include "mx-inlines.cc"
+#include "f77-uscore.h"
+#include "lo-error.h"
+
+// Fortran functions we call.
+
+extern "C"
+{
+  int F77_FCN (zgemm) (const char*, const char*, const int*,
+		       const int*, const int*, const Complex*,
+		       const Complex*, const int*, const Complex*,
+		       const int*, const Complex*, Complex*, const int*,
+		       long, long);
+}
+
+/*
+ * Complex Column Vector class
+ */
+
+#define KLUDGE_VECTORS
+#define TYPE Complex
+#define KL_VEC_TYPE ComplexColumnVector
+#include "mx-kludge.cc"
+#undef KLUDGE_VECTORS
+#undef TYPE
+#undef KL_VEC_TYPE
+
+ComplexColumnVector::ComplexColumnVector (const ColumnVector& a)
+   : Array<Complex> (a.length ())
+{
+  for (int i = 0; i < length (); i++)
+    elem (i) = a.elem (i);
+}
+
+#if 0
+ComplexColumnVector&
+ComplexColumnVector::resize (int n)
+{
+  if (n < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimension");
+      return *this;
+    }
+
+  Complex *new_data = (Complex *) NULL;
+  if (n > 0)
+    {
+      new_data = new Complex [n];
+      int min_len = len < n ? len : n;
+
+      for (int i = 0; i < min_len; i++)
+	new_data[i] = data[i];
+    }
+
+  delete [] data;
+  len = n;
+  data = new_data;
+
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::resize (int n, double val)
+{
+  int old_len = len;
+  resize (n);
+  for (int i = old_len; i < len; i++)
+    data[i] = val;
+
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::resize (int n, const Complex& val)
+{
+  int old_len = len;
+  resize (n);
+  for (int i = old_len; i < len; i++)
+    data[i] = val;
+
+  return *this;
+}
+#endif
+
+int
+ComplexColumnVector::operator == (const ComplexColumnVector& a) const
+{
+  int len = length ();
+  if (len != a.length ())
+    return 0;
+  return equal (data (), a.data (), len);
+}
+
+int
+ComplexColumnVector::operator != (const ComplexColumnVector& a) const
+{
+  return !(*this == a);
+}
+
+// destructive insert/delete/reorder operations
+
+ComplexColumnVector&
+ComplexColumnVector::insert (const ColumnVector& a, int r)
+{
+  int a_len = a.length ();
+  if (r < 0 || r + a_len - 1 > length ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (r+i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::insert (const ComplexColumnVector& a, int r)
+{
+  int a_len = a.length ();
+  if (r < 0 || r + a_len - 1 > length ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (r+i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::fill (double val)
+{
+  int len = length ();
+  if (len > 0)
+    for (int i = 0; i < len; i++)
+      elem (i) = val;
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::fill (const Complex& val)
+{
+  int len = length ();
+  if (len > 0)
+    for (int i = 0; i < len; i++)
+      elem (i) = val;
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::fill (double val, int r1, int r2)
+{
+  int len = length ();
+  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+
+  for (int i = r1; i <= r2; i++)
+    elem (i) = val;
+
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::fill (const Complex& val, int r1, int r2)
+{
+  int len = length ();
+  if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+
+  for (int i = r1; i <= r2; i++)
+    elem (i) = val;
+
+  return *this;
+}
+
+ComplexColumnVector
+ComplexColumnVector::stack (const ColumnVector& a) const
+{
+  int len = length ();
+  int nr_insert = len;
+  ComplexColumnVector retval (len + a.length ());
+  retval.insert (*this, 0);
+  retval.insert (a, nr_insert);
+  return retval;
+}
+
+ComplexColumnVector
+ComplexColumnVector::stack (const ComplexColumnVector& a) const
+{
+  int len = length ();
+  int nr_insert = len;
+  ComplexColumnVector retval (len + a.length ());
+  retval.insert (*this, 0);
+  retval.insert (a, nr_insert);
+  return retval;
+}
+
+ComplexRowVector
+ComplexColumnVector::hermitian (void) const
+{
+  int len = length ();
+  return ComplexRowVector (conj_dup (data (), len), len);
+}
+
+ComplexRowVector
+ComplexColumnVector::transpose (void) const
+{
+  int len = length ();
+  return ComplexRowVector (dup (data (), len), len);
+}
+
+ColumnVector
+real (const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  ColumnVector retval;
+  if (a_len > 0)
+    retval = ColumnVector (real_dup (a.data (), a_len), a_len);
+  return retval;
+}
+
+ColumnVector
+imag (const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  ColumnVector retval;
+  if (a_len > 0)
+    retval = ColumnVector (imag_dup (a.data (), a_len), a_len);
+  return retval;
+}
+
+ComplexColumnVector
+conj (const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  ComplexColumnVector retval;
+  if (a_len > 0)
+    retval = ComplexColumnVector (conj_dup (a.data (), a_len), a_len);
+  return retval;
+}
+
+// resize is the destructive equivalent for this one
+
+ComplexColumnVector
+ComplexColumnVector::extract (int r1, int r2) const
+{
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+
+  int new_r = r2 - r1 + 1;
+
+  ComplexColumnVector result (new_r);
+
+  for (int i = 0; i < new_r; i++)
+    result.elem (i) = elem (r1+i);
+
+  return result;
+}
+
+// column vector by column vector -> column vector operations
+
+ComplexColumnVector&
+ComplexColumnVector::operator += (const ColumnVector& a)
+{
+  int len = length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector += operation attempted");
+      return *this;
+    }
+
+  if (len == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  add2 (d, a.data (), len);
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::operator -= (const ColumnVector& a)
+{
+  int len = length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector -= operation attempted");
+      return *this;
+    }
+
+  if (len == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  subtract2 (d, a.data (), len);
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::operator += (const ComplexColumnVector& a)
+{
+  int len = length ();
+
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector += operation attempted");
+      return *this;
+    }
+
+  if (len == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  add2 (d, a.data (), len);
+  return *this;
+}
+
+ComplexColumnVector&
+ComplexColumnVector::operator -= (const ComplexColumnVector& a)
+{
+  int len = length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector -= operation attempted");
+      return *this;
+    }
+
+  if (len == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  subtract2 (d, a.data (), len);
+  return *this;
+}
+
+// column vector by scalar -> column vector operations
+
+ComplexColumnVector
+operator + (const ComplexColumnVector& v, double s)
+{
+  int len = v.length ();
+  return ComplexColumnVector (add (v.data (), len, s), len);
+}
+
+ComplexColumnVector
+operator - (const ComplexColumnVector& v, double s)
+{
+  int len = v.length ();
+  return ComplexColumnVector (subtract (v.data (), len, s), len);
+}
+
+ComplexColumnVector
+operator * (const ComplexColumnVector& v, double s)
+{
+  int len = v.length ();
+  return ComplexColumnVector (multiply (v.data (), len, s), len);
+}
+
+ComplexColumnVector
+operator / (const ComplexColumnVector& v, double s)
+{
+  int len = v.length ();
+  return ComplexColumnVector (divide (v.data (), len, s), len);
+}
+
+// scalar by column vector -> column vector operations
+
+ComplexColumnVector
+operator + (double s, const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  return ComplexColumnVector (add (a.data (), a_len, s), a_len);
+}
+
+ComplexColumnVector
+operator - (double s, const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  return ComplexColumnVector (subtract (s, a.data (), a_len), a_len);
+}
+
+ComplexColumnVector
+operator * (double s, const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  return ComplexColumnVector (multiply (a.data (), a_len, s), a_len);
+}
+
+ComplexColumnVector
+operator / (double s, const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  return ComplexColumnVector (divide (s, a.data (), a_len), a_len);
+}
+
+// column vector by row vector -> matrix operations
+
+ComplexMatrix
+operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
+{
+  int len = v.length ();
+  int a_len = a.length ();
+  if (len != a_len)
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector multiplication attempted");
+      return ComplexMatrix ();
+    }
+
+  if (len == 0)
+    return ComplexMatrix (len, len, 0.0);
+
+  char transa = 'N';
+  char transb = 'N';
+  Complex alpha (1.0);
+  Complex beta (0.0);
+  int anr = 1;
+
+  Complex *c = new Complex [len * a_len];
+
+  F77_FCN (zgemm) (&transa, &transb, &len, &a_len, &anr, &alpha,
+		   v.data (), &len, a.data (), &anr, &beta, c, &len,
+		   1L, 1L);
+
+  return ComplexMatrix (c, len, a_len);
+}
+
+// column vector by column vector -> column vector operations
+
+ComplexColumnVector
+operator + (const ComplexColumnVector& v, const ColumnVector& a)
+{
+  int len = v.length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector addition attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (len == 0)
+    return ComplexColumnVector (0);
+
+  return ComplexColumnVector (add (v.data (), a.data (), len), len);
+}
+
+ComplexColumnVector
+operator - (const ComplexColumnVector& v, const ColumnVector& a)
+{
+  int len = v.length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector subtraction attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (len == 0)
+    return ComplexColumnVector (0);
+
+  return ComplexColumnVector (subtract (v.data (), a.data (), len), len);
+}
+
+ComplexColumnVector
+product (const ComplexColumnVector& v, const ColumnVector& a)
+{
+  int len = v.length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector product attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (len == 0)
+    return ComplexColumnVector (0);
+
+  return ComplexColumnVector (multiply (v.data (), a.data (), len), len);
+}
+
+ComplexColumnVector
+quotient (const ComplexColumnVector& v, const ColumnVector& a)
+{
+  int len = v.length ();
+  if (len != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant vector quotient attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (len == 0)
+    return ComplexColumnVector (0);
+
+  return ComplexColumnVector (divide (v.data (), a.data (), len), len);
+}
+
+// other operations
+
+ComplexColumnVector
+map (c_c_Mapper f, const ComplexColumnVector& a)
+{
+  ComplexColumnVector b (a);
+  b.map (f);
+  return b;
+}
+
+ColumnVector
+map (d_c_Mapper f, const ComplexColumnVector& a)
+{
+  int a_len = a.length ();
+  ColumnVector b (a_len);
+  for (int i = 0; i < a_len; i++)
+    b.elem (i) = f (a.elem (i));
+  return b;
+}
+
+void
+ComplexColumnVector::map (c_c_Mapper f)
+{
+  for (int i = 0; i < length (); i++)
+    elem (i) = f (elem (i));
+}
+
+Complex
+ComplexColumnVector::min (void) const
+{
+  int len = length ();
+  if (len == 0)
+    return 0.0;
+
+  Complex res = elem (0);
+  double absres = abs (res);
+
+  for (int i = 1; i < len; i++)
+    if (abs (elem (i)) < absres)
+      {
+	res = elem (i);
+	absres = abs (res);
+      }
+
+  return res;
+}
+
+Complex
+ComplexColumnVector::max (void) const
+{
+  int len = length ();
+  if (len == 0)
+    return 0.0;
+
+  Complex res = elem (0);
+  double absres = abs (res);
+
+  for (int i = 1; i < len; i++)
+    if (abs (elem (i)) > absres)
+      {
+	res = elem (i);
+	absres = abs (res);
+      }
+
+  return res;
+}
+
+// i/o
+
+ostream&
+operator << (ostream& os, const ComplexColumnVector& a)
+{
+//  int field_width = os.precision () + 7;
+  for (int i = 0; i < a.length (); i++)
+    os << /* setw (field_width) << */ a.elem (i) << "\n";
+  return os;
+}
+
+istream&
+operator >> (istream& is, ComplexColumnVector& a)
+{
+  int len = a.length();
+
+  if (len < 1)
+    is.clear (ios::badbit);
+  else
+    {
+      double tmp;
+      for (int i = 0; i < len; i++)
+        {
+          is >> tmp;
+          if (is)
+            a.elem (i) = tmp;
+          else
+            break;
+        }
+    }
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/CColVector.h
@@ -0,0 +1,171 @@
+//                                  -*- C++ -*-
+/*
+
+Copyright (C) 1992, 1993, 1994 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#if !defined (octave_ComplexColumnVector_h)
+#define octave_ComplexColumnVector_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include "Array.h"
+
+#include "mx-defs.h"
+
+extern "C++" {
+
+class ComplexColumnVector : public Array<Complex>
+{
+friend class ComplexRowVector;
+
+public:
+
+  ComplexColumnVector (void) : Array<Complex> () { }
+  ComplexColumnVector (int n) : Array<Complex> (n) { }
+  ComplexColumnVector (int n, const Complex& val)
+    : Array<Complex> (n, val) { }
+  ComplexColumnVector (const ColumnVector& a);
+  ComplexColumnVector (const Array<Complex>& a) : Array<Complex> (a) { }
+  ComplexColumnVector (const ComplexColumnVector& a) : Array<Complex> (a) { }
+//  ComplexColumnVector (double a) : Array<Complex> (1, a) { }
+//  ComplexColumnVector (const Complex& a) : Array<Complex> (1, a) { }
+
+  ComplexColumnVector& operator = (const ComplexColumnVector& a)
+    {
+      Array<Complex>::operator = (a);
+      return *this;
+    }
+
+//  operator Array<Complex>& () const { return *this; }
+
+  int operator == (const ComplexColumnVector& a) const;
+  int operator != (const ComplexColumnVector& a) const;
+
+// destructive insert/delete/reorder operations
+
+  ComplexColumnVector& insert (const ColumnVector& a, int r);
+  ComplexColumnVector& insert (const ComplexColumnVector& a, int r);
+
+  ComplexColumnVector& fill (double val);
+  ComplexColumnVector& fill (const Complex& val);
+  ComplexColumnVector& fill (double val, int r1, int r2);
+  ComplexColumnVector& fill (const Complex& val, int r1, int r2);
+
+  ComplexColumnVector stack (const ColumnVector& a) const;
+  ComplexColumnVector stack (const ComplexColumnVector& a) const;
+
+  ComplexRowVector hermitian (void) const;  // complex conjugate transpose.
+  ComplexRowVector transpose (void) const;
+
+  friend ColumnVector real (const ComplexColumnVector& a);
+  friend ColumnVector imag (const ComplexColumnVector& a);
+  friend ComplexColumnVector conj (const ComplexColumnVector& a);
+
+// resize is the destructive equivalent for this one
+
+  ComplexColumnVector extract (int r1, int r2) const;
+
+// column vector by column vector -> column vector operations
+
+  ComplexColumnVector& operator += (const ColumnVector& a);
+  ComplexColumnVector& operator -= (const ColumnVector& a);
+
+  ComplexColumnVector& operator += (const ComplexColumnVector& a);
+  ComplexColumnVector& operator -= (const ComplexColumnVector& a);
+
+// column vector by scalar -> column vector operations
+
+  friend ComplexColumnVector operator + (const ComplexColumnVector& a,
+					 double s);
+  friend ComplexColumnVector operator - (const ComplexColumnVector& a,
+					 double s);
+  friend ComplexColumnVector operator * (const ComplexColumnVector& a,
+					 double s);
+  friend ComplexColumnVector operator / (const ComplexColumnVector& a,
+					 double s);
+
+// scalar by column vector -> column vector operations
+
+  friend ComplexColumnVector operator + (double s,
+					 const ComplexColumnVector& a); 
+  friend ComplexColumnVector operator - (double s,
+					 const ComplexColumnVector& a);
+  friend ComplexColumnVector operator * (double s,
+					 const ComplexColumnVector& a);
+  friend ComplexColumnVector operator / (double s,
+					 const ComplexColumnVector& a);
+
+// column vector by row vector -> matrix operations
+
+  friend ComplexMatrix operator * (const ComplexColumnVector& a,
+				   const ComplexRowVector& b);
+
+// column vector by column vector -> column vector operations
+
+  friend ComplexColumnVector operator + (const ComplexColumnVector& a,
+					 const ColumnVector& b);
+  friend ComplexColumnVector operator - (const ComplexColumnVector& a,
+					 const ColumnVector& b);
+
+  friend ComplexColumnVector product (const ComplexColumnVector& a,
+				      const ColumnVector& b);
+  friend ComplexColumnVector quotient (const ComplexColumnVector& a,
+				       const ColumnVector& b);
+
+// other operations
+
+  friend ComplexColumnVector map (c_c_Mapper f, const ComplexColumnVector& a);
+  friend ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a);
+  void map (c_c_Mapper f);
+
+  Complex min (void) const;
+  Complex max (void) const;
+
+// i/o
+
+  friend ostream& operator << (ostream& os, const ComplexColumnVector& a);
+  friend ostream& operator >> (ostream& is, ComplexColumnVector& a);
+
+#define KLUDGE_VECTORS
+#define TYPE Complex
+#define KL_VEC_TYPE ComplexColumnVector
+#include "mx-kludge.h"
+#undef KLUDGE_VECTORS
+#undef TYPE
+#undef KL_VEC_TYPE
+
+private:
+
+  ComplexColumnVector (Complex *d, int l) : Array<Complex> (d, l) { }
+};
+
+} // extern "C++"
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/CDiagMatrix.cc
@@ -0,0 +1,1168 @@
+// DiagMatrix manipulations.                             -*- C++ -*-
+/*
+
+Copyright (C) 1992, 1993, 1994 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#include <iostream.h>
+
+#include <Complex.h>
+
+#include "mx-base.h"
+#include "mx-inlines.cc"
+#include "lo-error.h"
+
+/*
+ * Complex Diagonal Matrix class
+ */
+
+#define KLUDGE_DIAG_MATRICES
+#define TYPE Complex
+#define KL_DMAT_TYPE ComplexDiagMatrix
+#include "mx-kludge.cc"
+#undef KLUDGE_DIAG_MATRICES
+#undef TYPE
+#undef KL_DMAT_TYPE
+
+ComplexDiagMatrix::ComplexDiagMatrix (const RowVector& a)
+  : DiagArray<Complex> (a.length ())
+{
+  for (int i = 0; i < length (); i++)
+    elem (i, i) = a.elem (i);
+}
+
+ComplexDiagMatrix::ComplexDiagMatrix (const ColumnVector& a)
+  : DiagArray<Complex> (a.length ())
+{
+  for (int i = 0; i < length (); i++)
+    elem (i, i) = a.elem (i);
+}
+
+ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a)
+  : DiagArray<Complex> (a.rows (), a.cols ())
+{
+  for (int i = 0; i < length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+#if 0
+ComplexDiagMatrix&
+ComplexDiagMatrix::resize (int r, int c)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimensions");
+      return *this;
+    }
+
+  int new_len = r < c ? r : c;
+  Complex *new_data = (Complex *) NULL;
+  if (new_len > 0)
+    {
+      new_data = new Complex [new_len];
+
+      int min_len = new_len < len ? new_len : len;
+
+      for (int i = 0; i < min_len; i++)
+	new_data[i] = data[i];
+    }
+
+  delete [] data;
+  nr = r;
+  nc = c;
+  len = new_len;
+  data = new_data;
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::resize (int r, int c, double val)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimensions");
+      return *this;
+    }
+
+  int new_len = r < c ? r : c;
+  Complex *new_data = (Complex *) NULL;
+  if (new_len > 0)
+    {
+      new_data = new Complex [new_len];
+
+      int min_len = new_len < len ? new_len : len;
+
+      for (int i = 0; i < min_len; i++)
+	new_data[i] = data[i];
+
+      for (i = min_len; i < new_len; i++)
+	new_data[i] = val;
+    }
+
+  delete [] data;
+  nr = r;
+  nc = c;
+  len = new_len;
+  data = new_data;
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::resize (int r, int c, const Complex& val)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimensions");
+      return *this;
+    }
+
+  int new_len = r < c ? r : c;
+  Complex *new_data = (Complex *) NULL;
+  if (new_len > 0)
+    {
+      new_data = new Complex [new_len];
+
+      int min_len = new_len < len ? new_len : len;
+
+      for (int i = 0; i < min_len; i++)
+	new_data[i] = data[i];
+
+      for (i = min_len; i < new_len; i++)
+	new_data[i] = val;
+    }
+
+  delete [] data;
+  nr = r;
+  nc = c;
+  len = new_len;
+  data = new_data;
+
+  return *this;
+}
+#endif
+
+int
+ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const
+{
+  if (rows () != a.rows () || cols () != a.cols ())
+    return 0;
+
+  return equal (data (), a.data (), length ());
+}
+
+int
+ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const
+{
+  return !(*this == a);
+}
+
+ComplexDiagMatrix
+ComplexDiagMatrix::hermitian (void) const
+{
+  return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ());
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (double val)
+{
+  for (int i = 0; i < length (); i++)
+    elem (i, i) = val;
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const Complex& val)
+{
+  for (int i = 0; i < length (); i++)
+    elem (i, i) = val;
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (double val, int beg, int end)
+{
+  if (beg < 0 || end >= length () || end < beg)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = beg; i < end; i++)
+    elem (i, i) = val;
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const Complex& val, int beg, int end)
+{
+  if (beg < 0 || end >= length () || end < beg)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = beg; i < end; i++)
+    elem (i, i) = val;
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const ColumnVector& a)
+{
+  int len = length ();
+  if (a.length () != len)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < len; i++)
+    elem (i, i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const ComplexColumnVector& a)
+{
+  int len = length ();
+  if (a.length () != len)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < len; i++)
+    elem (i, i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const RowVector& a)
+{
+  int len = length ();
+  if (a.length () != len)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < len; i++)
+    elem (i, i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const ComplexRowVector& a)
+{
+  int len = length ();
+  if (a.length () != len)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < len; i++)
+    elem (i, i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const ColumnVector& a, int beg)
+{
+  int a_len = a.length ();
+  if (beg < 0 || beg + a_len >= length ())
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (i+beg, i+beg) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const ComplexColumnVector& a, int beg)
+{
+  int a_len = a.length ();
+  if (beg < 0 || beg + a_len >= length ())
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (i+beg, i+beg) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const RowVector& a, int beg)
+{
+  int a_len = a.length ();
+  if (beg < 0 || beg + a_len >= length ())
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (i+beg, i+beg) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::fill (const ComplexRowVector& a, int beg)
+{
+  int a_len = a.length ();
+  if (beg < 0 || beg + a_len >= length ())
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (i+beg, i+beg) = a.elem (i);
+
+  return *this;
+}
+
+ComplexDiagMatrix
+ComplexDiagMatrix::transpose (void) const
+{
+  return ComplexDiagMatrix (dup (data (), length ()), cols (), rows ());
+}
+
+DiagMatrix
+real (const ComplexDiagMatrix& a)
+{
+  DiagMatrix retval;
+  int a_len = a.length ();
+  if (a_len > 0)
+    retval = DiagMatrix (real_dup (a.data (), a_len), a.rows (),
+			 a.cols ());
+  return retval;
+}
+
+DiagMatrix
+imag (const ComplexDiagMatrix& a)
+{
+  DiagMatrix retval;
+  int a_len = a.length ();
+  if (a_len > 0)
+    retval = DiagMatrix (imag_dup (a.data (), a_len), a.rows (),
+			 a.cols ());
+  return retval;
+}
+
+ComplexDiagMatrix
+conj (const ComplexDiagMatrix& a)
+{
+  ComplexDiagMatrix retval;
+  int a_len = a.length ();
+  if (a_len > 0)
+    retval = ComplexDiagMatrix (conj_dup (a.data (), a_len),
+				a.rows (), a.cols ());
+  return retval;
+}
+
+// resize is the destructive analog for this one
+
+ComplexMatrix
+ComplexDiagMatrix::extract (int r1, int c1, int r2, int c2) const
+{
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
+
+  int new_r = r2 - r1 + 1;
+  int new_c = c2 - c1 + 1;
+
+  ComplexMatrix result (new_r, new_c);
+
+  for (int j = 0; j < new_c; j++)
+    for (int i = 0; i < new_r; i++)
+      result.elem (i, j) = elem (r1+i, c1+j);
+
+  return result;
+}
+
+// extract row or column i.
+
+ComplexRowVector
+ComplexDiagMatrix::row (int i) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (i < 0 || i >= nr)
+    {
+      (*current_liboctave_error_handler) ("invalid row selection");
+      return RowVector (); 
+    }
+
+  ComplexRowVector retval (nc, 0.0);
+  if (nr <= nc || (nr > nc && i < nc))
+    retval.elem (i) = elem (i, i);
+
+  return retval;
+}
+
+ComplexRowVector
+ComplexDiagMatrix::row (char *s) const
+{
+  if (s == (char *) NULL)
+    {
+      (*current_liboctave_error_handler) ("invalid row selection");
+      return ComplexRowVector (); 
+    }
+
+  char c = *s;
+  if (c == 'f' || c == 'F')
+    return row (0);
+  else if (c == 'l' || c == 'L')
+    return row (rows () - 1);
+  else
+    {
+      (*current_liboctave_error_handler) ("invalid row selection");
+      return ComplexRowVector ();
+    }
+}
+
+ComplexColumnVector
+ComplexDiagMatrix::column (int i) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (i < 0 || i >= nc)
+    {
+      (*current_liboctave_error_handler) ("invalid column selection");
+      return ColumnVector (); 
+    }
+
+  ComplexColumnVector retval (nr, 0.0);
+  if (nr >= nc || (nr < nc && i < nr))
+    retval.elem (i) = elem (i, i);
+
+  return retval;
+}
+
+ComplexColumnVector
+ComplexDiagMatrix::column (char *s) const
+{
+  if (s == (char *) NULL)
+    {
+      (*current_liboctave_error_handler) ("invalid column selection");
+      return ColumnVector (); 
+    }
+
+  char c = *s;
+  if (c == 'f' || c == 'F')
+    return column (0);
+  else if (c == 'l' || c == 'L')
+    return column (cols () - 1);
+  else
+    {
+      (*current_liboctave_error_handler) ("invalid column selection");
+      return ColumnVector (); 
+    }
+}
+
+ComplexDiagMatrix
+ComplexDiagMatrix::inverse (void) const
+{
+  int info;
+  return inverse (info);
+}
+
+ComplexDiagMatrix
+ComplexDiagMatrix::inverse (int& info) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != nc)
+    {
+      (*current_liboctave_error_handler) ("inverse requires square matrix");
+      return DiagMatrix ();
+    }
+
+  ComplexDiagMatrix retval (nr, nc);
+
+  info = 0;
+  for (int i = 0; i < length (); i++)
+    {
+      if (elem (i, i) == 0.0)
+	{
+	  info = -1;
+	  return *this;
+	}
+      else
+	retval.elem (i, i) = 1.0 / elem (i, i);
+    }
+
+  return *this;
+}
+
+// diagonal matrix by diagonal matrix -> diagonal matrix operations
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::operator += (const DiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix += operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  add2 (d, a.data (), length ());
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::operator -= (const DiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix -= operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  subtract2 (d, a.data (), length ());
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::operator += (const ComplexDiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix += operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  add2 (d, a.data (), length ());
+  return *this;
+}
+
+ComplexDiagMatrix&
+ComplexDiagMatrix::operator -= (const ComplexDiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix -= operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  subtract2 (d, a.data (), length ());
+  return *this;
+}
+
+// diagonal matrix by scalar -> matrix operations
+
+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;
+}
+
+// diagonal matrix by scalar -> diagonal matrix operations
+
+ComplexDiagMatrix
+operator * (const ComplexDiagMatrix& a, double s)
+{
+  return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
+			    a.rows (), a.cols ());
+}
+
+ComplexDiagMatrix
+operator / (const ComplexDiagMatrix& a, double s)
+{
+  return ComplexDiagMatrix (divide (a.data (), a.length (), s),
+			    a.rows (), a.cols ());
+}
+
+// scalar by diagonal matrix -> matrix operations
+
+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;
+}
+
+// scalar by diagonal matrix -> diagonal matrix operations
+
+ComplexDiagMatrix
+operator * (double s, const ComplexDiagMatrix& a)
+{
+  return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
+			    a.rows (), a.cols ());
+}
+
+// diagonal matrix by column vector -> column vector operations
+
+ComplexColumnVector
+operator * (const ComplexDiagMatrix& m, const ColumnVector& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  int a_len = a.length ();
+  if (nc != a_len)
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix muliplication attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (nc == 0 || nr == 0)
+    return ComplexColumnVector (0);
+
+  ComplexColumnVector result (nr);
+
+  for (int i = 0; i < a_len; i++)
+    result.elem (i) = a.elem (i) * m.elem (i, i);
+
+  for (i = a_len; i < nr; i++)
+    result.elem (i) = 0.0;
+
+  return result;
+}
+
+ComplexColumnVector
+operator * (const ComplexDiagMatrix& m, const ComplexColumnVector& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  int a_len = a.length ();
+  if (nc != a_len)
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix muliplication attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (nc == 0 || nr == 0)
+    return ComplexColumnVector (0);
+
+  ComplexColumnVector result (nr);
+
+  for (int i = 0; i < a_len; i++)
+    result.elem (i) = a.elem (i) * m.elem (i, i);
+
+  for (i = a_len; i < nr; i++)
+    result.elem (i) = 0.0;
+
+  return result;
+}
+
+// diagonal matrix by diagonal matrix -> diagonal matrix operations
+
+ComplexDiagMatrix
+operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
+{
+  int nr_a = a.rows ();
+  int nc_a = a.cols ();
+  int nr_b = b.rows ();
+  int nc_b = b.cols ();
+  if (nc_a != nr_b)
+    {
+      (*current_liboctave_error_handler)
+        ("nonconformant matrix multiplication attempted");
+      return ComplexDiagMatrix ();
+    }
+
+  if (nr_a == 0 || nc_a == 0 || nc_b == 0)
+    return ComplexDiagMatrix (nr_a, nc_a, 0.0);
+
+  ComplexDiagMatrix c (nr_a, nc_b);
+
+  int len = nr_a < nc_b ? nr_a : nc_b;
+
+  for (int i = 0; i < len; i++)
+    {
+      Complex a_element = a.elem (i, i);
+      Complex b_element = b.elem (i, i);
+
+      if (a_element == 0.0 || b_element == 0.0)
+        c.elem (i, i) = 0.0;
+      else if (a_element == 1.0)
+        c.elem (i, i) = b_element;
+      else if (b_element == 1.0)
+        c.elem (i, i) = a_element;
+      else
+        c.elem (i, i) = a_element * b_element;
+    }
+
+  return c;
+}
+
+ComplexDiagMatrix
+operator + (const ComplexDiagMatrix& m, const DiagMatrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix addition attempted");
+      return ComplexDiagMatrix ();
+    }
+
+  if (nr == 0 || nc == 0)
+    return ComplexDiagMatrix (nr, nc);
+
+  return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
+}
+
+ComplexDiagMatrix
+operator - (const ComplexDiagMatrix& m, const DiagMatrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix subtraction attempted");
+      return ComplexDiagMatrix ();
+    }
+
+  if (nr == 0 || nc == 0)
+    return ComplexDiagMatrix (nr, nc);
+
+  return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()),
+			    nr, nc);
+}
+
+ComplexDiagMatrix
+operator * (const ComplexDiagMatrix& a, const DiagMatrix& b)
+{
+  int nr_a = a.rows ();
+  int nc_a = a.cols ();
+  int nr_b = b.rows ();
+  int nc_b = b.cols ();
+  if (nc_a != nr_b)
+    {
+      (*current_liboctave_error_handler)
+        ("nonconformant matrix multiplication attempted");
+      return ComplexDiagMatrix ();
+    }
+
+  if (nr_a == 0 || nc_a == 0 || nc_b == 0)
+    return ComplexDiagMatrix (nr_a, nc_a, 0.0);
+
+  ComplexDiagMatrix c (nr_a, nc_b);
+
+  int len = nr_a < nc_b ? nr_a : nc_b;
+
+  for (int i = 0; i < len; i++)
+    {
+      Complex a_element = a.elem (i, i);
+      double b_element = b.elem (i, i);
+
+      if (a_element == 0.0 || b_element == 0.0)
+        c.elem (i, i) = 0.0;
+      else if (a_element == 1.0)
+        c.elem (i, i) = b_element;
+      else if (b_element == 1.0)
+        c.elem (i, i) = a_element;
+      else
+        c.elem (i, i) = a_element * b_element;
+    }
+
+  return c;
+}
+
+ComplexDiagMatrix
+product (const ComplexDiagMatrix& m, const DiagMatrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix product attempted");
+      return ComplexDiagMatrix ();
+    }
+
+  if (nr == 0 || nc == 0)
+    return ComplexDiagMatrix (nr, nc);
+
+  return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
+			    nr, nc);
+}
+
+// diagonal matrix by matrix -> matrix operations
+
+ComplexMatrix
+operator + (const ComplexDiagMatrix& m, const Matrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix addition attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix subtraction attempted");
+      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)
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix multiplication attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix addition attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix subtraction attempted");
+      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)
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix multiplication attempted");
+      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;
+}
+
+// other operations
+
+ComplexColumnVector
+ComplexDiagMatrix::diag (void) const
+{
+  return diag (0);
+}
+
+// Could be optimized...
+
+ComplexColumnVector
+ComplexDiagMatrix::diag (int k) const
+{
+  int nnr = rows ();
+  int nnc = cols ();
+  if (k > 0)
+    nnc -= k;
+  else if (k < 0)
+    nnr += k;
+
+  ComplexColumnVector d;
+
+  if (nnr > 0 && nnc > 0)
+    {
+      int ndiag = (nnr < nnc) ? nnr : nnc;
+
+      d.resize (ndiag);
+
+      if (k > 0)
+	{
+	  for (int i = 0; i < ndiag; i++)
+	    d.elem (i) = elem (i, i+k);
+	}
+      else if ( k < 0)
+	{
+	  for (int i = 0; i < ndiag; i++)
+	    d.elem (i) = elem (i-k, i);
+	}
+      else
+	{
+	  for (int i = 0; i < ndiag; i++)
+	    d.elem (i) = elem (i, i);
+	}
+    }
+  else
+    cerr << "diag: requested diagonal out of range\n";
+
+  return d;
+}
+
+// i/o
+
+ostream&
+operator << (ostream& os, const ComplexDiagMatrix& a)
+{
+  Complex ZERO (0.0);
+//  int field_width = os.precision () + 7;
+  for (int i = 0; i < a.rows (); i++)
+    {
+      for (int j = 0; j < a.cols (); j++)
+	{
+	  if (i == j)
+	    os << " " /* setw (field_width) */ << a.elem (i, i);
+	  else
+	    os << " " /* setw (field_width) */ << ZERO;
+	}
+      os << "\n";
+    }
+  return os;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/CDiagMatrix.h
@@ -0,0 +1,218 @@
+//                                  -*- C++ -*-
+/*
+
+Copyright (C) 1992, 1993, 1994 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#if !defined (octave_ComplexDiagMatrix_h)
+#define octave_ComplexDiagMatrix_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include "Array.h"
+
+#include "dRowVector.h"
+#include "CRowVector.h"
+#include "dColVector.h"
+#include "CColVector.h"
+
+#include "mx-defs.h"
+
+extern "C++" {
+
+class ComplexDiagMatrix : public DiagArray<Complex>
+{
+public:
+
+  ComplexDiagMatrix (void) : DiagArray<Complex> () { }
+  ComplexDiagMatrix (int n) : DiagArray<Complex> (n) { }
+  ComplexDiagMatrix (int n, const Complex& val)
+    : DiagArray<Complex> (n, val) { }
+  ComplexDiagMatrix (int r, int c) : DiagArray<Complex> (r, c) { }
+  ComplexDiagMatrix (int r, int c, const Complex& val)
+    : DiagArray<Complex> (r, c, val) { }
+  ComplexDiagMatrix (const RowVector& a);
+  ComplexDiagMatrix (const ComplexRowVector& a) : DiagArray<Complex> (a) { }
+  ComplexDiagMatrix (const ColumnVector& a);
+  ComplexDiagMatrix (const ComplexColumnVector& a)
+    : DiagArray<Complex> (a) { }
+  ComplexDiagMatrix (const DiagMatrix& a);
+  ComplexDiagMatrix (const DiagArray<Complex>& a)
+    : DiagArray<Complex> (a) { }
+  ComplexDiagMatrix (const ComplexDiagMatrix& a) : DiagArray<Complex> (a) { }
+//  ComplexDiagMatrix (const Complex& a) : DiagArray<Complex> (1, a) { }
+
+  ComplexDiagMatrix& operator = (const ComplexDiagMatrix& a)
+    {
+      DiagArray<Complex>::operator = (a);
+      return *this;
+    }
+
+//  operator DiagArray<Complex>& () const { return *this; }
+
+  int operator == (const ComplexDiagMatrix& a) const;
+  int operator != (const ComplexDiagMatrix& a) const;
+
+  ComplexDiagMatrix& fill (double val);
+  ComplexDiagMatrix& fill (const Complex& val);
+  ComplexDiagMatrix& fill (double val, int beg, int end);
+  ComplexDiagMatrix& fill (const Complex& val, int beg, int end);
+  ComplexDiagMatrix& fill (const ColumnVector& a);
+  ComplexDiagMatrix& fill (const ComplexColumnVector& a);
+  ComplexDiagMatrix& fill (const RowVector& a);
+  ComplexDiagMatrix& fill (const ComplexRowVector& a);
+  ComplexDiagMatrix& fill (const ColumnVector& a, int beg);
+  ComplexDiagMatrix& fill (const ComplexColumnVector& a, int beg);
+  ComplexDiagMatrix& fill (const RowVector& a, int beg);
+  ComplexDiagMatrix& fill (const ComplexRowVector& a, int beg);
+
+  ComplexDiagMatrix hermitian (void) const;  // complex conjugate transpose
+  ComplexDiagMatrix transpose (void) const;
+
+  friend DiagMatrix real (const ComplexDiagMatrix& a);
+  friend DiagMatrix imag (const ComplexDiagMatrix& a);
+  friend ComplexDiagMatrix conj (const ComplexDiagMatrix& a);
+
+// resize is the destructive analog for this one
+
+  ComplexMatrix extract (int r1, int c1, int r2, int c2) const;
+
+// extract row or column i.
+
+  ComplexRowVector row (int i) const;
+  ComplexRowVector row (char *s) const;
+
+  ComplexColumnVector column (int i) const;
+  ComplexColumnVector column (char *s) const;
+
+  ComplexDiagMatrix inverse (int& info) const;
+  ComplexDiagMatrix inverse (void) const;
+
+// diagonal matrix by diagonal matrix -> diagonal matrix operations
+
+  ComplexDiagMatrix& operator += (const DiagMatrix& a);
+  ComplexDiagMatrix& operator -= (const DiagMatrix& a);
+
+  ComplexDiagMatrix& operator += (const ComplexDiagMatrix& a);
+  ComplexDiagMatrix& operator -= (const ComplexDiagMatrix& a);
+
+// diagonal matrix by scalar -> matrix operations
+
+  friend ComplexMatrix operator + (const ComplexDiagMatrix& a, double s);
+  friend ComplexMatrix operator - (const ComplexDiagMatrix& a, double s);
+
+  friend ComplexMatrix operator + (const ComplexDiagMatrix& a,
+				   const Complex& s);
+  friend ComplexMatrix operator - (const ComplexDiagMatrix& a,
+				   const Complex& s);
+
+// diagonal matrix by scalar -> diagonal matrix operations
+
+  friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, double s);
+  friend ComplexDiagMatrix operator / (const ComplexDiagMatrix& a, double s);
+
+// scalar by diagonal matrix -> matrix operations
+
+  friend ComplexMatrix operator + (double s, const ComplexDiagMatrix& a);
+  friend ComplexMatrix operator - (double s, const ComplexDiagMatrix& a);
+
+  friend ComplexMatrix operator + (const Complex& s,
+				   const ComplexDiagMatrix& a);
+  friend ComplexMatrix operator - (const Complex& s,
+				   const ComplexDiagMatrix& a);
+
+// scalar by diagonal matrix -> diagonal matrix operations
+
+  friend ComplexDiagMatrix operator * (double s, const ComplexDiagMatrix& a);
+
+// diagonal matrix by column vector -> column vector operations
+
+  friend ComplexColumnVector operator * (const ComplexDiagMatrix& a,
+					 const ColumnVector& b);
+
+  friend ComplexColumnVector operator * (const ComplexDiagMatrix& a,
+					 const ComplexColumnVector& b);
+
+// diagonal matrix by diagonal matrix -> diagonal matrix operations
+
+  friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a,
+				       const ComplexDiagMatrix& b);
+
+  friend ComplexDiagMatrix operator + (const ComplexDiagMatrix& a,
+				       const DiagMatrix& b);
+  friend ComplexDiagMatrix operator - (const ComplexDiagMatrix& a,
+				       const DiagMatrix& b);
+  friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a,
+				       const DiagMatrix& b);
+
+  friend ComplexDiagMatrix product (const ComplexDiagMatrix& a,
+				    const DiagMatrix& b); 
+
+// diagonal matrix by matrix -> matrix operations
+
+  friend ComplexMatrix operator + (const ComplexDiagMatrix& a,
+				   const Matrix& b); 
+  friend ComplexMatrix operator - (const ComplexDiagMatrix& a,
+				   const Matrix& b);
+  friend ComplexMatrix operator * (const ComplexDiagMatrix& a,
+				   const Matrix& b);
+
+  friend ComplexMatrix operator + (const ComplexDiagMatrix& a,
+				   const ComplexMatrix& b);
+  friend ComplexMatrix operator - (const ComplexDiagMatrix& a,
+				   const ComplexMatrix& b);
+  friend ComplexMatrix operator * (const ComplexDiagMatrix& a,
+				   const ComplexMatrix& b);
+
+// other operations
+
+  ComplexColumnVector diag (void) const;
+  ComplexColumnVector diag (int k) const;
+
+// i/o
+
+  friend ostream& operator << (ostream& os, const ComplexDiagMatrix& a);
+
+#define KLUDGE_DIAG_MATRICES
+#define TYPE Complex
+#define KL_DMAT_TYPE ComplexDiagMatrix
+#include "mx-kludge.h"
+#undef KLUDGE_DIAG_MATRICES
+#undef TYPE
+#undef KL_DMAT_TYPE
+
+private:
+
+  ComplexDiagMatrix (Complex *d, int nr, int nc)
+    : DiagArray<Complex> (d, nr, nc) { }
+};
+
+} // extern "C++"
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/CMatrix.cc
@@ -0,0 +1,2613 @@
+// Matrix manipulations.                              -*- C++ -*-
+/*
+
+Copyright (C) 1992, 1993, 1994 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#include <sys/types.h>
+#include <iostream.h>
+
+#include <Complex.h>
+
+#include "mx-base.h"
+#include "CmplxDET.h"
+#include "mx-inlines.cc"
+#include "lo-error.h"
+#include "f77-uscore.h"
+
+// Fortran functions we call.
+
+extern "C"
+{
+  int F77_FCN (zgemm) (const char*, const char*, const int*,
+		       const int*, const int*, const Complex*,
+		       const Complex*, const int*, const Complex*,
+		       const int*, const Complex*, Complex*, const int*,
+		       long, long);
+
+  int F77_FCN (zgemv) (const char*, const int*, const int*,
+		       const Complex*, const Complex*, const int*,
+		       const Complex*, const int*, const Complex*,
+		       Complex*, const int*, long);
+
+  int F77_FCN (zgeco) (Complex*, const int*, const int*, int*,
+		       double*, Complex*);
+
+  int F77_FCN (zgedi) (Complex*, const int*, const int*, int*,
+		       Complex*, Complex*, const int*);
+
+  int F77_FCN (zgesl) (Complex*, const int*, const int*, int*,
+		       Complex*, const int*);
+
+  int F77_FCN (zgelss) (const int*, const int*, const int*, Complex*,
+			const int*, Complex*, const int*, double*,
+			const double*, int*, Complex*, const int*,
+			double*, int*);
+
+// Note that the original complex fft routines were not written for
+// double complex arguments.  They have been modified by adding an
+// implicit double precision (a-h,o-z) statement at the beginning of
+// each subroutine.
+
+  int F77_FCN (cffti) (const int*, Complex*);
+
+  int F77_FCN (cfftf) (const int*, Complex*, Complex*);
+
+  int F77_FCN (cfftb) (const int*, Complex*, Complex*);
+}
+
+#define KLUDGE_MATRICES
+#define TYPE Complex
+#define KL_MAT_TYPE ComplexMatrix
+#include "mx-kludge.cc"
+#undef KLUDGE_MATRICES
+#undef TYPE
+#undef KL_MAT_TYPE
+
+/*
+ * Complex Matrix class
+ */
+
+ComplexMatrix::ComplexMatrix (const Matrix& a)
+  : Array2<Complex> (a.rows (), a.cols ())
+{
+  for (int j = 0; j < cols (); j++)
+    for (int i = 0; i < rows (); i++)
+      elem (i, j) = a.elem (i, j);
+}
+
+ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
+  : Array2<Complex> (a.rows (), a.cols (), 0.0)
+{
+  for (int i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
+  : Array2<Complex> (a.rows (), a.cols (), 0.0)
+{
+  for (int i = 0; i < a.length (); i++)
+    elem (i, i) = a.elem (i, i);
+}
+
+#if 0
+ComplexMatrix&
+ComplexMatrix::resize (int r, int c)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimensions");
+      return *this;
+    }
+
+  int new_len = r * c;
+  Complex* new_data = (Complex *) NULL;
+  if (new_len > 0)
+    {
+      new_data = new Complex [new_len];
+
+      int min_r = nr < r ? nr : r;
+      int min_c = nc < c ? nc : c;
+
+      for (int j = 0; j < min_c; j++)
+	for (int i = 0; i < min_r; i++)
+	  new_data[r*j+i] = elem (i, j);
+    }
+
+  delete [] data;
+  nr = r;
+  nc = c;
+  len = new_len;
+  data = new_data;
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::resize (int r, int c, double val)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimensions");
+      return *this;
+    }
+
+  int new_len = r * c;
+  Complex *new_data = (Complex *) NULL;
+  if (new_len > 0)
+    {
+      new_data = new Complex [new_len];
+
+// There may be faster or cleaner ways to do this.
+
+      if (r > nr || c > nc)
+	copy (new_data, new_len, val);
+
+      int min_r = nr < r ? nr : r;
+      int min_c = nc < c ? nc : c;
+
+      for (int j = 0; j < min_c; j++)
+	for (int i = 0; i < min_r; i++)
+	  new_data[r*j+i] = elem (i, j);
+    }
+
+  delete [] data;
+  nr = r;
+  nc = c;
+  len = new_len;
+  data = new_data;
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::resize (int r, int c, const Complex& val)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler)
+	("can't resize to negative dimensions");
+      return *this;
+    }
+
+  int new_len = r * c;
+  Complex *new_data = (Complex *) NULL;
+  if (new_len > 0)
+    {
+      new_data = new Complex [new_len];
+
+// There may be faster or cleaner ways to do this.
+
+      if (r > nr || c > nc)
+	copy (new_data, new_len, val);
+
+      int min_r = nr < r ? nr : r;
+      int min_c = nc < c ? nc : c;
+
+      for (int j = 0; j < min_c; j++)
+	for (int i = 0; i < min_r; i++)
+	  new_data[r*j+i] = elem (i, j);
+    }
+
+  delete [] data;
+  nr = r;
+  nc = c;
+  len = new_len;
+  data = new_data;
+
+  return *this;
+}
+#endif
+
+int
+ComplexMatrix::operator == (const ComplexMatrix& a) const
+{
+  if (rows () != a.rows () || cols () != a.cols ())
+    return 0;
+
+  return equal (data (), a.data (), length ());
+}
+
+int
+ComplexMatrix::operator != (const ComplexMatrix& a) const
+{
+  return !(*this == a);
+}
+
+// destructive insert/delete/reorder operations
+
+ComplexMatrix&
+ComplexMatrix::insert (const Matrix& a, int r, int c)
+{
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+  if (r < 0 || r + a_nr - 1 > rows () || c < 0 || c + a_nc - 1 > cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int j = 0; j < a_nc; j++)
+    for (int i = 0; i < a_nr; i++)
+      elem (r+i, c+j) = a.elem (i, j);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const RowVector& a, int r, int c)
+{
+  int a_len = a.length ();
+  if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (r, c+i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const ColumnVector& a, int r, int c)
+{
+  int a_len = a.length ();
+  if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (r+i, c) = a.elem (i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const DiagMatrix& a, int r, int c)
+{
+  if (r < 0 || r + a.rows () - 1 > rows ()
+      || c < 0 || c + a.cols () - 1 > cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a.length (); i++)
+    elem (r+i, c+i) = a.elem (i, i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const ComplexMatrix& a, int r, int c)
+{
+  int a_nr = a.rows ();
+  int a_nc = a.cols ();
+  if (r < 0 || r + a_nr - 1 > rows () || c < 0 || c + a_nc - 1 > cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int j = 0; j < a_nc; j++)
+    for (int i = 0; i < a_nr; i++)
+      elem (r+i, c+j) = a.elem (i, j);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const ComplexRowVector& a, int r, int c)
+{
+  int a_len = a.length ();
+  if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (r, c+i) = a.elem (i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c)
+{
+  int a_len = a.length ();
+  if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a_len; i++)
+    elem (r+i, c) = a.elem (i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c)
+{
+  if (r < 0 || r + a.rows () - 1 > rows ()
+      || c < 0 || c + a.cols () - 1 > cols ())
+    {
+      (*current_liboctave_error_handler) ("range error for insert");
+      return *this;
+    }
+
+  for (int i = 0; i < a.length (); i++)
+    elem (r+i, c+i) = a.elem (i, i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::fill (double val)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr > 0 && nc > 0)
+    for (int j = 0; j < nc; j++)
+      for (int i = 0; i < nr; i++)
+	elem (i, j) = val;
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::fill (const Complex& val)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr > 0 && nc > 0)
+    for (int j = 0; j < nc; j++)
+      for (int i = 0; i < nr; i++)
+	elem (i, j) = val;
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
+      || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
+
+  for (int j = c1; j <= c2; j++)
+    for (int i = r1; i <= r2; i++)
+      elem (i, j) = val;
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
+      || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
+    {
+      (*current_liboctave_error_handler) ("range error for fill");
+      return *this;
+    }
+
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
+
+  for (int j = c1; j <= c2; j++)
+    for (int i = r1; i <= r2; i++)
+      elem (i, j) = val;
+
+  return *this;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const Matrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows ())
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + a.cols ());
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const RowVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != 1)
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + a.length ());
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const ColumnVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.length ())
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + 1);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const DiagMatrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows ())
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + a.cols ());
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const ComplexMatrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows ())
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + a.cols ());
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const ComplexRowVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != 1)
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + a.length ());
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const ComplexColumnVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.length ())
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + 1);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::append (const ComplexDiagMatrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows ())
+    {
+      (*current_liboctave_error_handler) ("row dimension mismatch for append");
+      return *this;
+    }
+
+  int nc_insert = nc;
+  ComplexMatrix retval (nr, nc + a.cols ());
+  retval.insert (*this, 0, 0);
+  retval.insert (a, 0, nc_insert);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const Matrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + a.rows (), nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const RowVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + 1, nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const ColumnVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != 1)
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + a.length (), nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const DiagMatrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + a.rows (), nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const ComplexMatrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + a.rows (), nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const ComplexRowVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + 1, nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const ComplexColumnVector& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != 1)
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + a.length (), nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::stack (const ComplexDiagMatrix& a) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("column dimension mismatch for stack");
+      return *this;
+    }
+
+  int nr_insert = nr;
+  ComplexMatrix retval (nr + a.rows (), nc);
+  retval.insert (*this, 0, 0);
+  retval.insert (a, nr_insert, 0);
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::hermitian (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix result;
+  if (length () > 0)
+    {
+      result.resize (nc, nr);
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  result.elem (j, i) = conj (elem (i, j));
+    }
+  return result;
+}
+
+ComplexMatrix
+ComplexMatrix::transpose (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix result (nc, nr);
+  if (length () > 0)
+    {
+      for (int j = 0; j < nc; j++)
+	for (int i = 0; i < nr; i++)
+	  result.elem (j, i) = elem (i, j);
+    }
+  return result;
+}
+
+Matrix
+real (const ComplexMatrix& a)
+{
+  int a_len = a.length ();
+  Matrix retval;
+  if (a_len > 0)
+    retval = Matrix (real_dup (a.data (), a_len), a.rows (), a.cols ());
+  return retval;
+}
+
+Matrix
+imag (const ComplexMatrix& a)
+{
+  int a_len = a.length ();
+  Matrix retval;
+  if (a_len > 0)
+    retval = Matrix (imag_dup (a.data (), a_len), a.rows (), a.cols ());
+  return retval;
+}
+
+ComplexMatrix
+conj (const ComplexMatrix& a)
+{
+  int a_len = a.length ();
+  ComplexMatrix retval;
+  if (a_len > 0)
+    retval = ComplexMatrix (conj_dup (a.data (), a_len), a.rows (),
+			    a.cols ());
+  return retval;
+}
+
+// resize is the destructive equivalent for this one
+
+ComplexMatrix
+ComplexMatrix::extract (int r1, int c1, int r2, int c2) const
+{
+  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
+  if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
+
+  int new_r = r2 - r1 + 1;
+  int new_c = c2 - c1 + 1;
+
+  ComplexMatrix result (new_r, new_c);
+
+  for (int j = 0; j < new_c; j++)
+    for (int i = 0; i < new_r; i++)
+      result.elem (i, j) = elem (r1+i, c1+j);
+
+  return result;
+}
+
+// extract row or column i.
+
+ComplexRowVector
+ComplexMatrix::row (int i) const
+{
+  int nc = cols ();
+  if (i < 0 || i >= rows ())
+    {
+      (*current_liboctave_error_handler) ("invalid row selection");
+      return ComplexRowVector ();
+    }
+
+  ComplexRowVector retval (nc);
+  for (int j = 0; j < cols (); j++)
+    retval.elem (j) = elem (i, j);
+
+  return retval;
+}
+
+ComplexRowVector
+ComplexMatrix::row (char *s) const
+{
+  if (s == (char *) NULL)
+    {
+      (*current_liboctave_error_handler) ("invalid row selection");
+      return ComplexRowVector ();
+    }
+
+  char c = *s;
+  if (c == 'f' || c == 'F')
+    return row (0);
+  else if (c == 'l' || c == 'L')
+    return row (rows () - 1);
+  else
+    {
+      (*current_liboctave_error_handler) ("invalid row selection");
+      return ComplexRowVector ();
+    }
+}
+
+ComplexColumnVector
+ComplexMatrix::column (int i) const
+{
+  int nr = rows ();
+  if (i < 0 || i >= cols ())
+    {
+      (*current_liboctave_error_handler) ("invalid column selection");
+      return ComplexColumnVector ();
+    }
+
+  ComplexColumnVector retval (nr);
+  for (int j = 0; j < nr; j++)
+    retval.elem (j) = elem (j, i);
+
+  return retval;
+}
+
+ComplexColumnVector
+ComplexMatrix::column (char *s) const
+{
+  if (s == (char *) NULL)
+    {
+      (*current_liboctave_error_handler) ("invalid column selection");
+      return ComplexColumnVector ();
+    }
+
+  char c = *s;
+  if (c == 'f' || c == 'F')
+    return column (0);
+  else if (c == 'l' || c == 'L')
+    return column (cols () - 1);
+  else
+    {
+      (*current_liboctave_error_handler) ("invalid column selection");
+      return ComplexColumnVector ();
+    }
+}
+
+ComplexMatrix
+ComplexMatrix::inverse (void) const
+{
+  int info;
+  double rcond; return inverse (info, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::inverse (int& info) const
+{
+  double rcond;
+  return inverse (info, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::inverse (int& info, double& rcond) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int len = length ();
+  if (nr != nc)
+    {
+      (*current_liboctave_error_handler) ("inverse requires square matrix");
+      return ComplexMatrix ();
+    }
+
+  info = 0;
+
+  int *ipvt = new int [nr];
+  Complex *z = new Complex [nr];
+  Complex *tmp_data = dup (data (), len);
+
+  F77_FCN (zgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z);
+
+  if (rcond + 1.0 == 1.0)
+    {
+      info = -1;
+      copy (tmp_data, data (), len);  // Restore contents.
+    }
+  else
+    {
+      int job = 1;
+      Complex dummy;
+
+      F77_FCN (zgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job);
+    }
+
+  delete [] ipvt;
+  delete [] z;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
+ComplexMatrix
+ComplexMatrix::fourier (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int npts, nsamples;
+  if (nr == 1 || nc == 1)
+    {
+      npts = nr > nc ? nr : nc;
+      nsamples = 1;
+    }
+  else
+    {
+      npts = nr;
+      nsamples = nc;
+    }
+
+  int nn = 4*npts+15;
+  Complex *wsave = new Complex [nn];
+  Complex *tmp_data = dup (data (), length ());
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (int j = 0; j < nsamples; j++)
+    F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave);
+
+  delete [] wsave;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
+ComplexMatrix
+ComplexMatrix::ifourier (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  int npts, nsamples;
+  if (nr == 1 || nc == 1)
+    {
+      npts = nr > nc ? nr : nc;
+      nsamples = 1;
+    }
+  else
+    {
+      npts = nr;
+      nsamples = nc;
+    }
+
+  int nn = 4*npts+15;
+  Complex *wsave = new Complex [nn];
+  Complex *tmp_data = dup (data (), length ());
+
+  F77_FCN (cffti) (&npts, wsave);
+
+  for (int j = 0; j < nsamples; j++)
+    F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave);
+
+  for (j = 0; j < npts*nsamples; j++)
+    tmp_data[j] = tmp_data[j] / (double) npts;
+
+  delete [] wsave;
+
+  return ComplexMatrix (tmp_data, nr, nc);
+}
+
+ComplexDET
+ComplexMatrix::determinant (void) const
+{
+  int info;
+  double rcond;
+  return determinant (info, rcond);
+}
+
+ComplexDET
+ComplexMatrix::determinant (int& info) const
+{
+  double rcond;
+  return determinant (info, rcond);
+}
+
+ComplexDET
+ComplexMatrix::determinant (int& info, double& rcond) const
+{
+  ComplexDET retval;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr == 0 || nc == 0)
+    {
+      Complex d[2];
+      d[0] = 1.0;
+      d[1] = 0.0;
+      retval = ComplexDET (d);
+    }
+  else
+    {
+      info = 0;
+      int *ipvt = new int [nr];
+
+      Complex *z = new Complex [nr];
+      Complex *tmp_data = dup (data (), length ());
+
+      F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
+
+      if (rcond + 1.0 == 1.0)
+	{
+	  info = -1;
+	  retval = ComplexDET ();
+	}
+      else
+	{
+	  int job = 10;
+	  Complex d[2];
+	  F77_FCN (zgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job);
+	  retval = ComplexDET (d);
+	}
+
+      delete [] tmp_data;
+      delete [] ipvt;
+      delete [] z;
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::solve (const Matrix& b) const
+{
+  int info;
+  double rcond;
+  return solve (b, info, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::solve (const Matrix& b, int& info) const
+{
+  double rcond;
+  return solve (b, info, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const
+{
+  ComplexMatrix tmp (b);
+  return solve (tmp, info, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::solve (const ComplexMatrix& b) const
+{
+  int info;
+  double rcond;
+  return solve (b, info, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::solve (const ComplexMatrix& b, int& info) const
+{
+  double rcond;
+  return solve (b, info, rcond);
+}
+ComplexMatrix
+ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
+{
+  ComplexMatrix retval;
+
+  int nr = rows ();
+  int nc = cols ();
+  int b_nr = b.rows ();
+  int b_nc = b.cols ();
+  if (nr == 0 || nc == 0 || nr != nc || nr != b_nr)
+    {
+      (*current_liboctave_error_handler)
+	("matrix dimension mismatch in solution of linear equations");
+      return ComplexMatrix ();
+    }
+
+  info = 0;
+  int *ipvt = new int [nr];
+
+  Complex *z = new Complex [nr];
+  Complex *tmp_data = dup (data (), length ());
+
+  F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
+
+  if (rcond + 1.0 == 1.0)
+    {
+      info = -2;
+    }
+  else
+    {
+      int job = 0;
+
+      Complex *result = dup (b.data (), b.length ());
+
+      for (int j = 0; j < b_nc; j++)
+	F77_FCN (zgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job);
+
+      retval = ComplexMatrix (result, b_nr, b_nc);
+    }
+
+  delete [] tmp_data;
+  delete [] ipvt;
+  delete [] z;
+
+  return retval;
+}
+
+ComplexColumnVector
+ComplexMatrix::solve (const ComplexColumnVector& b) const
+{
+  int info;
+  double rcond;
+  return solve (b, info, rcond);
+}
+
+ComplexColumnVector
+ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
+{
+  double rcond;
+  return solve (b, info, rcond);
+}
+
+ComplexColumnVector
+ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
+		      double& rcond) const
+{
+  ComplexColumnVector retval;
+
+  int nr = rows ();
+  int nc = cols ();
+  int b_len = b.length ();
+  if (nr == 0 || nc == 0 || nr != nc || nr != b_len)
+    {
+      (*current_liboctave_error_handler)
+	("matrix dimension mismatch in solution of linear equations");
+      return ComplexColumnVector ();
+    }
+
+  info = 0;
+  int *ipvt = new int [nr];
+
+  Complex *z = new Complex [nr];
+  Complex *tmp_data = dup (data (), length ());
+
+  F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
+
+  if (rcond + 1.0 == 1.0)
+    {
+      info = -2;
+    }
+  else
+    {
+      int job = 0;
+
+      Complex *result = dup (b.data (), b_len);
+
+      F77_FCN (zgesl) (tmp_data, &nr, &nr, ipvt, result, &job);
+
+      retval = ComplexColumnVector (result, b_len);
+    }
+
+  delete [] tmp_data;
+  delete [] ipvt;
+  delete [] z;
+
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::lssolve (const ComplexMatrix& b) const
+{
+  int info;
+  int rank;
+  return lssolve (b, info, rank);
+}
+
+ComplexMatrix
+ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const
+{
+  int rank;
+  return lssolve (b, info, rank);
+}
+
+ComplexMatrix
+ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
+{
+  int nrhs = b.cols ();
+
+  int m = rows ();
+  int n = cols ();
+
+  if (m == 0 || n == 0 || m != b.rows ())
+    {
+      (*current_liboctave_error_handler)
+	("matrix dimension mismatch solution of linear equations");
+      return Matrix ();
+    }
+
+  Complex *tmp_data = dup (data (), length ());
+
+  int nrr = m > n ? m : n;
+  ComplexMatrix result (nrr, nrhs);
+
+  int i, j;
+  for (j = 0; j < nrhs; j++)
+    for (i = 0; i < m; i++)
+      result.elem (i, j) = b.elem (i, j);
+
+  Complex *presult = result.fortran_vec ();
+
+  int len_s = m < n ? m : n;
+  double *s = new double [len_s];
+  double rcond = -1.0;
+  int lwork;
+  if (m < n)
+    lwork = 2*m + (nrhs > n ? nrhs : n);
+  else
+    lwork = 2*n + (nrhs > m ? nrhs : m);
+
+  Complex *work = new Complex [lwork];
+
+  int lrwork = (5 * (m < n ? m : n)) - 4;
+  lrwork = lrwork > 1 ? lrwork : 1;
+  double *rwork = new double [lrwork];
+
+  F77_FCN (zgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s,
+		    &rcond, &rank, work, &lwork, rwork, &info);
+
+  ComplexMatrix retval (n, nrhs);
+  for (j = 0; j < nrhs; j++)
+    for (i = 0; i < n; i++)
+      retval.elem (i, j) = result.elem (i, j);
+
+  delete [] tmp_data;
+  delete [] s;
+  delete [] work;
+  delete [] rwork;
+
+  return retval;
+}
+
+ComplexColumnVector
+ComplexMatrix::lssolve (const ComplexColumnVector& b) const
+{
+  int info;
+  int rank;
+  return lssolve (b, info, rank);
+}
+
+ComplexColumnVector
+ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const
+{
+  int rank;
+  return lssolve (b, info, rank);
+}
+
+ComplexColumnVector
+ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info,
+			int& rank) const
+{
+  int nrhs = 1;
+
+  int m = rows ();
+  int n = cols ();
+
+  if (m == 0 || n == 0 || m != b.length ())
+    {
+      (*current_liboctave_error_handler)
+	("matrix dimension mismatch solution of least squares problem");
+      return ComplexColumnVector ();
+    }
+
+  Complex *tmp_data = dup (data (), length ());
+
+  int nrr = m > n ? m : n;
+  ComplexColumnVector result (nrr);
+
+  int i;
+  for (i = 0; i < m; i++)
+    result.elem (i) = b.elem (i);
+
+  Complex *presult = result.fortran_vec ();
+
+  int len_s = m < n ? m : n;
+  double *s = new double [len_s];
+  double rcond = -1.0;
+  int lwork;
+  if (m < n)
+    lwork = 2*m + (nrhs > n ? nrhs : n);
+  else
+    lwork = 2*n + (nrhs > m ? nrhs : m);
+
+  Complex *work = new Complex [lwork];
+
+  int lrwork = (5 * (m < n ? m : n)) - 4;
+  lrwork = lrwork > 1 ? lrwork : 1;
+  double *rwork = new double [lrwork];
+
+  F77_FCN (zgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s,
+		    &rcond, &rank, work, &lwork, rwork, &info);
+
+  ComplexColumnVector retval (n);
+  for (i = 0; i < n; i++)
+    retval.elem (i) = result.elem (i);
+
+  delete [] tmp_data;
+  delete [] s;
+  delete [] work;
+  delete [] rwork;
+
+  return retval;
+}
+
+// matrix by diagonal matrix -> matrix operations
+
+ComplexMatrix&
+ComplexMatrix::operator += (const DiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix += operation attempted");
+      return ComplexMatrix ();
+    }
+
+  for (int i = 0; i < a.length (); i++)
+    elem (i, i) += a.elem (i, i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::operator -= (const DiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix -= operation attempted");
+      return ComplexMatrix ();
+    }
+
+  for (int i = 0; i < a.length (); i++)
+    elem (i, i) -= a.elem (i, i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::operator += (const ComplexDiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix += operation attempted");
+      return ComplexMatrix ();
+    }
+
+  for (int i = 0; i < a.length (); i++)
+    elem (i, i) += a.elem (i, i);
+
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix -= operation attempted");
+      return ComplexMatrix ();
+    }
+
+  for (int i = 0; i < a.length (); i++)
+    elem (i, i) -= a.elem (i, i);
+
+  return *this;
+}
+
+// matrix by matrix -> matrix operations
+
+ComplexMatrix&
+ComplexMatrix::operator += (const Matrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix += operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  add2 (d, a.data (), length ());
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::operator -= (const Matrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix -= operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  subtract2 (d, a.data (), length ());
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::operator += (const ComplexMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix += operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  add2 (d, a.data (), length ());
+  return *this;
+}
+
+ComplexMatrix&
+ComplexMatrix::operator -= (const ComplexMatrix& a)
+{
+  int nr = rows ();
+  int nc = cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix -= operation attempted");
+      return *this;
+    }
+
+  if (nr == 0 || nc == 0)
+    return *this;
+
+  Complex *d = fortran_vec (); // Ensures only one reference to my privates!
+
+  subtract2 (d, a.data (), length ());
+  return *this;
+}
+
+// unary operations
+
+Matrix
+ComplexMatrix::operator ! (void) const
+{
+  return Matrix (not (data (), length ()), rows (), cols ());
+}
+
+// matrix by scalar -> matrix operations
+
+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 ());
+}
+
+// matrix by column vector -> column vector operations
+
+ComplexColumnVector
+operator * (const ComplexMatrix& m, const ColumnVector& a)
+{
+  ComplexColumnVector tmp (a);
+  return m * tmp;
+}
+
+ComplexColumnVector
+operator * (const ComplexMatrix& m, const ComplexColumnVector& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nc != a.length ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix multiplication attempted");
+      return ComplexColumnVector ();
+    }
+
+  if (nc == 0 || nr == 0)
+    return ComplexColumnVector (0);
+
+  char trans = 'N';
+  int ld = nr;
+  Complex alpha (1.0);
+  Complex beta (0.0);
+  int i_one = 1;
+
+  Complex *y = new Complex [nr];
+
+  F77_FCN (zgemv) (&trans, &nr, &nc, &alpha, m.data (), &ld, a.data (),
+		   &i_one, &beta, y, &i_one, 1L); 
+
+  return ComplexColumnVector (y, nr);
+}
+
+// matrix by diagonal matrix -> matrix operations
+
+ComplexMatrix
+operator + (const ComplexMatrix& m, const DiagMatrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix addition attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix subtraction attempted");
+      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_nc = a.cols ();
+  if (nc != a.rows ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix multiplication attempted");
+      return ComplexMatrix ();
+    }
+
+  if (nr == 0 || nc == 0 || a_nc == 0)
+    return ComplexMatrix (nr, nc, 0.0);
+
+  Complex *c = new Complex [nr*a_nc];
+  Complex *ctmp = (Complex *) NULL;
+
+  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 ComplexMatrix (c, nr, a_nc);
+}
+
+ComplexMatrix
+operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix addition attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix subtraction attempted");
+      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_nc = a.cols ();
+  if (nc != a.rows ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix multiplication attempted");
+      return ComplexMatrix ();
+    }
+
+  if (nr == 0 || nc == 0 || a_nc == 0)
+    return ComplexMatrix (nr, nc, 0.0);
+
+  Complex *c = new Complex [nr*a_nc];
+  Complex *ctmp = (Complex *) NULL;
+
+  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 ComplexMatrix (c, nr, a_nc);
+}
+
+// matrix by matrix -> matrix operations
+
+ComplexMatrix
+operator + (const ComplexMatrix& m, const Matrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix addition attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix subtraction attempted");
+      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 ComplexMatrix& m, const ComplexMatrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  int a_nc = a.cols ();
+  if (nc != a.rows ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix multiplication attempted");
+      return ComplexMatrix ();
+    }
+
+  if (nr == 0 || nc == 0 || a_nc == 0)
+    return ComplexMatrix (nr, nc, 0.0);
+
+  char trans  = 'N';
+  char transa = 'N';
+
+  int ld  = nr;
+  int lda = a.rows ();
+
+  Complex alpha (1.0);
+  Complex beta (0.0);
+
+  Complex *c = new Complex [nr*a_nc];
+
+  F77_FCN (zgemm) (&trans, &transa, &nr, &a_nc, &nc, &alpha, m.data (),
+		   &ld, a.data (), &lda, &beta, c, &nr, 1L, 1L);
+
+  return ComplexMatrix (c, nr, a_nc);
+}
+
+ComplexMatrix
+product (const ComplexMatrix& m, const Matrix& a)
+{
+  int nr = m.rows ();
+  int nc = m.cols ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix product attempted");
+      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 ();
+  if (nr != a.rows () || nc != a.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant matrix quotient attempted");
+      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
+map (c_c_Mapper f, const ComplexMatrix& a)
+{
+  ComplexMatrix b (a);
+  b.map (f);
+  return b;
+}
+
+Matrix
+map (d_c_Mapper f, const ComplexMatrix& a)
+{
+  int a_nc = a.cols ();
+  int a_nr = a.rows ();
+  Matrix b (a_nr, a_nc);
+  for (int j = 0; j < a_nc; j++)
+    for (int i = 0; i < a_nr; i++)
+      b.elem (i, j) = f (a.elem (i, j));
+  return b;
+}
+
+void
+ComplexMatrix::map (c_c_Mapper f)
+{
+  for (int j = 0; j < cols (); j++)
+    for (int i = 0; i < rows (); i++)
+      elem (i, j) = f (elem (i, j));
+}
+
+Matrix
+ComplexMatrix::all (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  Matrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 1.0;
+	  for (int j = 0; j < nc; j++)
+	    {
+	      if (elem (0, j) == 0.0)
+		{
+		  retval.elem (0, 0) = 0.0;
+		  break;
+		}
+	    }
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 1.0;
+	  for (int i = 0; i < nr; i++)
+	    {
+	      if (elem (i, 0) == 0.0)
+		{
+		  retval.elem (0, 0) = 0.0;
+		  break;
+		}
+	    }
+	}
+      else
+	{
+	  retval.resize (1, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = 1.0;
+	      for (int i = 0; i < nr; i++)
+		{
+		  if (elem (i, j) == 0.0)
+		    {
+		      retval.elem (0, j) = 0.0;
+		      break;
+		    }
+		}
+	    }
+	}
+    }
+  return retval;
+}
+
+Matrix
+ComplexMatrix::any (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  Matrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 0.0;
+	  for (int j = 0; j < nc; j++)
+	    {
+	      if (elem (0, j) != 0.0)
+		{
+		  retval.elem (0, 0) = 1.0;
+		  break;
+		}
+	    }
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 0.0;
+	  for (int i = 0; i < nr; i++)
+	    {
+	      if (elem (i, 0) != 0.0)
+		{
+		  retval.elem (0, 0) = 1.0;
+		  break;
+		}
+	    }
+	}
+      else
+	{
+	  retval.resize (1, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = 0.0;
+	      for (int i = 0; i < nr; i++)
+		{
+		  if (elem (i, j) != 0.0)
+		    {
+		      retval.elem (0, j) = 1.0;
+		      break;
+		    }
+		}
+	    }
+	}
+    }
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::cumprod (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, nc);
+	  Complex prod = elem (0, 0);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = prod;
+	      if (j < nc - 1)
+		prod *= elem (0, j+1);
+	    }
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (nr, 1);
+	  Complex prod = elem (0, 0);
+	  for (int i = 0; i < nr; i++)
+	    {
+	      retval.elem (i, 0) = prod;
+	      if (i < nr - 1)
+		prod *= elem (i+1, 0);
+	    }
+	}
+      else
+	{
+	  retval.resize (nr, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      Complex prod = elem (0, j);
+	      for (int i = 0; i < nr; i++)
+		{
+		  retval.elem (i, j) = prod;
+		  if (i < nr - 1)
+		    prod *= elem (i+1, j);
+		}
+	    }
+	}
+    }
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::cumsum (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, nc);
+	  Complex sum = elem (0, 0);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = sum;
+	      if (j < nc - 1)
+		sum += elem (0, j+1);
+	    }
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (nr, 1);
+	  Complex sum = elem (0, 0);
+	  for (int i = 0; i < nr; i++)
+	    {
+	      retval.elem (i, 0) = sum;
+	      if (i < nr - 1)
+		sum += elem (i+1, 0);
+	    }
+	}
+      else
+	{
+	  retval.resize (nr, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      Complex sum = elem (0, j);
+	      for (int i = 0; i < nr; i++)
+		{
+		  retval.elem (i, j) = sum;
+		  if (i < nr - 1)
+		    sum += elem (i+1, j);
+		}
+	    }
+	}
+    }
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::prod (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 1.0;
+	  for (int j = 0; j < nc; j++)
+	    retval.elem (0, 0) *= elem (0, j);
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 1.0;
+	  for (int i = 0; i < nr; i++)
+	    retval.elem (0, 0) *= elem (i, 0);
+	}
+      else
+	{
+	  retval.resize (1, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = 1.0;
+	      for (int i = 0; i < nr; i++)
+		retval.elem (0, j) *= elem (i, j);
+	    }
+	}
+    }
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::sum (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 0.0;
+	  for (int j = 0; j < nc; j++)
+	    retval.elem (0, 0) += elem (0, j);
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 0.0;
+	  for (int i = 0; i < nr; i++)
+	    retval.elem (0, 0) += elem (i, 0);
+	}
+      else
+	{
+	  retval.resize (1, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = 0.0;
+	      for (int i = 0; i < nr; i++)
+		retval.elem (0, j) += elem (i, j);
+	    }
+	}
+    }
+  return retval;
+}
+
+ComplexMatrix
+ComplexMatrix::sumsq (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+  ComplexMatrix retval;
+  if (nr > 0 && nc > 0)
+    {
+      if (nr == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 0.0;
+	  for (int j = 0; j < nc; j++)
+	    {
+	      Complex d = elem (0, j);
+	      retval.elem (0, 0) += d * d;
+	    }
+	}
+      else if (nc == 1)
+	{
+	  retval.resize (1, 1);
+	  retval.elem (0, 0) = 0.0;
+	  for (int i = 0; i < nr; i++)
+	    {
+	      Complex d = elem (i, 0);
+	      retval.elem (0, 0) += d * d;
+	    }
+	}
+      else
+	{
+	  retval.resize (1, nc);
+	  for (int j = 0; j < nc; j++)
+	    {
+	      retval.elem (0, j) = 0.0;
+	      for (int i = 0; i < nr; i++)
+		{
+		  Complex d = elem (i, j);
+		  retval.elem (0, j) += d * d;
+		}
+	    }
+	}
+    }
+  return retval;
+}
+
+ComplexColumnVector
+ComplexMatrix::diag (void) const
+{
+  return diag (0);
+}
+
+ComplexColumnVector
+ComplexMatrix::diag (int k) const
+{
+  int nnr = rows ();
+  int nnc = cols ();
+  if (k > 0)
+    nnc -= k;
+  else if (k < 0)
+    nnr += k;
+
+  ComplexColumnVector d;
+
+  if (nnr > 0 && nnc > 0)
+    {
+      int ndiag = (nnr < nnc) ? nnr : nnc;
+
+      d.resize (ndiag);
+
+      if (k > 0)
+	{
+	  for (int i = 0; i < ndiag; i++)
+	    d.elem (i) = elem (i, i+k);
+	}
+      else if ( k < 0)
+	{
+	  for (int i = 0; i < ndiag; i++)
+	    d.elem (i) = elem (i-k, i);
+	}
+      else
+	{
+	  for (int i = 0; i < ndiag; i++)
+	    d.elem (i) = elem (i, i);
+	}
+    }
+  else
+    cerr << "diag: requested diagonal out of range\n";
+
+  return d;
+}
+
+ComplexColumnVector
+ComplexMatrix::row_min (void) const
+{
+  ComplexColumnVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nr);
+
+      for (int i = 0; i < nr; i++)
+	{
+	  Complex res = elem (i, 0);
+	  double absres = abs (res);
+	  for (int j = 1; j < nc; j++)
+	    if (abs (elem (i, j)) < absres)
+	      {
+		res = elem (i, j);
+		absres = abs (res);
+	      }
+	  result.elem (i) = res;
+	}
+    }
+
+  return result;
+}
+
+ComplexColumnVector
+ComplexMatrix::row_min_loc (void) const
+{
+  ComplexColumnVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nr);
+
+      for (int i = 0; i < nr; i++)
+        {
+          Complex res = 0;
+          double absres = abs (elem (i, 0));
+          for (int j = 0; j < nc; j++)
+            if (abs (elem (i, j)) < absres)
+              {
+                res = j;
+                absres = abs (elem (i, j));
+              }
+          result.elem (i) = res + 1;
+        }
+    }
+
+  return result;
+}
+
+ComplexColumnVector
+ComplexMatrix::row_max (void) const
+{
+  ComplexColumnVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nr);
+
+      for (int i = 0; i < nr; i++)
+	{
+	  Complex res = elem (i, 0);
+	  double absres = abs (res);
+	  for (int j = 1; j < nc; j++)
+	    if (abs (elem (i, j)) > absres)
+	      {
+		res = elem (i, j);
+		absres = abs (res);
+	      }
+	  result.elem (i) = res;
+	}
+    }
+
+  return result;
+}
+
+ComplexColumnVector
+ComplexMatrix::row_max_loc (void) const
+{
+  ComplexColumnVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nr);
+
+      for (int i = 0; i < nr; i++)
+        {
+          Complex res = 0;
+          double absres = abs (elem (i, 0));
+          for (int j = 0; j < nc; j++)
+            if (abs (elem (i, j)) > absres)
+              {
+                res = j;
+                absres = abs (elem (i, j));
+              }
+          result.elem (i) = res + 1;
+        }
+    }
+
+  return result;
+}
+
+ComplexRowVector
+ComplexMatrix::column_min (void) const
+{
+  ComplexRowVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nc);
+
+      for (int j = 0; j < nc; j++)
+	{
+	  Complex res = elem (0, j);
+	  double absres = abs (res);
+	  for (int i = 1; i < nr; i++)
+	    if (abs (elem (i, j)) < absres)
+	      {
+		res = elem (i, j);
+		absres = abs (res);
+	      }
+	  result.elem (j) = res;
+	}
+    }
+
+  return result;
+}
+
+ComplexRowVector
+ComplexMatrix::column_min_loc (void) const
+{
+  ComplexRowVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nc);
+
+      for (int j = 0; j < nc; j++)
+        {
+          Complex res = 0;
+          double absres = abs (elem (0, j));
+          for (int i = 0; i < nr; i++)
+            if (abs (elem (i, j)) < absres)
+              {
+                res = i;
+                absres = abs (elem (i, j));
+              }
+          result.elem (j) = res + 1;
+        }
+    }
+
+  return result;
+}
+
+ComplexRowVector
+ComplexMatrix::column_max (void) const
+{
+  ComplexRowVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nc);
+
+      for (int j = 0; j < nc; j++)
+	{
+	  Complex res = elem (0, j);
+	  double absres = abs (res);
+	  for (int i = 1; i < nr; i++)
+	    if (abs (elem (i, j)) > absres)
+	      {
+		res = elem (i, j);
+		absres = abs (res);
+	      }
+	  result.elem (j) = res;
+	}
+    }
+
+  return result;
+}
+
+ComplexRowVector
+ComplexMatrix::column_max_loc (void) const
+{
+  ComplexRowVector result;
+
+  int nr = rows ();
+  int nc = cols ();
+
+  if (nr > 0 && nc > 0)
+    {
+      result.resize (nc);
+
+      for (int j = 0; j < nc; j++)
+        {
+          Complex res = 0;
+          double absres = abs (elem (0, j));
+          for (int i = 0; i < nr; i++)
+            if (abs (elem (i, j)) > absres)
+              {
+                res = i;
+                absres = abs (elem (i, j));
+              }
+          result.elem (j) = res + 1;
+        }
+    }
+
+  return result;
+}
+
+// i/o
+
+ostream&
+operator << (ostream& os, const ComplexMatrix& a)
+{
+//  int field_width = os.precision () + 7;
+  for (int i = 0; i < a.rows (); i++)
+    {
+      for (int j = 0; j < a.cols (); j++)
+	os << " " /* setw (field_width) */ << a.elem (i, j);
+      os << "\n";
+    }
+  return os;
+}
+
+istream&
+operator >> (istream& is, ComplexMatrix& a)
+{
+  int nr = a.rows ();
+  int nc = a.cols ();
+
+  if (nr < 1 || nc < 1)
+    is.clear (ios::badbit);
+  else
+    {
+      Complex tmp;
+      for (int i = 0; i < nr; i++)
+	for (int j = 0; j < nc; j++)
+	  {
+	    is >> tmp;
+	    if (is)
+	      a.elem (i, j) = tmp;
+	    else
+	      break;
+	  }
+    }
+
+  return is;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/