deleted file mode 100644
--- a/liboctave/ColVector.cc
+++ /dev/null
@@ -1,1092 +0,0 @@
-// 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
-
-#include <iostream.h>
-
-#include "Matrix.h"
-#include "mx-inlines.cc"
-#include "f77-uscore.h"
-#include "lo-error.h"
-
-// Fortran functions we call.
-
-extern "C"
-{
-  int F77_FCN (dgemm) (const char*, const char*, const int*,
-		       const int*, const int*, const double*,
-		       const double*, const int*, const double*,
-		       const int*, const double*, double*, const int*,
-		       long, long);
-
-/*
- * f2c translates complex*16 as
- *
- *   typedef struct { doublereal re, im; } doublecomplex;
- *
- * and Complex.h from libg++ uses
- *
- *   protected:
- *     double re;
- *     double im;
- *
- * as the only data members, so this should work (fingers crossed that
- * things don't change).
- */
-
-  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);
-}
-
-/*
- * Column Vector class.
- */
-
-#if 0
-ColumnVector&
-ColumnVector::resize (int n)
-{
-  if (n < 0)
-    {
-      (*current_liboctave_error_handler)
-	("can't resize to negative dimension");
-      return *this;
-    }
-
-  double *new_data = (double *) NULL;
-  if (n > 0)
-    {
-      new_data = new double [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;
-}
-
-ColumnVector&
-ColumnVector::resize (int n, double val)
-{
-  int old_len = len;
-  resize (n);
-  for (int i = old_len; i < len; i++)
-    data[i] = val;
-
-  return *this;
-}
-#endif
-
-int
-ColumnVector::operator == (const ColumnVector& a) const
-{
-  int len = length ();
-  if (len != a.length ())
-    return 0;
-  return equal (data (), a.data (), len);
-}
-
-int
-ColumnVector::operator != (const ColumnVector& a) const
-{
-  return !(*this == a);
-}
-
-ColumnVector&
-ColumnVector::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;
-}
-
-ColumnVector&
-ColumnVector::fill (double val)
-{
-  int len = length ();
-  if (len > 0)
-    for (int i = 0; i < len; i++)
-      elem (i) = val;
-  return *this;
-}
-
-ColumnVector&
-ColumnVector::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;
-}
-
-ColumnVector
-ColumnVector::stack (const ColumnVector& a) const
-{
-  int len = length ();
-  int nr_insert = len;
-  ColumnVector retval (len + a.length ());
-  retval.insert (*this, 0);
-  retval.insert (a, nr_insert);
-  return retval;
-}
-
-RowVector
-ColumnVector::transpose (void) const
-{
-  int len = length ();
-  return RowVector (dup (data (), len), len);
-}
-
-// resize is the destructive equivalent for this one
-
-ColumnVector
-ColumnVector::extract (int r1, int r2) const
-{
-  if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
-
-  int new_r = r2 - r1 + 1;
-
-  ColumnVector 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
-
-ColumnVector&
-ColumnVector::operator += (const ColumnVector& a)
-{
-  int len = length ();
-  if (len != a.length ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant vector += operation attempted");
-      return ColumnVector ();
-    }
-
-  if (len == 0)
-    return *this;
-
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  add2 (d, a.data (), len);
-  return *this;
-}
-
-ColumnVector&
-ColumnVector::operator -= (const ColumnVector& a)
-{
-  int len = length ();
-  if (len != a.length ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant vector -= operation attempted");
-      return ColumnVector ();
-    }
-
-  if (len == 0)
-    return *this;
-
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  subtract2 (d, a.data (), len);
-  return *this;
-}
-
-// scalar by column vector -> column vector operations
-
-ComplexColumnVector
-operator + (const ColumnVector& a, const Complex& s)
-{
-  int len = a.length ();
-  return ComplexColumnVector (add (a.data (), len, s), len);
-}
-
-ComplexColumnVector
-operator - (const ColumnVector& a, const Complex& s)
-{
-  int len = a.length ();
-  return ComplexColumnVector (subtract (a.data (), len, s), len);
-}
-
-ComplexColumnVector
-operator * (const ColumnVector& a, const Complex& s)
-{
-  int len = a.length ();
-  return ComplexColumnVector (multiply (a.data (), len, s), len);
-}
-
-ComplexColumnVector
-operator / (const ColumnVector& a, const Complex& s)
-{
-  int len = a.length ();
-  return ComplexColumnVector (divide (a.data (), len, s), len);
-}
-
-// scalar by column vector -> column vector operations
-
-ComplexColumnVector
-operator + (const Complex& s, const ColumnVector& a)
-{
-  int a_len = a.length ();
-  return ComplexColumnVector (add (a.data (), a_len, s), a_len);
-}
-
-ComplexColumnVector
-operator - (const Complex& s, const ColumnVector& a)
-{
-  int a_len = a.length ();
-  return ComplexColumnVector (subtract (s, a.data (), a_len), a_len);
-}
-
-ComplexColumnVector
-operator * (const Complex& s, const ColumnVector& a)
-{
-  int a_len = a.length ();
-  return ComplexColumnVector (multiply (a.data (), a_len, s), a_len);
-}
-
-ComplexColumnVector
-operator / (const Complex& s, const ColumnVector& a)
-{
-  int a_len = a.length ();
-  return ComplexColumnVector (divide (s, a.data (), a_len), a_len);
-}
-
-// column vector by row vector -> matrix operations
-
-Matrix
-operator * (const ColumnVector& v, const RowVector& a)
-{
-  int len = v.length ();
-  int a_len = a.length ();
-  if (len != a_len)
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant vector multiplication attempted");
-      return Matrix ();
-    }
-
-  if (len == 0)
-    return Matrix (len, len, 0.0);
-
-  char transa = 'N';
-  char transb = 'N';
-  double alpha = 1.0;
-  double beta  = 0.0;
-  int anr = 1;
-
-  double *c = new double [len * a_len];
-
-  F77_FCN (dgemm) (&transa, &transb, &len, &a_len, &anr, &alpha,
-		   v.data (), &len, a.data (), &anr, &beta, c, &len,
-		   1L, 1L); 
-
-  return Matrix (c, len, a_len);
-}
-
-ComplexMatrix
-operator * (const ColumnVector& v, const ComplexRowVector& a)
-{
-  ComplexColumnVector tmp (v);
-  return tmp * a;
-}
-
-ComplexColumnVector
-operator + (const ColumnVector& v, const ComplexColumnVector& 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 (add (v.data (), a.data (), len), len);
-}
-
-ComplexColumnVector
-operator - (const ColumnVector& v, const ComplexColumnVector& 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 ColumnVector& v, const ComplexColumnVector& a)
-{
-  int len = v.length ();
-  if (len != a.length ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant vector product attempted");
-      return ColumnVector ();
-    }
-
-  if (len == 0)
-    return ComplexColumnVector (0);
-
-  return ComplexColumnVector (multiply (v.data (), a.data (), len), len);
-}
-
-ComplexColumnVector
-quotient (const ColumnVector& v, const ComplexColumnVector& a)
-{
-  int len = v.length ();
-  if (len != a.length ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant vector quotient attempted");
-      return ColumnVector ();
-    }
-
-  if (len == 0)
-    return ComplexColumnVector (0);
-
-  return ComplexColumnVector (divide (v.data (), a.data (), len), len);
-}
-
-// other operations
-
-ColumnVector
-map (d_d_Mapper f, const ColumnVector& a)
-{
-  ColumnVector b (a);
-  b.map (f);
-  return b;
-}
-
-void
-ColumnVector::map (d_d_Mapper f)
-{
-  for (int i = 0; i < length (); i++)
-    elem (i) = f (elem (i));
-}
-
-double
-ColumnVector::min (void) const
-{
-  int len = length ();
-  if (len == 0)
-    return 0.0;
-
-  double res = elem (0);
-
-  for (int i = 1; i < len; i++)
-    if (elem (i) < res)
-      res = elem (i);
-
-  return res;
-}
-
-double
-ColumnVector::max (void) const
-{
-  int len = length ();
-  if (len == 0)
-    return 0.0;
-
-  double res = elem (0);
-
-  for (int i = 1; i < len; i++)
-    if (elem (i) > res)
-      res = elem (i);
-
-  return res;
-}
-
-ostream&
-operator << (ostream& os, const ColumnVector& 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, ColumnVector& 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;
-        }
-    }
-}
-
-/*
- * Complex Column Vector class
- */
-
-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: ***
-*/

deleted file mode 100644
--- a/liboctave/DiagMatrix.cc
+++ /dev/null
@@ -1,2044 +0,0 @@
-// 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
-
-#include <iostream.h>
-
-#include "Matrix.h"
-#include "mx-inlines.cc"
-#include "lo-error.h"
-
-/*
- * Diagonal Matrix class.
- */
-
-#if 0
-DiagMatrix&
-DiagMatrix::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;
-  double *new_data = (double *) NULL;
-  if (new_len > 0)
-    {
-      new_data = new double [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;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-  double *new_data = (double *) NULL;
-  if (new_len > 0)
-    {
-      new_data = new double [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
-DiagMatrix::operator == (const DiagMatrix& a) const
-{
-  if (rows () != a.rows () || cols () != a.cols ())
-    return 0;
-
-  return equal (data (), a.data (), length ());
-}
-
-int
-DiagMatrix::operator != (const DiagMatrix& a) const
-{
-  return !(*this == a);
-}
-
-DiagMatrix&
-DiagMatrix::fill (double val)
-{
-  for (int i = 0; i < length (); i++)
-    elem (i, i) = val;
-  return *this;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-}
-
-DiagMatrix
-DiagMatrix::transpose (void) const
-{
-  return DiagMatrix (dup (data (), length ()), cols (), rows ());
-}
-
-Matrix
-DiagMatrix::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;
-
-  Matrix 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.
-
-RowVector
-DiagMatrix::row (int i) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  if (i < 0 || i >= nr)
-    {
-      (*current_liboctave_error_handler) ("invalid row selection");
-      return RowVector (); 
-    }
-
-  RowVector retval (nc, 0.0);
-  if (nr <= nc || (nr > nc && i < nc))
-    retval.elem (i) = elem (i, i);
-
-  return retval;
-}
-
-RowVector
-DiagMatrix::row (char *s) const
-{
-  if (s == (char *) NULL)
-    {
-      (*current_liboctave_error_handler) ("invalid row selection");
-      return RowVector (); 
-    }
-
-  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 RowVector (); 
-    }
-}
-
-ColumnVector
-DiagMatrix::column (int i) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  if (i < 0 || i >= nc)
-    {
-      (*current_liboctave_error_handler) ("invalid column selection");
-      return ColumnVector (); 
-    }
-
-  ColumnVector retval (nr, 0.0);
-  if (nr >= nc || (nr < nc && i < nr))
-    retval.elem (i) = elem (i, i);
-
-  return retval;
-}
-
-ColumnVector
-DiagMatrix::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 (); 
-    }
-}
-
-DiagMatrix
-DiagMatrix::inverse (void) const
-{
-  int info;
-  return inverse (info);
-}
-
-DiagMatrix
-DiagMatrix::inverse (int &info) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  int len = length ();
-  if (nr != nc)
-    {
-      (*current_liboctave_error_handler) ("inverse requires square matrix");
-      return DiagMatrix ();
-    }
-
-  info = 0;
-  double *tmp_data = dup (data (), len);
-  for (int i = 0; i < len; i++)
-    {
-      if (elem (i, i) == 0.0)
-	{
-	  info = -1;
-	  copy (tmp_data, data (), len); // Restore contents.
-	  break;
-	}
-      else
-	{
-	  tmp_data[i] = 1.0 / elem (i, i);
-	}
-    }
-
-  return DiagMatrix (tmp_data, nr, nc);
-}
-
-// diagonal matrix by diagonal matrix -> diagonal matrix operations
-
-DiagMatrix&
-DiagMatrix::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 (nc == 0 || nr == 0)
-    return *this;
-
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  add2 (d, a.data (), length ());
-  return *this;
-}
-
-DiagMatrix&
-DiagMatrix::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;
-
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  subtract2 (d, a.data (), length ());
-  return *this;
-}
-
-// diagonal matrix by scalar -> matrix operations
-
-Matrix
-operator + (const DiagMatrix& a, double s)
-{
-  Matrix tmp (a.rows (), a.cols (), s);
-  return a + tmp;
-}
-
-Matrix
-operator - (const DiagMatrix& a, double s)
-{
-  Matrix tmp (a.rows (), a.cols (), -s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator + (const DiagMatrix& a, const Complex& s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return a + tmp;
-}
-
-ComplexMatrix
-operator - (const DiagMatrix& a, const Complex& s)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), -s);
-  return a + tmp;
-}
-
-// diagonal matrix by scalar -> diagonal matrix operations
-
-ComplexDiagMatrix
-operator * (const DiagMatrix& a, const Complex& s)
-{
-  return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
-			    a.rows (), a.cols ());
-}
-
-ComplexDiagMatrix
-operator / (const DiagMatrix& a, const Complex& s)
-{
-  return ComplexDiagMatrix (divide (a.data (), a.length (), s),
-			    a.rows (), a.cols ());
-}
-
-// scalar by diagonal matrix -> matrix operations
-
-Matrix
-operator + (double s, const DiagMatrix& a)
-{
-  Matrix tmp (a.rows (), a.cols (), s);
-  return tmp + a;
-}
-
-Matrix
-operator - (double s, const DiagMatrix& a)
-{
-  Matrix tmp (a.rows (), a.cols (), s);
-  return tmp - a;
-}
-
-ComplexMatrix
-operator + (const Complex& s, const DiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp + a;
-}
-
-ComplexMatrix
-operator - (const Complex& s, const DiagMatrix& a)
-{
-  ComplexMatrix tmp (a.rows (), a.cols (), s);
-  return tmp - a;
-}
-
-// scalar by diagonal matrix -> diagonal matrix operations
-
-ComplexDiagMatrix
-operator * (const Complex& s, const DiagMatrix& a)
-{
-  return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
-			    a.rows (), a.cols ());
-}
-
-// diagonal matrix by column vector -> column vector operations
-
-ColumnVector
-operator * (const DiagMatrix& 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 multiplication attempted");
-      return ColumnVector ();
-    }
-
-  if (nc == 0 || nr == 0)
-    return ColumnVector (0);
-
-  ColumnVector 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 DiagMatrix& 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 multiplication attempted");
-      return ColumnVector ();
-    }
-
-  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
-
-DiagMatrix
-operator * (const DiagMatrix& 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 DiagMatrix ();
-    }
-
-  if (nr_a == 0 || nc_a == 0 || nc_b == 0)
-    return DiagMatrix (nr_a, nc_a, 0.0);
-
-  DiagMatrix c (nr_a, nc_b);
-
-  int len = nr_a < nc_b ? nr_a : nc_b;
-
-  for (int i = 0; i < len; i++)
-    {
-      double 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
-operator + (const DiagMatrix& 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 ComplexDiagMatrix ();
-    }
-
-  if (nc == 0 || nr == 0)
-    return ComplexDiagMatrix (nr, nc);
-
-  return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()),  nr, nc);
-}
-
-ComplexDiagMatrix
-operator - (const DiagMatrix& 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 ComplexDiagMatrix ();
-    }
-
-  if (nc == 0 || nr == 0)
-    return ComplexDiagMatrix (nr, nc);
-
-  return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()),
-			    nr, nc);
-}
-
-ComplexDiagMatrix
-operator * (const DiagMatrix& 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++)
-    {
-      double 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
-product (const DiagMatrix& 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 product attempted");
-      return ComplexDiagMatrix ();
-    }
-
-  if (nc == 0 || nr == 0)
-    return ComplexDiagMatrix (nr, nc);
-
-  return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
-			    nr, nc);
-}
-
-// diagonal matrix by matrix -> matrix operations
-
-Matrix
-operator + (const DiagMatrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return Matrix (nr, nc);
-
-  Matrix result (a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-Matrix
-operator - (const DiagMatrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return Matrix (nr, nc);
-
-  Matrix result (-a);
-  for (int i = 0; i < m.length (); i++)
-    result.elem (i, i) += m.elem (i, i);
-
-  return result;
-}
-
-Matrix
-operator * (const DiagMatrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return Matrix (nr, a_nc, 0.0);
-
-  Matrix 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 DiagMatrix& 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 DiagMatrix& 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 DiagMatrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-  int a_nr = a.rows ();
-  int a_nc = a.cols ();
-  if (nc != a_nr)
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant matrix multiplication attempted");
-      return ComplexMatrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return ComplexMatrix (nr, nc, 0.0);
-
-  ComplexMatrix c (nr, a_nc);
-
-  for (int i = 0; i < m.length (); i++)
-    {
-      if (m.elem (i, i) == 1.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = a.elem (i, j);
-	}
-      else if (m.elem (i, i) == 0.0)
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = 0.0;
-	}
-      else
-	{
-	  for (int j = 0; j < a_nc; j++)
-	    c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
-	}
-    }
-
-  if (nr > nc)
-    {
-      for (int j = 0; j < a_nc; j++)
-	for (int i = a_nr; i < nr; i++)
-	  c.elem (i, j) = 0.0;
-    }
-
-  return c;
-}
-
-// other operations
-
-ColumnVector
-DiagMatrix::diag (void) const
-{
-  return diag (0);
-}
-
-// Could be optimized...
-
-ColumnVector
-DiagMatrix::diag (int k) const
-{
-  int nnr = rows ();
-  int nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  ColumnVector 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;
-}
-
-ostream&
-operator << (ostream& os, const DiagMatrix& a)
-{
-//  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) */ << 0.0;
-	}
-      os << "\n";
-    }
-  return os;
-}
-
-/*
- * Complex Diagonal Matrix class
- */
-
-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: ***
-*/

deleted file mode 100644
--- a/liboctave/Matrix.cc
+++ /dev/null
@@ -1,5032 +0,0 @@
-// 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
-
-#include <sys/types.h>
-#include <iostream.h>
-
-#include "Matrix.h"
-#include "mx-inlines.cc"
-#include "lo-error.h"
-#include "f77-uscore.h"
-
-// Fortran functions we call.
-
-extern "C"
-{
-  int F77_FCN (dgemm) (const char*, const char*, const int*,
-		       const int*, const int*, const double*,
-		       const double*, const int*, const double*,
-		       const int*, const double*, double*, const int*,
-		       long, long);
-
-  int F77_FCN (dgemv) (const char*, const int*, const int*,
-		       const double*, const double*, const int*,
-		       const double*, const int*, const double*,
-		       double*, const int*, long);
-
-  int F77_FCN (dgeco) (double*, const int*, const int*, int*, double*,
-		       double*);
-
-  int F77_FCN (dgesl) (const double*, const int*, const int*,
-		       const int*, double*, const int*); 
-
-  int F77_FCN (dgedi) (double*, const int*, const int*, const int*,
-		       double*, double*, const int*);
-
-  int F77_FCN (dgelss) (const int*, const int*, const int*, double*,
-			const int*, double*, const int*, double*,
-			const double*, int*, double*, const int*,
-			int*);
-
-/*
- * f2c translates complex*16 as
- *
- *   typedef struct { doublereal re, im; } doublecomplex;
- *
- * and Complex.h from libg++ uses
- *
- *   protected:
- *     double re;
- *     double im;
- *
- * as the only data members, so this should work (fingers crossed that
- * things don't change).
- */
-
-  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*);
-}
-
-// Since this is only temporary, put all of this here, rather than
-// putting each type where it logically belongs.  This way, it will be
-// easier to delete.
-
-#define KLUDGE_MATRICES
-#define TYPE double
-#define KL_MAT_TYPE Matrix
-#include "mx-kludge.cc"
-#undef KLUDGE_MATRICES
-#undef TYPE
-#undef KL_MAT_TYPE
-
-#define KLUDGE_VECTORS
-#define TYPE double
-#define KL_VEC_TYPE ColumnVector
-#include "mx-kludge.cc"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-#define KLUDGE_VECTORS
-#define TYPE double
-#define KL_VEC_TYPE RowVector
-#include "mx-kludge.cc"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-#define KLUDGE_DIAG_MATRICES
-#define TYPE double
-#define KL_DMAT_TYPE DiagMatrix
-#include "mx-kludge.cc"
-#undef KLUDGE_DIAG_MATRICES
-#undef TYPE
-#undef KL_DMAT_TYPE
-
-#define KLUDGE_MATRICES
-#define TYPE Complex
-#define KL_MAT_TYPE ComplexMatrix
-#include "mx-kludge.cc"
-#undef KLUDGE_MATRICES
-#undef TYPE
-#undef KL_MAT_TYPE
-
-#define KLUDGE_VECTORS
-#define TYPE Complex
-#define KL_VEC_TYPE ComplexColumnVector
-#include "mx-kludge.cc"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-#define KLUDGE_VECTORS
-#define TYPE Complex
-#define KL_VEC_TYPE ComplexRowVector
-#include "mx-kludge.cc"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-#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
-
-/*
- * Matrix class.
- */
-
-Matrix::Matrix (const DiagMatrix& a)
-  : Array2<double> (a.rows (), a.cols (), 0.0)
-{
-  for (int i = 0; i < a.length (); i++)
-    elem (i, i) = a.elem (i, i);
-}
-
-#if 0
-Matrix&
-Matrix::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;
-  double* new_data = (double *) NULL;
-  if (new_len > 0)
-    {
-      new_data = new double [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;
-}
-
-Matrix&
-Matrix::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;
-  double *new_data = (double *) NULL;
-  if (new_len > 0)
-    {
-      new_data = new double [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
-Matrix::operator == (const Matrix& a) const
-{
-  if (rows () != a.rows () || cols () != a.cols ())
-    return 0;
-
-  return equal (data (), a.data (), length ());
-}
-
-int
-Matrix::operator != (const Matrix& a) const
-{
-  return !(*this == a);
-}
-
-Matrix&
-Matrix::insert (const Matrix& a, int r, int c)
-{
-  int a_rows = a.rows ();
-  int a_cols = a.cols ();
-  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 j = 0; j < a_cols; j++)
-    for (int i = 0; i < a_rows; i++)
-      elem (r+i, c+j) = a.elem (i, j);
-
-  return *this;
-}
-
-Matrix&
-Matrix::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;
-}
-
-Matrix&
-Matrix::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;
-}
-
-Matrix&
-Matrix::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;
-}
-
-Matrix&
-Matrix::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;
-}
-
-Matrix&
-Matrix::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;
-}
-
-Matrix
-Matrix::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 Matrix ();
-    }
-
-  int nc_insert = nc;
-  Matrix retval (nr, nc + a.cols ());
-  retval.insert (*this, 0, 0);
-  retval.insert (a, 0, nc_insert);
-  return retval;
-}
-
-Matrix
-Matrix::append (const RowVector& a) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  if (nr != 1)
-    {
-      (*current_liboctave_error_handler) ("row dimension mismatch for append");
-      return Matrix ();
-    }
-
-  int nc_insert = nc;
-  Matrix retval (nr, nc + a.length ());
-  retval.insert (*this, 0, 0);
-  retval.insert (a, 0, nc_insert);
-  return retval;
-}
-
-Matrix
-Matrix::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 Matrix ();
-    }
-
-  int nc_insert = nc;
-  Matrix retval (nr, nc + 1);
-  retval.insert (*this, 0, 0);
-  retval.insert (a, 0, nc_insert);
-  return retval;
-}
-
-Matrix
-Matrix::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;
-  Matrix retval (nr, nc + a.cols ());
-  retval.insert (*this, 0, 0);
-  retval.insert (a, 0, nc_insert);
-  return retval;
-}
-
-Matrix
-Matrix::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 Matrix ();
-    }
-
-  int nr_insert = nr;
-  Matrix retval (nr + a.rows (), nc);
-  retval.insert (*this, 0, 0);
-  retval.insert (a, nr_insert, 0);
-  return retval;
-}
-
-Matrix
-Matrix::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 Matrix ();
-    }
-
-  int nr_insert = nr;
-  Matrix retval (nr + 1, nc);
-  retval.insert (*this, 0, 0);
-  retval.insert (a, nr_insert, 0);
-  return retval;
-}
-
-Matrix
-Matrix::stack (const ColumnVector& a) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  if (nc != 1)
-    {
-      (*current_liboctave_error_handler)
-	("column dimension mismatch for stack");
-      return Matrix ();
-    }
-
-  int nr_insert = nr;
-  Matrix retval (nr + a.length (), nc);
-  retval.insert (*this, 0, 0);
-  retval.insert (a, nr_insert, 0);
-  return retval;
-}
-
-Matrix
-Matrix::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 Matrix ();
-    }
-
-  int nr_insert = nr;
-  Matrix retval (nr + a.rows (), nc);
-  retval.insert (*this, 0, 0);
-  retval.insert (a, nr_insert, 0);
-  return retval;
-}
-
-Matrix
-Matrix::transpose (void) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  Matrix 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
-Matrix::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;
-
-  Matrix 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.
-
-RowVector
-Matrix::row (int i) const
-{
-  int nc = cols ();
-  if (i < 0 || i >= rows ())
-    {
-      (*current_liboctave_error_handler) ("invalid row selection");
-      return RowVector ();
-    }
-
-  RowVector retval (nc);
-  for (int j = 0; j < nc; j++)
-    retval.elem (j) = elem (i, j);
-
-  return retval;
-}
-
-RowVector
-Matrix::row (char *s) const
-{
-  if (s == (char *) NULL)
-    {
-      (*current_liboctave_error_handler) ("invalid row selection");
-      return RowVector ();
-    }
-
-  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 RowVector ();
-    }
-}
-
-ColumnVector
-Matrix::column (int i) const
-{
-  int nr = rows ();
-  if (i < 0 || i >= cols ())
-    {
-      (*current_liboctave_error_handler) ("invalid column selection");
-      return ColumnVector ();
-    }
-
-  ColumnVector retval (nr);
-  for (int j = 0; j < nr; j++)
-    retval.elem (j) = elem (j, i);
-
-  return retval;
-}
-
-ColumnVector
-Matrix::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 ();
-    }
-}
-
-Matrix
-Matrix::inverse (void) const
-{
-  int info;
-  double rcond;
-  return inverse (info, rcond);
-}
-
-Matrix
-Matrix::inverse (int& info) const
-{
-  double rcond;
-  return inverse (info, rcond);
-}
-
-Matrix
-Matrix::inverse (int& info, double& rcond) const
-{
-  int nr = rows ();
-  int nc = cols ();
-  int len = length ();
-  if (nr != nc || nr == 0 || nc == 0)
-    {
-      (*current_liboctave_error_handler) ("inverse requires square matrix");
-      return Matrix ();
-    }
-
-  info = 0;
-
-  int *ipvt = new int [nr];
-  double *z = new double [nr];
-  double *tmp_data = dup (data (), len);
-
-  F77_FCN (dgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z);
-
-  if (rcond + 1.0 == 1.0)
-    {
-      info = -1;
-      copy (tmp_data, data (), len);  // Restore matrix contents.
-    }
-  else
-    {
-      int job = 1;
-      double dummy;
-
-      F77_FCN (dgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job);
-    }
-
-  delete [] ipvt;
-  delete [] z;
-
-  return Matrix (tmp_data, nr, nc);
-}
-
-ComplexMatrix
-Matrix::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 = make_complex (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
-Matrix::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 = make_complex (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);
-}
-
-DET
-Matrix::determinant (void) const
-{
-  int info;
-  double rcond;
-  return determinant (info, rcond);
-}
-
-DET
-Matrix::determinant (int& info) const
-{
-  double rcond;
-  return determinant (info, rcond);
-}
-
-DET
-Matrix::determinant (int& info, double& rcond) const
-{
-  DET retval;
-
-  int nr = rows ();
-  int nc = cols ();
-  if (nr == 0 || nc == 0)
-    {
-      double d[2];
-      d[0] = 1.0;
-      d[1] = 0.0;
-      return DET (d);
-    }
-
-  info = 0;
-  int *ipvt = new int [nr];
-
-  double *z = new double [nr];
-  double *tmp_data = dup (data (), length ());
-
-  F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
-
-  if (rcond + 1.0 == 1.0)
-    {
-      info = -1;
-    }
-  else
-    {
-      int job = 10;
-      double d[2];
-      F77_FCN (dgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job);
-      retval = DET (d);
-    }
-
-  delete [] tmp_data;
-  delete [] ipvt;
-  delete [] z;
-
-  return retval;
-}
-
-Matrix
-Matrix::solve (const Matrix& b) const
-{
-  int info;
-  double rcond;
-  return solve (b, info, rcond);
-}
-
-Matrix
-Matrix::solve (const Matrix& b, int& info) const
-{
-  double rcond;
-  return solve (b, info, rcond);
-}
-
-Matrix
-Matrix::solve (const Matrix& b, int& info, double& rcond) const
-{
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-  if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
-    {
-      (*current_liboctave_error_handler)
-	("matrix dimension mismatch solution of linear equations");
-      return Matrix ();
-    }
-
-  info = 0;
-  int *ipvt = new int [nr];
-
-  double *z = new double [nr];
-  double *tmp_data = dup (data (), length ());
-
-  F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
-
-  if (rcond + 1.0 == 1.0)
-    {
-      info = -2;
-    }
-  else
-    {
-      int job = 0;
-
-      double *result = dup (b.data (), b.length ());
-
-      int b_nc = b.cols ();
-      for (int j = 0; j < b_nc; j++)
-	F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job);
-
-      retval = Matrix (result, b.rows (), b_nc);
-    }
-
-  delete [] tmp_data;
-  delete [] ipvt;
-  delete [] z;
-
-  return retval;
-}
-
-ComplexMatrix
-Matrix::solve (const ComplexMatrix& b) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.solve (b);
-}
-
-ComplexMatrix
-Matrix::solve (const ComplexMatrix& b, int& info) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.solve (b, info);
-}
-
-ComplexMatrix
-Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.solve (b, info, rcond);
-}
-
-ColumnVector
-Matrix::solve (const ColumnVector& b) const
-{
-  int info; double rcond;
-  return solve (b, info, rcond);
-}
-
-ColumnVector
-Matrix::solve (const ColumnVector& b, int& info) const
-{
-  double rcond;
-  return solve (b, info, rcond);
-}
-
-ColumnVector
-Matrix::solve (const ColumnVector& b, int& info, double& rcond) const
-{
-  ColumnVector retval;
-
-  int nr = rows ();
-  int nc = cols ();
-  if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
-    {
-      (*current_liboctave_error_handler)
-	("matrix dimension mismatch solution of linear equations");
-      return ColumnVector ();
-    }
-
-  info = 0;
-  int *ipvt = new int [nr];
-
-  double *z = new double [nr];
-  double *tmp_data = dup (data (), length ());
-
-  F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z);
-
-  if (rcond + 1.0 == 1.0)
-    {
-      info = -2;
-    }
-  else
-    {
-      int job = 0;
-
-      int b_len = b.length ();
-
-      double *result = dup (b.data (), b_len);
-
-      F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, result, &job);
-
-      retval = ColumnVector (result, b_len);
-    }
-
-  delete [] tmp_data;
-  delete [] ipvt;
-  delete [] z;
-
-  return retval;
-}
-
-ComplexColumnVector
-Matrix::solve (const ComplexColumnVector& b) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.solve (b);
-}
-
-ComplexColumnVector
-Matrix::solve (const ComplexColumnVector& b, int& info) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.solve (b, info);
-}
-
-ComplexColumnVector
-Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.solve (b, info, rcond);
-}
-
-Matrix
-Matrix::lssolve (const Matrix& b) const
-{
-  int info;
-  int rank;
-  return lssolve (b, info, rank);
-}
-
-Matrix
-Matrix::lssolve (const Matrix& b, int& info) const
-{
-  int rank;
-  return lssolve (b, info, rank);
-}
-
-Matrix
-Matrix::lssolve (const Matrix& 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 in solution of least squares problem");
-      return Matrix ();
-    }
-
-  double *tmp_data = dup (data (), length ());
-
-  int nrr = m > n ? m : n;
-  Matrix 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);
-
-  double *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 = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n));
-  else
-    lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m));
-
-  double *work = new double [lwork];
-
-  F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s,
-		    &rcond, &rank, work, &lwork, &info);
-
-  Matrix 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;
-
-  return retval;
-}
-
-ComplexMatrix
-Matrix::lssolve (const ComplexMatrix& b) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.lssolve (b);
-}
-
-ComplexMatrix
-Matrix::lssolve (const ComplexMatrix& b, int& info) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.lssolve (b);
-}
-
-ComplexMatrix
-Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.lssolve (b);
-}
-
-ColumnVector
-Matrix::lssolve (const ColumnVector& b) const
-{
-  int info;
-  int rank; return lssolve (b, info, rank);
-}
-
-ColumnVector
-Matrix::lssolve (const ColumnVector& b, int& info) const
-{
-  int rank;
-  return lssolve (b, info, rank);
-}
-
-ColumnVector
-Matrix::lssolve (const ColumnVector& 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 in solution of least squares problem");
-      return ColumnVector ();
-    }
-
-  double *tmp_data = dup (data (), length ());
-
-  int nrr = m > n ? m : n;
-  ColumnVector result (nrr);
-
-  int i;
-  for (i = 0; i < m; i++)
-    result.elem (i) = b.elem (i);
-
-  double *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 = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n));
-  else
-    lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m));
-
-  double *work = new double [lwork];
-
-  F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s,
-		    &rcond, &rank, work, &lwork, &info);
-
-  ColumnVector retval (n);
-  for (i = 0; i < n; i++)
-    retval.elem (i) = result.elem (i);
-
-  delete [] tmp_data;
-  delete [] s;
-  delete [] work;
-
-  return retval;
-}
-
-ComplexColumnVector
-Matrix::lssolve (const ComplexColumnVector& b) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.lssolve (b);
-}
-
-ComplexColumnVector
-Matrix::lssolve (const ComplexColumnVector& b, int& info) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.lssolve (b, info);
-}
-
-ComplexColumnVector
-Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const
-{
-  ComplexMatrix tmp (*this);
-  return tmp.lssolve (b, info, rank);
-}
-
-Matrix&
-Matrix::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;
-
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  add2 (d, a.data (), length ());
-
-  return *this;
-}
-
-Matrix&
-Matrix::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;
-
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  subtract2 (d, a.data (), length ());
-
-  return *this;
-}
-
-Matrix&
-Matrix::operator += (const DiagMatrix& a)
-{
-  if (rows () != a.rows () || cols () != a.cols ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant matrix += operation attempted");
-      return *this;
-    }
-
-  for (int i = 0; i < a.length (); i++)
-    elem (i, i) += a.elem (i, i);
-
-  return *this;
-}
-
-Matrix&
-Matrix::operator -= (const DiagMatrix& a)
-{
-  if (rows () != a.rows () || cols () != a.cols ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant matrix += operation attempted");
-      return *this;
-    }
-
-  for (int i = 0; i < a.length (); i++)
-    elem (i, i) -= a.elem (i, i);
-
-  return *this;
-}
-
-// unary operations
-
-Matrix
-Matrix::operator ! (void) const
-{
-  int nr = rows ();
-  int nc = cols ();
-
-  Matrix b (nr, nc);
-
-  for (int j = 0; j < nc; j++)
-    for (int i = 0; i < nr; i++)
-      b.elem (i, j) = ! elem (i, j);
-
-  return b;
-}
-
-// matrix by scalar -> matrix operations.
-
-ComplexMatrix
-operator + (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (add (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator - (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (subtract (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator * (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (multiply (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator / (const Matrix& a, const Complex& s)
-{
-  return ComplexMatrix (divide (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-// scalar by matrix -> matrix operations.
-
-ComplexMatrix
-operator + (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (add (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator - (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (subtract (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator * (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (multiply (a.data (), a.length (), s),
-			a.rows (), a.cols ());
-}
-
-ComplexMatrix
-operator / (const Complex& s, const Matrix& a)
-{
-  return ComplexMatrix (divide (s, a.data (), a.length ()),
-			a.rows (), a.cols ());
-}
-
-// matrix by column vector -> column vector operations
-
-ColumnVector
-operator * (const Matrix& m, const ColumnVector& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-  if (nc != a.length ())
-    {
-      (*current_liboctave_error_handler)
-	("nonconformant matrix multiplication attempted");
-      return ColumnVector ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return ColumnVector (0);
-
-  char trans = 'N';
-  int ld = nr;
-  double alpha = 1.0;
-  double beta  = 0.0;
-  int i_one = 1;
-
-  double *y = new double [nr];
-
-  F77_FCN (dgemv) (&trans, &nr, &nc, &alpha, m.data (), &ld, a.data (),
-		   &i_one, &beta, y, &i_one, 1L); 
-
-  return ColumnVector (y, nr);
-}
-
-ComplexColumnVector
-operator * (const Matrix& m, const ComplexColumnVector& a)
-{
-  ComplexMatrix tmp (m);
-  return tmp * a;
-}
-
-// matrix by diagonal matrix -> matrix operations
-
-Matrix
-operator + (const Matrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return Matrix (nr, nc);
-
-  Matrix result (m);
-  int a_len = a.length ();
-  for (int i = 0; i < a_len; i++)
-    result.elem (i, i) += a.elem (i, i);
-
-  return result;
-}
-
-Matrix
-operator - (const Matrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0)
-    return Matrix (nr, nc);
-
-  Matrix result (m);
-  int a_len = a.length ();
-  for (int i = 0; i < a_len; i++)
-    result.elem (i, i) -= a.elem (i, i);
-
-  return result;
-}
-
-Matrix
-operator * (const Matrix& m, const DiagMatrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return Matrix (nr, a_nc, 0.0);
-
-  double *c = new double [nr*a_nc];
-  double *ctmp = (double *) NULL;
-
-  int a_len = a.length ();
-  for (int j = 0; j < a_len; j++)
-    {
-      int idx = j * nr;
-      ctmp = c + idx;
-      if (a.elem (j, j) == 1.0)
-	{
-	  for (int i = 0; i < nr; i++)
-	    ctmp[i] = m.elem (i, j);
-	}
-      else if (a.elem (j, j) == 0.0)
-	{
-	  for (int i = 0; i < nr; i++)
-	    ctmp[i] = 0.0;
-	}
-      else
-	{
-	  for (int i = 0; i < nr; i++)
-	    ctmp[i] = a.elem (j, j) * m.elem (i, j);
-	}
-    }
-
-  if (a_nr < a_nc)
-    {
-      for (int i = nr * nc; i < nr * a_nc; i++)
-	ctmp[i] = 0.0;
-    }
-
-  return Matrix (c, nr, a_nc);
-}
-
-ComplexMatrix
-operator + (const Matrix& 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 Matrix& 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 Matrix& m, const ComplexDiagMatrix& 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);
-
-  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_nr < 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
-
-Matrix
-operator * (const Matrix& 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 Matrix ();
-    }
-
-  if (nr == 0 || nc == 0 || a_nc == 0)
-    return Matrix (nr, a_nc, 0.0);
-
-  char trans  = 'N';
-  char transa = 'N';
-
-  int ld  = nr;
-  int lda = a_nr;
-
-  double alpha = 1.0;
-  double beta  = 0.0;
-
-  double *c = new double [nr*a_nc];
-
-  F77_FCN (dgemm) (&trans, &transa, &nr, &a_nc, &nc, &alpha, m.data (),
-		   &ld, a.data (), &lda, &beta, c, &nr, 1L, 1L);
-
-  return Matrix (c, nr, a_nc);
-}
-
-ComplexMatrix
-operator * (const Matrix& m, const ComplexMatrix& a)
-{
-  ComplexMatrix tmp (m);
-  return tmp * a;
-}
-
-ComplexMatrix
-operator + (const Matrix& 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 ();
-    }
-
-  return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
-}
-
-ComplexMatrix
-operator - (const Matrix& m, const ComplexMatrix& a)
-{
-  int nr = m.rows ();
-  int nc = m.cols ();
-  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
-product (const Matrix& 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 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 Matrix& 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 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.
-
-Matrix
-map (d_d_Mapper f, const Matrix& a)
-{
-  Matrix b (a);
-  b.map (f);
-  return b;
-}
-
-void
-Matrix::map (d_d_Mapper f)
-{
-  double *d = fortran_vec (); // Ensures only one reference to my privates!
-
-  for (int i = 0; i < length (); i++)
-    d[i] = f (d[i]);
-}
-
-// XXX FIXME XXX Do these really belong here?  They should maybe be
-// cleaned up a bit, no?  What about corresponding functions for the
-// Vectors?
-
-Matrix
-Matrix::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
-Matrix::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;
-}
-
-Matrix
-Matrix::cumprod (void) const
-{
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr == 1)
-    {
-      retval.resize (1, nc);
-      if (nc > 0)
-	{
-	  double 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);
-      if (nr > 0)
-	{
-	  double 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);
-      if (nr > 0 && nc > 0)
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      double 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;
-}
-
-Matrix
-Matrix::cumsum (void) const
-{
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr == 1)
-    {
-      retval.resize (1, nc);
-      if (nc > 0)
-	{
-	  double 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);
-      if (nr > 0)
-	{
-	  double 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);
-      if (nr > 0 && nc > 0)
-	{
-	  for (int j = 0; j < nc; j++)
-	    {
-	      double 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;
-}
-
-Matrix
-Matrix::prod (void) const
-{
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  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
-    {
-      if (nc == 0)
-	{
-	  retval.resize (1, 1);
-	  retval.elem (0, 0) = 1.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;
-}
-
-Matrix
-Matrix::sum (void) const
-{
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  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
-    {
-      if (nc == 0)
-	{
-	  retval.resize (1, 1);
-	  retval.elem (0, 0) = 0.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;
-}
-
-Matrix
-Matrix::sumsq (void) const
-{
-  Matrix retval;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr == 1)
-    {
-      retval.resize (1, 1);
-      retval.elem (0, 0) = 0.0;
-      for (int j = 0; j < nc; j++)
-	{
-	  double 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++)
-	{
-	  double 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++)
-	    {
-	      double d = elem (i, j);
-	      retval.elem (0, j) += d * d;
-	    }
-	}
-    }
-  return retval;
-}
-
-ColumnVector
-Matrix::diag (void) const
-{
-  return diag (0);
-}
-
-ColumnVector
-Matrix::diag (int k) const
-{
-  int nnr = rows ();
-  int nnc = cols ();
-  if (k > 0)
-    nnc -= k;
-  else if (k < 0)
-    nnr += k;
-
-  ColumnVector 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;
-}
-
-ColumnVector
-Matrix::row_min (void) const
-{
-  ColumnVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nr);
-
-      for (int i = 0; i < nr; i++)
-	{
-	  double res = elem (i, 0);
-	  for (int j = 1; j < nc; j++)
-	    if (elem (i, j) < res)
-	      res = elem (i, j);
-	  result.elem (i) = res;
-	}
-    }
-
-  return result;
-}
-
-ColumnVector
-Matrix::row_min_loc (void) const
-{
-  ColumnVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nr);
-
-      for (int i = 0; i < nr; i++)
-        {
-          int res = 0;
-          for (int j = 0; j < nc; j++)
-            if (elem (i, j) < elem (i, res))
-              res = j;
-          result.elem (i) = (double) (res + 1);
-        }
-    }
-
-  return result;
-}
-
-ColumnVector
-Matrix::row_max (void) const
-{
-  ColumnVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nr);
-
-      for (int i = 0; i < nr; i++)
-	{
-	  double res = elem (i, 0);
-	  for (int j = 1; j < nc; j++)
-	    if (elem (i, j) > res)
-	      res = elem (i, j);
-	  result.elem (i) = res;
-	}
-    }
-
-  return result;
-}
-
-ColumnVector
-Matrix::row_max_loc (void) const
-{
-  ColumnVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nr);
-
-      for (int i = 0; i < nr; i++)
-        {
-          int res = 0;
-          for (int j = 0; j < nc; j++)
-            if (elem (i, j) > elem (i, res))
-              res = j;
-          result.elem (i) = (double) (res + 1);
-        }
-    }
-
-  return result;
-}
-
-RowVector
-Matrix::column_min (void) const
-{
-  RowVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nc);
-
-      for (int j = 0; j < nc; j++)
-	{
-	  double res = elem (0, j);
-	  for (int i = 1; i < nr; i++)
-	    if (elem (i, j) < res)
-	      res = elem (i, j);
-	  result.elem (j) = res;
-	}
-    }
-
-  return result;
-}
-RowVector
-Matrix::column_min_loc (void) const
-{
-  RowVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nc);
-
-      for (int j = 0; j < nc; j++)
-        {
-          int res = 0;
-          for (int i = 0; i < nr; i++)
-            if (elem (i, j) < elem (res, j))
-              res = i;
-          result.elem (j) = (double) (res + 1);
-        }
-    }
-
-  return result;
-}
-
-
-RowVector
-Matrix::column_max (void) const
-{
-  RowVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nc);
-
-      for (int j = 0; j < nc; j++)
-	{
-	  double res = elem (0, j);
-	  for (int i = 1; i < nr; i++)
-	    if (elem (i, j) > res)
-	      res = elem (i, j);
-	  result.elem (j) = res;
-	}
-    }
-
-  return result;
-}
-
-RowVector
-Matrix::column_max_loc (void) const
-{
-  RowVector result;
-
-  int nr = rows ();
-  int nc = cols ();
-
-  if (nr > 0 && nc > 0)
-    {
-      result.resize (nc);
-
-      for (int j = 0; j < nc; j++)
-        {
-          int res = 0;
-          for (int i = 0; i < nr; i++)
-            if (elem (i, j) > elem (res, j))
-              res = i;
-          result.elem (j) = (double) (res + 1);
-        }
-    }
-
-  return result;
-}
-
-ostream&
-operator << (ostream& os, const Matrix& 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, Matrix& a)
-{
-  int nr = a.rows ();
-  int nc = a.cols ();
-
-  if (nr < 1 || nc < 1)
-    is.clear (ios::badbit);
-  else
-    {
-      double 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;
-}
-
-/*
- * 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;
-      return ComplexDET (d);
-    }
-
-  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;
-    }
-  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;
-}
-
-/*
- * Read an array of data froma file in binary format.
- */
-int
-Matrix::read (FILE *fptr, int size, Matrix::conversion conv)
-{
-// Allocate buffer pointers.
-
-  union
-    {
-      void *vd;
-      char *ch;
-      u_char *uc;
-//    s_char *sc; // Some systems may need this?
-      short *sh;
-      u_short *us;
-      int *in;
-      u_int *ui;
-      long *ln;
-      u_long *ul;
-      float *fl;
-      double *db;
-    }
-  buf;
-
-  buf.db = fortran_vec ();
-
-// Read data directly into matrix data array.
-
-  int count = fread (buf.ch, size, length (), fptr);
-
-// Convert data to double.
-
-  int k;
-
-  switch (conv)
-    {
-    case CNV_DOUBLE:
-      break;
-
-    case CNV_CHAR:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.ch[k];
-      break;
-
-    case CNV_UCHAR:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.uc[k];
-      break;
-
-// Some systems may need this??
-//    case CNV_SCHAR:
-//      for (k = count - 1; k > -1; k--)
-//	buf.db[k] = buf.sc[k];
-//      break;
-
-    case CNV_SHORT:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.sh[k];
-      break;
-
-    case CNV_USHORT:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.us[k];
-      break;
-
-    case CNV_INT:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.in[k];
-      break;
-
-    case CNV_UINT:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.ui[k];
-      break;
-
-    case CNV_LONG:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.ln[k];
-      break;
-
-    case CNV_ULONG:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.ul[k];
-      break;
-
-    case CNV_FLOAT:
-      for (k = count - 1; k > -1; k--)
-	buf.db[k] = buf.fl[k];
-      break;
-
-    default:
-      return 0;
-    }
-
-  return count;
-}
-
-/*
- * Write the data array to a file in binary format.
- */
-int
-Matrix::write (FILE *fptr, int size, Matrix::conversion conv)
-{
-// Allocate buffer pointers.
-
-  union
-    {
-      void *vd;
-      char *ch;
-      u_char *uc;
-//    s_char *sc; // Some systems may need this?
-      short *sh;
-      u_short *us;
-      int *in;
-      u_int *ui;
-      long *ln;
-      u_long *ul;
-      float *fl;
-      double *db;
-    }
-  buf;
-
-  int len = length ();
-
-  if (conv != CNV_DOUBLE)
-    buf.db = new double [len];
-
-  double *bufi = fortran_vec ();
-
-// Convert from double to correct size.
-
-  int k;
-
-  switch (conv)
-    {
-    case CNV_DOUBLE:
-      buf.db = bufi;
-      break;
-
-    case CNV_CHAR:
-      for (k = 0; k < len; k++)
-	buf.ch[k] = (char) bufi[k];
-      break;
-
-    case CNV_UCHAR:
-      for (k = 0; k < len; k++)
-	buf.uc[k] = (u_char) bufi[k];
-      break;
-
-// Some systems may need this?
-//    case CNV_SCHAR:
-//      for (k = 0; k < len; k++)
-//	buf.uc[k] = (s_char) bufi[k];
-//      break;
-
-    case CNV_SHORT:
-      for (k = 0; k < len; k++)
-	buf.sh[k] = (short) bufi[k];
-      break;
-
-    case CNV_USHORT:
-      for (k = 0; k < len; k++)
-	buf.us[k] = (u_short) bufi[k];
-      break;
-
-    case CNV_INT:
-      for (k = 0; k < len; k++)
-	buf.in[k] = (int) bufi[k];
-      break;
-
-    case CNV_UINT:
-      for (k = 0; k < len; k++)
-	buf.ui[k] = (u_int) bufi[k];
-      break;
-
-    case CNV_LONG:
-      for (k = 0; k < len; k++)
-	buf.ln[k] = (long) bufi[k];
-      break;
-
-    case CNV_ULONG:
-      for (k = 0; k < len; k++)
-	buf.ul[k] = (u_long) bufi[k];
-      break;
-
-    case CNV_FLOAT:
-      for (k = 0; k < len; k++)
-	buf.fl[k] = (float) bufi[k];
-      break;
-
-    default:
-      return 0;
-  }
-
-// Write data from converted matrix data array.
-
-  int count = fwrite (buf.ch, size, length (), fptr);
-
-  if (conv != CNV_DOUBLE)
-    delete [] buf.db;
-
-  return count;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; page-delimiter: "^/\\*" ***
-;;; End: ***
-*/