# HG changeset patch
# User jwe
# Date 770864274 0
# Node ID 00f8b2242a18b6c5e75d729fbbda710408a05527
# Parent  c45e50c4ee94551b404e4dea709606dae2ad3deb
[project @ 1994-06-06 00:45:57 by jwe]

diff --git a/liboctave/Bounds.cc b/liboctave/Bounds.cc
--- a/liboctave/Bounds.cc
+++ b/liboctave/Bounds.cc
@@ -1,7 +1,7 @@
 // Bounds.cc                                              -*- C++ -*-
 /*
 
-Copyright (C) 1992, 1993 John W. Eaton
+Copyright (C) 1992, 1993, 1994 John W. Eaton
 
 This file is part of Octave.
 
@@ -25,6 +25,10 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include <iostream.h>
 
 #include "Bounds.h"
diff --git a/liboctave/Bounds.h b/liboctave/Bounds.h
--- a/liboctave/Bounds.h
+++ b/liboctave/Bounds.h
@@ -24,12 +24,16 @@
 #if !defined (octave_Bounds_h)
 #define octave_Bounds_h 1
 
-#include "Matrix.h"
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+class ostream;
+
+#include "dColVector.h"
 
 extern "C++" {
 
-class ostream;
-
 #ifndef Vector
 #define Vector ColumnVector
 #endif
diff --git a/liboctave/LSODE.cc b/liboctave/LSODE.cc
--- a/liboctave/LSODE.cc
+++ b/liboctave/LSODE.cc
@@ -25,8 +25,13 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include <math.h>
 #include <float.h>
+#include <iostream.h>
 
 #include "ODE.h"
 #include "f77-uscore.h"
@@ -73,8 +78,8 @@
       rwork[i] = 0.0;
     }
 
-  fun = NULL;
-  jac = NULL;
+  fun = 0;
+  jac = 0;
 }
 
 ODE::ODE (int size)
@@ -104,8 +109,8 @@
       rwork[i] = 0.0;
     }
 
-  fun = NULL;
-  jac = NULL;
+  fun = 0;
+  jac = 0;
 }
 
 ODE::ODE (const ColumnVector& state, double time, const ODEFunc& f)
@@ -194,10 +199,10 @@
 ColumnVector
 ODE::integrate (double tout)
 {
-  if (jac == NULL)
+  if (jac)
+    method_flag = 21;
+  else
     method_flag = 22;
-  else
-    method_flag = 21;
 
   integration_error = 0;
 
diff --git a/liboctave/Matrix.h b/liboctave/Matrix.h
--- a/liboctave/Matrix.h
+++ b/liboctave/Matrix.h
@@ -31,2448 +31,9 @@
 #if !defined (octave_Matrix_h)
 #define octave_Matrix_h 1
 
-#include <stdlib.h>
-#include <stddef.h>
-#include <math.h>
-#include <assert.h>
-#include <iostream.h>
-#include <stdio.h>
-// #include <iomanip.h>  // We don\'t use this yet.
-#include <Complex.h>
-
-extern "C++" {
-
-class ostream;
-class istream;
-
-#ifndef MAPPER_FCN_TYPEDEFS
-#define MAPPER_FCN_TYPEDEFS 1
-
-typedef double (*d_d_Mapper)(double);
-typedef double (*d_c_Mapper)(const Complex&);
-typedef Complex (*c_c_Mapper)(const Complex&);
-
-#endif
-
-#include "Array.h"
-
-// Classes we declare.
-
-class Matrix;
-class ColumnVector;
-class RowVector;
-class DiagMatrix;
-class ComplexMatrix;
-class ComplexColumnVector;
-class ComplexRowVector;
-class ComplexDiagMatrix;
-class AEPBALANCE;
-class ComplexAEPBALANCE;
-class GEPBALANCE;
-class CHOL;
-class ComplexCHOL;
-class DET;
-class ComplexDET;
-class EIG;
-class HESS;
-class ComplexHESS;
-class SCHUR;
-class ComplexSCHUR;
-class SVD;
-class ComplexSVD;
-class LU;
-class ComplexLU;
-class QR;
-class ComplexQR;
-
-/*
- * Matrix class
- */
-
-class Matrix : public Array2<double>
-{
-friend class LU;
-friend class SVD;
-
-public:
-
-  Matrix (void) : Array2<double> () { }
-  Matrix (int r, int c) : Array2<double> (r, c) { }
-  Matrix (int r, int c, double val) : Array2<double> (r, c, val) { }
-  Matrix (const Array2<double>& a) : Array2<double> (a) { }
-  Matrix (const Matrix& a) : Array2<double> (a) { }
-  Matrix (const DiagArray<double>& a) : Array2<double> (a) { }
-  Matrix (const DiagMatrix& a);
-//  Matrix (double a) : Array2<double> (1, 1, a) { }
-
-  Matrix& operator = (const Matrix& a)
-    {
-      Array2<double>::operator = (a);
-      return *this;
-    }
-
-  int operator == (const Matrix& a) const;
-  int operator != (const Matrix& a) const;
-
-// destructive insert/delete/reorder operations
-
-  Matrix& insert (const Matrix& a, int r, int c);
-  Matrix& insert (const RowVector& a, int r, int c);
-  Matrix& insert (const ColumnVector& a, int r, int c);
-  Matrix& insert (const DiagMatrix& a, int r, int c);
-
-  Matrix& fill (double val);
-  Matrix& fill (double val, int r1, int c1, int r2, int c2);
-
-  Matrix append (const Matrix& a) const;
-  Matrix append (const RowVector& a) const;
-  Matrix append (const ColumnVector& a) const;
-  Matrix append (const DiagMatrix& a) const;
-
-  Matrix stack (const Matrix& a) const;
-  Matrix stack (const RowVector& a) const;
-  Matrix stack (const ColumnVector& a) const;
-  Matrix stack (const DiagMatrix& a) const;
-
-  Matrix transpose (void) const;
-
-// resize is the destructive equivalent for this one
-
-  Matrix extract (int r1, int c1, int r2, int c2) const;
-
-// extract row or column i.
-
-  RowVector row (int i) const;
-  RowVector row (char *s) const;
-
-  ColumnVector column (int i) const;
-  ColumnVector column (char *s) const;
-
-  Matrix inverse (void) const;
-  Matrix inverse (int& info) const;
-  Matrix inverse (int& info, double& rcond) const;
-
-  ComplexMatrix fourier (void) const;
-  ComplexMatrix ifourier (void) const;
-
-  DET determinant (void) const;
-  DET determinant (int& info) const;
-  DET determinant (int& info, double& rcond) const;
-
-  Matrix solve (const Matrix& b) const;
-  Matrix solve (const Matrix& b, int& info) const;
-  Matrix solve (const Matrix& b, int& info, double& rcond) const;
-
-  ComplexMatrix solve (const ComplexMatrix& b) const;
-  ComplexMatrix solve (const ComplexMatrix& b, int& info) const;
-  ComplexMatrix solve (const ComplexMatrix& b, int& info, double& rcond) const;
-
-  ColumnVector solve (const ColumnVector& b) const;
-  ColumnVector solve (const ColumnVector& b, int& info) const;
-  ColumnVector solve (const ColumnVector& b, int& info, double& rcond) const;
-
-  ComplexColumnVector solve (const ComplexColumnVector& b) const;
-  ComplexColumnVector solve (const ComplexColumnVector& b, int& info) const;
-  ComplexColumnVector solve (const ComplexColumnVector& b, int& info,
-			     double& rcond) const;
-
-  Matrix lssolve (const Matrix& b) const;
-  Matrix lssolve (const Matrix& b, int& info) const;
-  Matrix lssolve (const Matrix& b, int& info, int& rank) const;
-
-  ComplexMatrix lssolve (const ComplexMatrix& b) const;
-  ComplexMatrix lssolve (const ComplexMatrix& b, int& info) const;
-  ComplexMatrix lssolve (const ComplexMatrix& b, int& info,
-			 int& rank) const;
-
-  ColumnVector lssolve (const ColumnVector& b) const;
-  ColumnVector lssolve (const ColumnVector& b, int& info) const;
-  ColumnVector lssolve (const ColumnVector& b, int& info, int& rank) const;
-
-  ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info,
-			       int& rank) const;
-
-  Matrix& operator += (const Matrix& a);
-  Matrix& operator -= (const Matrix& a);
-
-  Matrix& operator += (const DiagMatrix& a);
-  Matrix& operator -= (const DiagMatrix& a);
-
-// unary operations
-
-  Matrix operator ! (void) const;
-
-// matrix by scalar -> matrix operations
-
-  friend ComplexMatrix operator + (const Matrix& a, const Complex& s);
-  friend ComplexMatrix operator - (const Matrix& a, const Complex& s);
-  friend ComplexMatrix operator * (const Matrix& a, const Complex& s);
-  friend ComplexMatrix operator / (const Matrix& a, const Complex& s);
-
-// scalar by matrix -> matrix operations
-
-  friend ComplexMatrix operator + (const Complex& s, const Matrix& a);
-  friend ComplexMatrix operator - (const Complex& s, const Matrix& a);
-  friend ComplexMatrix operator * (const Complex& s, const Matrix& a);
-  friend ComplexMatrix operator / (const Complex& s, const Matrix& a);
-
-// matrix by column vector -> column vector operations
-
-  friend ColumnVector operator * (const Matrix& a, const ColumnVector& b);
-  friend ComplexColumnVector operator * (const Matrix& a,
-					 const ComplexColumnVector& b);
-
-// matrix by diagonal matrix -> matrix operations
-
-  friend Matrix operator + (const Matrix& a, const DiagMatrix& b);
-  friend Matrix operator - (const Matrix& a, const DiagMatrix& b);
-  friend Matrix operator * (const Matrix& a, const DiagMatrix& b);
-
-  friend ComplexMatrix operator + (const Matrix& a,
-				   const ComplexDiagMatrix& b); 
-  friend ComplexMatrix operator - (const Matrix& a,
-				   const ComplexDiagMatrix& b);
-  friend ComplexMatrix operator * (const Matrix& a,
-				   const ComplexDiagMatrix& b);
-
-// matrix by matrix -> matrix operations
-
-  friend Matrix operator * (const Matrix& a, const Matrix& b);
-  friend ComplexMatrix operator * (const Matrix& a, const ComplexMatrix& b);
-
-  friend ComplexMatrix operator + (const Matrix& a, const ComplexMatrix& b);
-  friend ComplexMatrix operator - (const Matrix& a, const ComplexMatrix& b);
-
-  friend ComplexMatrix product (const Matrix& a, const ComplexMatrix& b);
-  friend ComplexMatrix quotient (const Matrix& a, const ComplexMatrix& b);
-
-// other operations
-
-  friend Matrix map (d_d_Mapper f, const Matrix& a);
-  void map (d_d_Mapper f);
-
-  Matrix all (void) const;
-  Matrix any (void) const;
-
-  Matrix cumprod (void) const;
-  Matrix cumsum (void) const;
-  Matrix prod (void) const;
-  Matrix sum (void) const;
-  Matrix sumsq (void) const;
-
-  ColumnVector diag (void) const;
-  ColumnVector diag (int k) const;
-
-  ColumnVector row_min (void) const;
-  ColumnVector row_min_loc (void) const;
-
-  ColumnVector row_max (void) const;
-  ColumnVector row_max_loc (void) const;
-
-  RowVector column_min (void) const;
-  RowVector column_min_loc (void) const;
-
-  RowVector column_max (void) const;
-  RowVector column_max_loc (void) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const Matrix& a);
-  friend istream& operator >> (istream& is, Matrix& a);
-
-  enum conversion
-    {
-      CNV_UNKNOWN,
-      CNV_UCHAR,
-      CNV_CHAR,
-      CNV_SCHAR,
-      CNV_SHORT,
-      CNV_USHORT,
-      CNV_INT,
-      CNV_UINT,
-      CNV_LONG,
-      CNV_ULONG,
-      CNV_FLOAT,
-      CNV_DOUBLE,
-    };
-
-
-  int read (FILE *fptr, int size, Matrix::conversion);
-  int write (FILE *fptr, int size, Matrix::conversion);
-
-// Until templates really work with g++:
-
-#define KLUDGE_MATRICES
-#define TYPE double
-#define KL_MAT_TYPE Matrix
-#include "mx-kludge.h"
-#undef KLUDGE_MATRICES
-#undef TYPE
-#undef KL_MAT_TYPE
-
-private:
-
-  Matrix (double *d, int r, int c) : Array2<double> (d, r, c) { }
-};
-
-/*
- * Column Vector class
- */
-
-class ColumnVector : public Array<double>
-{
-friend class RowVector;
-
-public:
-
-  ColumnVector (void) : Array<double> () { }
-  ColumnVector (int n) : Array<double> (n) { }
-  ColumnVector (int n, double val) : Array<double> (n, val) { }
-  ColumnVector (const Array<double>& a) : Array<double> (a) { }
-  ColumnVector (const ColumnVector& a) : Array<double> (a) { }
-//  ColumnVector (double a) : Array<double> (1, a) { }
-
-  ColumnVector& operator = (const ColumnVector& a)
-    {
-      Array<double>::operator = (a);
-      return *this;
-    }
-
-//  operator Array<double>& () const { return *this; }
-
-  int operator == (const ColumnVector& a) const;
-  int operator != (const ColumnVector& a) const;
-
-// destructive insert/delete/reorder operations
-
-  ColumnVector& insert (const ColumnVector& a, int r);
-
-  ColumnVector& fill (double val);
-  ColumnVector& fill (double val, int r1, int r2);
-
-  ColumnVector stack (const ColumnVector& a) const;
-
-  RowVector transpose (void) const;
-
-// resize is the destructive equivalent for this one
-
-  ColumnVector extract (int r1, int r2) const;
-
-// column vector by column vector -> column vector operations
-
-  ColumnVector& operator += (const ColumnVector& a);
-  ColumnVector& operator -= (const ColumnVector& a);
-
-// column vector by scalar -> column vector operations
-
-  friend ComplexColumnVector operator + (const ColumnVector& a,
-					 const Complex& s);  
-  friend ComplexColumnVector operator - (const ColumnVector& a,
-					 const Complex& s);
-  friend ComplexColumnVector operator * (const ColumnVector& a,
-					 const Complex& s);
-  friend ComplexColumnVector operator / (const ColumnVector& a,
-					 const Complex& s);
-
-// scalar by column vector -> column vector operations
-
-  friend ComplexColumnVector operator + (const Complex& s,
-					 const ColumnVector& a); 
-  friend ComplexColumnVector operator - (const Complex& s,
-					 const ColumnVector& a);
-  friend ComplexColumnVector operator * (const Complex& s,
-					 const ColumnVector& a);
-  friend ComplexColumnVector operator / (const Complex& s,
-					 const ColumnVector& a);
-
-// column vector by row vector -> matrix operations
-
-  friend Matrix operator * (const ColumnVector& a, const RowVector& a);
-
-  friend ComplexMatrix operator * (const ColumnVector& a,
-				   const ComplexRowVector& b);
-
-// column vector by column vector -> column vector operations
-
-  friend ComplexColumnVector operator + (const ComplexColumnVector& a,
-					 const ComplexColumnVector& b);
-
-  friend ComplexColumnVector operator - (const ComplexColumnVector& a,
-					 const ComplexColumnVector& b); 
-
-  friend ComplexColumnVector product (const ComplexColumnVector& a,
-				      const ComplexColumnVector& b); 
-
-  friend ComplexColumnVector quotient (const ComplexColumnVector& a,
-				       const ComplexColumnVector& b); 
-
-// other operations
-
-  friend ColumnVector map (d_d_Mapper f, const ColumnVector& a);
-  void map (d_d_Mapper f);
-
-  double min (void) const;
-  double max (void) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const ColumnVector& a);
-  friend ostream& operator >> (ostream& is, ColumnVector& a);
-
-#define KLUDGE_VECTORS
-#define TYPE double
-#define KL_VEC_TYPE ColumnVector
-#include "mx-kludge.h"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-private:
-
-  ColumnVector (double *d, int l) : Array<double> (d, l) { }
-};
-
-/*
- * Row Vector class
- */
-
-class RowVector : public Array<double>
-{
-friend class ColumnVector;
-
-public:
-
-  RowVector (void) : Array<double> () { }
-  RowVector (int n) : Array<double> (n) { }
-  RowVector (int n, double val) : Array<double> (n, val) { }
-  RowVector (const Array<double>& a) : Array<double> (a) { }
-  RowVector (const RowVector& a) : Array<double> (a) { }
-//  RowVector (double a) : Array<double> (1, a) { }
-
-  RowVector& operator = (const RowVector& a)
-    {
-      Array<double>::operator = (a);
-      return *this;
-    }
-
-//  operator Array<double>& () const { return *this; }
-
-  int operator == (const RowVector& a) const;
-  int operator != (const RowVector& a) const;
-
-// destructive insert/delete/reorder operations
-
-  RowVector& insert (const RowVector& a, int c);
-
-  RowVector& fill (double val);
-  RowVector& fill (double val, int c1, int c2);
-
-  RowVector append (const RowVector& a) const;
-
-  ColumnVector transpose (void) const;
-
-// resize is the destructive equivalent for this one
-
-  RowVector extract (int c1, int c2) const;
-
-// row vector by row vector -> row vector operations
-
-  RowVector& operator += (const RowVector& a);
-  RowVector& operator -= (const RowVector& a);
-
-// row vector by scalar -> row vector operations
-
-  friend ComplexRowVector operator + (const RowVector& a, const Complex& s);
-  friend ComplexRowVector operator - (const RowVector& a, const Complex& s);
-  friend ComplexRowVector operator * (const RowVector& a, const Complex& s);
-  friend ComplexRowVector operator / (const RowVector& a, const Complex& s);
-
-// scalar by row vector -> row vector operations
-
-  friend ComplexRowVector operator + (const Complex& s, const RowVector& a);
-  friend ComplexRowVector operator - (const Complex& s, const RowVector& a);
-  friend ComplexRowVector operator * (const Complex& s, const RowVector& a);
-  friend ComplexRowVector operator / (const Complex& s, const RowVector& a);
-
-// row vector by column vector -> scalar
-
-  friend double operator * (const RowVector& a, const ColumnVector& b);
-
-  friend Complex operator * (const RowVector& a, const ComplexColumnVector& b);
-
-// row vector by matrix -> row vector
-
-  friend RowVector operator * (const RowVector& a, const Matrix& b);
-
-  friend ComplexRowVector operator * (const RowVector& a,
-				      const ComplexMatrix& b);
-
-// row vector by row vector -> row vector operations
-
-  friend ComplexRowVector operator + (const RowVector& a,
-				      const ComplexRowVector& b);
-  friend ComplexRowVector operator - (const RowVector& a,
-				      const ComplexRowVector& b);
-
-  friend ComplexRowVector product (const RowVector& a,
-				   const ComplexRowVector& b);
-  friend ComplexRowVector quotient (const RowVector& a,
-				    const ComplexRowVector& b);
-
-// other operations
-
-  friend RowVector map (d_d_Mapper f, const RowVector& a);
-  void map (d_d_Mapper f);
-
-  double min (void) const;
-  double max (void) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const RowVector& a);
-  friend ostream& operator >> (ostream& is, RowVector& a);
-
-#define KLUDGE_VECTORS
-#define TYPE double
-#define KL_VEC_TYPE RowVector
-#include "mx-kludge.h"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-private:
-
-  RowVector (double *d, int l) : Array<double> (d, l) { }
-};
-
-/*
- * Diagonal Matrix class
- */
-
-class DiagMatrix : public DiagArray<double>
-{
-friend class SVD;
-friend class ComplexSVD;
-
-public:
-
-  DiagMatrix (void) : DiagArray<double> () { }
-  DiagMatrix (int n) : DiagArray<double> (n) { }
-  DiagMatrix (int n, double val) : DiagArray<double> (n, val) { }
-  DiagMatrix (int r, int c) : DiagArray<double> (r, c) { }
-  DiagMatrix (int r, int c, double val) : DiagArray<double> (r, c, val) { }
-  DiagMatrix (const RowVector& a) : DiagArray<double> (a) { }
-  DiagMatrix (const ColumnVector& a) : DiagArray<double> (a) { }
-  DiagMatrix (const DiagArray<double>& a) : DiagArray<double> (a) { }
-  DiagMatrix (const DiagMatrix& a) : DiagArray<double> (a) { }
-//  DiagMatrix (double a) : DiagArray<double> (1, a) { }
-
-  DiagMatrix& operator = (const DiagMatrix& a)
-    {
-      DiagArray<double>::operator = (a);
-      return *this;
-    }
-
-//  operator DiagArray<double>& () const { return *this; }
-
-  int operator == (const DiagMatrix& a) const;
-  int operator != (const DiagMatrix& a) const;
-
-  DiagMatrix& fill (double val);
-  DiagMatrix& fill (double val, int beg, int end);
-  DiagMatrix& fill (const ColumnVector& a);
-  DiagMatrix& fill (const RowVector& a);
-  DiagMatrix& fill (const ColumnVector& a, int beg);
-  DiagMatrix& fill (const RowVector& a, int beg);
-
-  DiagMatrix transpose (void) const;
-
-// resize is the destructive analog for this one
-
-  Matrix extract (int r1, int c1, int r2, int c2) const;
-
-// extract row or column i.
-
-  RowVector row (int i) const;
-  RowVector row (char *s) const;
-
-  ColumnVector column (int i) const;
-  ColumnVector column (char *s) const;
-
-  DiagMatrix inverse (void) const;
-  DiagMatrix inverse (int& info) const;
-
-// diagonal matrix by diagonal matrix -> diagonal matrix operations
-
-  DiagMatrix& operator += (const DiagMatrix& a);
-  DiagMatrix& operator -= (const DiagMatrix& a);
-
-// diagonal matrix by scalar -> matrix operations
-
-  friend Matrix operator + (const DiagMatrix& a, double s);
-  friend Matrix operator - (const DiagMatrix& a, double s);
-
-  friend ComplexMatrix operator + (const DiagMatrix& a, const Complex& s);
-  friend ComplexMatrix operator - (const DiagMatrix& a, const Complex& s);
-
-// diagonal matrix by scalar -> diagonal matrix operations
-
-  friend ComplexDiagMatrix operator * (const DiagMatrix& a, const Complex& s);
-  friend ComplexDiagMatrix operator / (const DiagMatrix& a, const Complex& s);
-
-// scalar by diagonal matrix -> matrix operations
-
-  friend Matrix operator + (double s, const DiagMatrix& a);
-  friend Matrix operator - (double s, const DiagMatrix& a);
-
-  friend ComplexMatrix operator + (const Complex& s, const DiagMatrix& a);
-  friend ComplexMatrix operator - (const Complex& s, const DiagMatrix& a);
-
-// scalar by diagonal matrix -> diagonal matrix operations
-
-  friend ComplexDiagMatrix operator * (const Complex& s, const DiagMatrix& a);
-
-// diagonal matrix by column vector -> column vector operations
-
-  friend ColumnVector operator * (const DiagMatrix& a, const ColumnVector& b);
-
-  friend ComplexColumnVector operator * (const DiagMatrix& a,
-					 const ComplexColumnVector& b);
-
-// diagonal matrix by diagonal matrix -> diagonal matrix operations
-
-  friend DiagMatrix operator * (const DiagMatrix& a,
-				const DiagMatrix& b);
-
-  friend ComplexDiagMatrix operator + (const DiagMatrix& a,
-				       const ComplexDiagMatrix& b);
-  friend ComplexDiagMatrix operator - (const DiagMatrix& a,
-				       const ComplexDiagMatrix& b);
-  friend ComplexDiagMatrix operator * (const DiagMatrix& a,
-				       const ComplexDiagMatrix& b);
-
-  friend ComplexDiagMatrix product (const DiagMatrix& a,
-				    const ComplexDiagMatrix& b);
-
-// diagonal matrix by matrix -> matrix operations
-
-  friend Matrix operator + (const DiagMatrix& a, const Matrix& b);
-  friend Matrix operator - (const DiagMatrix& a, const Matrix& b);
-  friend Matrix operator * (const DiagMatrix& a, const Matrix& b);
-
-  friend ComplexMatrix operator + (const DiagMatrix& a,
-				   const ComplexMatrix& b);
-  friend ComplexMatrix operator - (const DiagMatrix& a,
-				   const ComplexMatrix& b);
-  friend ComplexMatrix operator * (const DiagMatrix& a,
-				   const ComplexMatrix& b);
-
-// other operations
-
-  ColumnVector diag (void) const;
-  ColumnVector diag (int k) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const DiagMatrix& a);
-
-#define KLUDGE_DIAG_MATRICES
-#define TYPE double
-#define KL_DMAT_TYPE DiagMatrix
-#include "mx-kludge.h"
-#undef KLUDGE_DIAG_MATRICES
-#undef TYPE
-#undef KL_DMAT_TYPE
-
-private:
-
-  DiagMatrix (double *d, int nr, int nc) : DiagArray<double> (d, nr, nc) { }
-};
-
-/*
- * Complex Matrix class
- */
-
-class ComplexMatrix : public Array2<Complex>
-{
-friend class ComplexLU;
-friend class ComplexSVD;
-
-public:
- 
-  ComplexMatrix (void) : Array2<Complex> () { }
-  ComplexMatrix (int r, int c) : Array2<Complex> (r, c) { }
-  ComplexMatrix (int r, int c, const Complex& val)
-    : Array2<Complex> (r, c, val) { }
-  ComplexMatrix (const Matrix& a);
-  ComplexMatrix (const Array2<Complex>& a) : Array2<Complex> (a) { }
-  ComplexMatrix (const ComplexMatrix& a) : Array2<Complex> (a) { }
-  ComplexMatrix (const DiagMatrix& a);
-  ComplexMatrix (const DiagArray<Complex>& a) : Array2<Complex> (a) { }
-  ComplexMatrix (const ComplexDiagMatrix& a);
-//  ComplexMatrix (double a) : Array2<Complex> (1, 1, a) { }
-//  ComplexMatrix (const Complex& a) : Array2<Complex> (1, 1, a) { }
-
-  ComplexMatrix& operator = (const ComplexMatrix& a)
-    {
-      Array2<Complex>::operator = (a);
-      return *this;
-    }
-
-//  operator Array2<Complex>& () const { return *this; }
-
-  int operator == (const ComplexMatrix& a) const;
-  int operator != (const ComplexMatrix& a) const;
-
-// destructive insert/delete/reorder operations
-
-  ComplexMatrix& insert (const Matrix& a, int r, int c);
-  ComplexMatrix& insert (const RowVector& a, int r, int c);
-  ComplexMatrix& insert (const ColumnVector& a, int r, int c);
-  ComplexMatrix& insert (const DiagMatrix& a, int r, int c);
-
-  ComplexMatrix& insert (const ComplexMatrix& a, int r, int c);
-  ComplexMatrix& insert (const ComplexRowVector& a, int r, int c);
-  ComplexMatrix& insert (const ComplexColumnVector& a, int r, int c);
-  ComplexMatrix& insert (const ComplexDiagMatrix& a, int r, int c);
-
-  ComplexMatrix& fill (double val);
-  ComplexMatrix& fill (const Complex& val);
-  ComplexMatrix& fill (double val, int r1, int c1, int r2, int c2);
-  ComplexMatrix& fill (const Complex& val, int r1, int c1, int r2, int c2);
-
-  ComplexMatrix append (const Matrix& a) const;
-  ComplexMatrix append (const RowVector& a) const;
-  ComplexMatrix append (const ColumnVector& a) const;
-  ComplexMatrix append (const DiagMatrix& a) const;
-
-  ComplexMatrix append (const ComplexMatrix& a) const;
-  ComplexMatrix append (const ComplexRowVector& a) const;
-  ComplexMatrix append (const ComplexColumnVector& a) const;
-  ComplexMatrix append (const ComplexDiagMatrix& a) const;
-
-  ComplexMatrix stack (const Matrix& a) const;
-  ComplexMatrix stack (const RowVector& a) const;
-  ComplexMatrix stack (const ColumnVector& a) const;
-  ComplexMatrix stack (const DiagMatrix& a) const;
-
-  ComplexMatrix stack (const ComplexMatrix& a) const;
-  ComplexMatrix stack (const ComplexRowVector& a) const;
-  ComplexMatrix stack (const ComplexColumnVector& a) const;
-  ComplexMatrix stack (const ComplexDiagMatrix& a) const;
-
-  ComplexMatrix hermitian (void) const;  // complex conjugate transpose
-  ComplexMatrix transpose (void) const;
-
-  friend Matrix real (const ComplexMatrix& a);
-  friend Matrix imag (const ComplexMatrix& a);
-  friend ComplexMatrix conj (const ComplexMatrix& a);
-
-// resize is the destructive equivalent for this one
-
-  ComplexMatrix extract (int r1, int c1, int r2, int c2) const;
-
-// extract row or column i.
-
-  ComplexRowVector row (int i) const;
-  ComplexRowVector row (char *s) const;
-
-  ComplexColumnVector column (int i) const;
-  ComplexColumnVector column (char *s) const;
-
-  ComplexMatrix inverse (void) const;
-  ComplexMatrix inverse (int& info) const;
-  ComplexMatrix inverse (int& info, double& rcond) const;
-
-  ComplexMatrix fourier (void) const;
-  ComplexMatrix ifourier (void) const;
-
-  ComplexDET determinant (void) const;
-  ComplexDET determinant (int& info) const;
-  ComplexDET determinant (int& info, double& rcond) const;
-
-  ComplexMatrix solve (const Matrix& b) const;
-  ComplexMatrix solve (const Matrix& b, int& info) const;
-  ComplexMatrix solve (const Matrix& b, int& info, double& rcond) const;
-
-  ComplexMatrix solve (const ComplexMatrix& b) const;
-  ComplexMatrix solve (const ComplexMatrix& b, int& info) const;
-  ComplexMatrix solve (const ComplexMatrix& b, int& info, double& rcond) const;
-
-  ComplexColumnVector solve (const ComplexColumnVector& b) const;
-  ComplexColumnVector solve (const ComplexColumnVector& b, int& info) const;
-  ComplexColumnVector solve (const ComplexColumnVector& b, int& info,
-			     double& rcond) const;
-
-  ComplexMatrix lssolve (const ComplexMatrix& b) const;
-  ComplexMatrix lssolve (const ComplexMatrix& b, int& info) const;
-  ComplexMatrix lssolve (const ComplexMatrix& b, int& info,
-			 int& rank) const;
-
-  ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info,
-			       int& rank) const;
-
-// matrix by diagonal matrix -> matrix operations
-
-  ComplexMatrix& operator += (const DiagMatrix& a);
-  ComplexMatrix& operator -= (const DiagMatrix& a);
-
-  ComplexMatrix& operator += (const ComplexDiagMatrix& a);
-  ComplexMatrix& operator -= (const ComplexDiagMatrix& a);
-
-// matrix by matrix -> matrix operations
-
-  ComplexMatrix& operator += (const Matrix& a);
-  ComplexMatrix& operator -= (const Matrix& a);
-
-  ComplexMatrix& operator += (const ComplexMatrix& a);
-  ComplexMatrix& operator -= (const ComplexMatrix& a);
-
-// unary operations
-
-  Matrix operator ! (void) const;
-
-// matrix by scalar -> matrix operations
-
-  friend ComplexMatrix operator + (const ComplexMatrix& a, double s);
-  friend ComplexMatrix operator - (const ComplexMatrix& a, double s);
-  friend ComplexMatrix operator * (const ComplexMatrix& a, double s);
-  friend ComplexMatrix operator / (const ComplexMatrix& a, double s);
-
-// scalar by matrix -> matrix operations
-
-  friend ComplexMatrix operator + (double s, const ComplexMatrix& a);
-  friend ComplexMatrix operator - (double s, const ComplexMatrix& a);
-  friend ComplexMatrix operator * (double s, const ComplexMatrix& a);
-  friend ComplexMatrix operator / (double s, const ComplexMatrix& a);
-
-// matrix by column vector -> column vector operations
-
-  friend ComplexColumnVector operator * (const ComplexMatrix& a,
-					 const ColumnVector& b);
-
-  friend ComplexColumnVector operator * (const ComplexMatrix& a,
-					 const ComplexColumnVector& b);
-
-// matrix by diagonal matrix -> matrix operations
-
-  friend ComplexMatrix operator + (const ComplexMatrix& a,
-				   const DiagMatrix& b);
-  friend ComplexMatrix operator - (const ComplexMatrix& a,
-				   const DiagMatrix& b);
-  friend ComplexMatrix operator * (const ComplexMatrix& a,
-				   const DiagMatrix& b);
-
-  friend ComplexMatrix operator + (const ComplexMatrix& a,
-				   const ComplexDiagMatrix& b);
-  friend ComplexMatrix operator - (const ComplexMatrix& a,
-				   const ComplexDiagMatrix& b);
-  friend ComplexMatrix operator * (const ComplexMatrix& a,
-				   const ComplexDiagMatrix& b);
-
-// matrix by matrix -> matrix operations
-
-  friend ComplexMatrix operator + (const ComplexMatrix& a, const Matrix& b);
-  friend ComplexMatrix operator - (const ComplexMatrix& a, const Matrix& b);
-
-  friend ComplexMatrix operator * (const ComplexMatrix& a, const Matrix& b);
-  friend ComplexMatrix operator * (const ComplexMatrix& a,
-				   const ComplexMatrix& b);
-
-  friend ComplexMatrix product (const ComplexMatrix& a, const Matrix& b);
-  friend ComplexMatrix quotient (const ComplexMatrix& a, const Matrix& b);
-
-// other operations
-
-  friend ComplexMatrix map (c_c_Mapper f, const ComplexMatrix& a);
-  friend Matrix map (d_c_Mapper f, const ComplexMatrix& a);
-  void map (c_c_Mapper f);
-
-  Matrix all (void) const;
-  Matrix any (void) const;
-
-  ComplexMatrix cumprod (void) const;
-  ComplexMatrix cumsum (void) const;
-  ComplexMatrix prod (void) const;
-  ComplexMatrix sum (void) const;
-  ComplexMatrix sumsq (void) const;
-
-  ComplexColumnVector diag (void) const;
-  ComplexColumnVector diag (int k) const;
-
-  ComplexColumnVector row_min (void) const;
-  ComplexColumnVector row_min_loc (void) const;
-
-  ComplexColumnVector row_max (void) const;
-  ComplexColumnVector row_max_loc (void) const;
-
-  ComplexRowVector column_min (void) const;
-  ComplexRowVector column_min_loc (void) const;
-
-  ComplexRowVector column_max (void) const;
-  ComplexRowVector column_max_loc (void) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const ComplexMatrix& a);
-  friend istream& operator >> (istream& is, ComplexMatrix& a);
-
-#define KLUDGE_MATRICES
-#define TYPE Complex
-#define KL_MAT_TYPE ComplexMatrix
-#include "mx-kludge.h"
-#undef KLUDGE_MATRICES
-#undef TYPE
-#undef KL_MAT_TYPE
-
-private:
-
-  ComplexMatrix (Complex *d, int r, int c) : Array2<Complex> (d, r, c) { }
-};
-
-/*
- * Complex Column Vector class
- */
-
-class ComplexColumnVector : public Array<Complex>
-{
-friend class ComplexRowVector;
-
-public:
-
-  ComplexColumnVector (void) : Array<Complex> () { }
-  ComplexColumnVector (int n) : Array<Complex> (n) { }
-  ComplexColumnVector (int n, const Complex& val)
-    : Array<Complex> (n, val) { }
-  ComplexColumnVector (const ColumnVector& a);
-  ComplexColumnVector (const Array<Complex>& a) : Array<Complex> (a) { }
-  ComplexColumnVector (const ComplexColumnVector& a) : Array<Complex> (a) { }
-//  ComplexColumnVector (double a) : Array<Complex> (1, a) { }
-//  ComplexColumnVector (const Complex& a) : Array<Complex> (1, a) { }
-
-  ComplexColumnVector& operator = (const ComplexColumnVector& a)
-    {
-      Array<Complex>::operator = (a);
-      return *this;
-    }
-
-//  operator Array<Complex>& () const { return *this; }
-
-  int operator == (const ComplexColumnVector& a) const;
-  int operator != (const ComplexColumnVector& a) const;
-
-// destructive insert/delete/reorder operations
-
-  ComplexColumnVector& insert (const ColumnVector& a, int r);
-  ComplexColumnVector& insert (const ComplexColumnVector& a, int r);
-
-  ComplexColumnVector& fill (double val);
-  ComplexColumnVector& fill (const Complex& val);
-  ComplexColumnVector& fill (double val, int r1, int r2);
-  ComplexColumnVector& fill (const Complex& val, int r1, int r2);
-
-  ComplexColumnVector stack (const ColumnVector& a) const;
-  ComplexColumnVector stack (const ComplexColumnVector& a) const;
-
-  ComplexRowVector hermitian (void) const;  // complex conjugate transpose.
-  ComplexRowVector transpose (void) const;
-
-  friend ColumnVector real (const ComplexColumnVector& a);
-  friend ColumnVector imag (const ComplexColumnVector& a);
-  friend ComplexColumnVector conj (const ComplexColumnVector& a);
-
-// resize is the destructive equivalent for this one
-
-  ComplexColumnVector extract (int r1, int r2) const;
-
-// column vector by column vector -> column vector operations
-
-  ComplexColumnVector& operator += (const ColumnVector& a);
-  ComplexColumnVector& operator -= (const ColumnVector& a);
-
-  ComplexColumnVector& operator += (const ComplexColumnVector& a);
-  ComplexColumnVector& operator -= (const ComplexColumnVector& a);
-
-// column vector by scalar -> column vector operations
-
-  friend ComplexColumnVector operator + (const ComplexColumnVector& a,
-					 double s);
-  friend ComplexColumnVector operator - (const ComplexColumnVector& a,
-					 double s);
-  friend ComplexColumnVector operator * (const ComplexColumnVector& a,
-					 double s);
-  friend ComplexColumnVector operator / (const ComplexColumnVector& a,
-					 double s);
-
-// scalar by column vector -> column vector operations
-
-  friend ComplexColumnVector operator + (double s,
-					 const ComplexColumnVector& a); 
-  friend ComplexColumnVector operator - (double s,
-					 const ComplexColumnVector& a);
-  friend ComplexColumnVector operator * (double s,
-					 const ComplexColumnVector& a);
-  friend ComplexColumnVector operator / (double s,
-					 const ComplexColumnVector& a);
-
-// column vector by row vector -> matrix operations
-
-  friend ComplexMatrix operator * (const ComplexColumnVector& a,
-				   const ComplexRowVector& b);
-
-// column vector by column vector -> column vector operations
-
-  friend ComplexColumnVector operator + (const ComplexColumnVector& a,
-					 const ColumnVector& b);
-  friend ComplexColumnVector operator - (const ComplexColumnVector& a,
-					 const ColumnVector& b);
-
-  friend ComplexColumnVector product (const ComplexColumnVector& a,
-				      const ColumnVector& b);
-  friend ComplexColumnVector quotient (const ComplexColumnVector& a,
-				       const ColumnVector& b);
-
-// other operations
-
-  friend ComplexColumnVector map (c_c_Mapper f, const ComplexColumnVector& a);
-  friend ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a);
-  void map (c_c_Mapper f);
-
-  Complex min (void) const;
-  Complex max (void) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const ComplexColumnVector& a);
-  friend ostream& operator >> (ostream& is, ComplexColumnVector& a);
-
-#define KLUDGE_VECTORS
-#define TYPE Complex
-#define KL_VEC_TYPE ComplexColumnVector
-#include "mx-kludge.h"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-private:
-
-  ComplexColumnVector (Complex *d, int l) : Array<Complex> (d, l) { }
-};
-
-/*
- * Complex Row Vector class
- */
-
-class ComplexRowVector : public Array<Complex>
-{
-friend class ComplexColumnVector;
-
-public:
-
-  ComplexRowVector (void) : Array<Complex> () { }
-  ComplexRowVector (int n) : Array<Complex> (n) { }
-  ComplexRowVector (int n, const Complex& val) : Array<Complex> (n, val) { }
-  ComplexRowVector (const RowVector& a);
-  ComplexRowVector (const Array<Complex>& a) : Array<Complex> (a) { }
-  ComplexRowVector (const ComplexRowVector& a) : Array<Complex> (a) { }
-//  ComplexRowVector (double a) : Array<Complex> (1, a) { }
-//  ComplexRowVector (const Complex& a) : Array<Complex> (1, a) { }
-
-  ComplexRowVector& operator = (const ComplexRowVector& a)
-    {
-      Array<Complex>::operator = (a);
-      return *this;
-    }
-
-//  operator Array<Complex>& () const { return *this; }
-
-  int operator == (const ComplexRowVector& a) const;
-  int operator != (const ComplexRowVector& a) const;
-
-// destructive insert/delete/reorder operations
-
-  ComplexRowVector& insert (const RowVector& a, int c);
-  ComplexRowVector& insert (const ComplexRowVector& a, int c);
-
-  ComplexRowVector& fill (double val);
-  ComplexRowVector& fill (const Complex& val);
-  ComplexRowVector& fill (double val, int c1, int c2);
-  ComplexRowVector& fill (const Complex& val, int c1, int c2);
-
-  ComplexRowVector append (const RowVector& a) const;
-  ComplexRowVector append (const ComplexRowVector& a) const;
-
-  ComplexColumnVector hermitian (void) const;  // complex conjugate transpose.
-  ComplexColumnVector transpose (void) const;
-
-  friend RowVector real (const ComplexRowVector& a);
-  friend RowVector imag (const ComplexRowVector& a);
-  friend ComplexRowVector conj (const ComplexRowVector& a);
-
-// resize is the destructive equivalent for this one
-
-  ComplexRowVector extract (int c1, int c2) const;
-
-// row vector by row vector -> row vector operations
-
-  ComplexRowVector& operator += (const RowVector& a);
-  ComplexRowVector& operator -= (const RowVector& a);
-
-  ComplexRowVector& operator += (const ComplexRowVector& a);
-  ComplexRowVector& operator -= (const ComplexRowVector& a);
-
-// row vector by scalar -> row vector operations
-
-  friend ComplexRowVector operator + (const ComplexRowVector& a, double s);
-  friend ComplexRowVector operator - (const ComplexRowVector& a, double s);
-  friend ComplexRowVector operator * (const ComplexRowVector& a, double s);
-  friend ComplexRowVector operator / (const ComplexRowVector& a, double s);
-
-// scalar by row vector -> row vector operations
-
-  friend ComplexRowVector operator + (double s, const ComplexRowVector& a);
-  friend ComplexRowVector operator - (double s, const ComplexRowVector& a);
-  friend ComplexRowVector operator * (double s, const ComplexRowVector& a);
-  friend ComplexRowVector operator / (double s, const ComplexRowVector& a);
-
-// row vector by column vector -> scalar
-
-  friend Complex operator * (const ComplexRowVector& a, const ColumnVector& b);
-
-  friend Complex operator * (const ComplexRowVector& a,
-			     const ComplexColumnVector& b);
-
-// row vector by matrix -> row vector
-
-  friend ComplexRowVector operator * (const ComplexRowVector& a,
-				      const ComplexMatrix& b);
-
-// row vector by row vector -> row vector operations
-
-  friend ComplexRowVector operator + (const ComplexRowVector& a,
-				      const RowVector& b);
-  friend ComplexRowVector operator - (const ComplexRowVector& a,
-				      const RowVector& b);
-
-  friend ComplexRowVector product (const ComplexRowVector& a,
-				   const RowVector& b);
-  friend ComplexRowVector quotient (const ComplexRowVector& a,
-				    const RowVector& b);
-
-// other operations
-
-  friend ComplexRowVector map (c_c_Mapper f, const ComplexRowVector& a);
-  friend RowVector map (d_c_Mapper f, const ComplexRowVector& a);
-  void map (c_c_Mapper f);
-
-  Complex min (void) const;
-  Complex max (void) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const ComplexRowVector& a);
-  friend ostream& operator >> (ostream& is, ComplexRowVector& a);
-
-#define KLUDGE_VECTORS
-#define TYPE Complex
-#define KL_VEC_TYPE ComplexRowVector
-#include "mx-kludge.h"
-#undef KLUDGE_VECTORS
-#undef TYPE
-#undef KL_VEC_TYPE
-
-private:
-
-  ComplexRowVector (Complex *d, int l) : Array<Complex> (d, l) { }
-};
-
-/*
- * Complex Diagonal Matrix class
- */
-
-class ComplexDiagMatrix : public DiagArray<Complex>
-{
-public:
-
-  ComplexDiagMatrix (void) : DiagArray<Complex> () { }
-  ComplexDiagMatrix (int n) : DiagArray<Complex> (n) { }
-  ComplexDiagMatrix (int n, const Complex& val)
-    : DiagArray<Complex> (n, val) { }
-  ComplexDiagMatrix (int r, int c) : DiagArray<Complex> (r, c) { }
-  ComplexDiagMatrix (int r, int c, const Complex& val)
-    : DiagArray<Complex> (r, c, val) { }
-  ComplexDiagMatrix (const RowVector& a);
-  ComplexDiagMatrix (const ComplexRowVector& a) : DiagArray<Complex> (a) { }
-  ComplexDiagMatrix (const ColumnVector& a);
-  ComplexDiagMatrix (const ComplexColumnVector& a)
-    : DiagArray<Complex> (a) { }
-  ComplexDiagMatrix (const DiagMatrix& a);
-  ComplexDiagMatrix (const DiagArray<Complex>& a)
-    : DiagArray<Complex> (a) { }
-  ComplexDiagMatrix (const ComplexDiagMatrix& a) : DiagArray<Complex> (a) { }
-//  ComplexDiagMatrix (const Complex& a) : DiagArray<Complex> (1, a) { }
-
-  ComplexDiagMatrix& operator = (const ComplexDiagMatrix& a)
-    {
-      DiagArray<Complex>::operator = (a);
-      return *this;
-    }
-
-//  operator DiagArray<Complex>& () const { return *this; }
-
-  int operator == (const ComplexDiagMatrix& a) const;
-  int operator != (const ComplexDiagMatrix& a) const;
-
-  ComplexDiagMatrix& fill (double val);
-  ComplexDiagMatrix& fill (const Complex& val);
-  ComplexDiagMatrix& fill (double val, int beg, int end);
-  ComplexDiagMatrix& fill (const Complex& val, int beg, int end);
-  ComplexDiagMatrix& fill (const ColumnVector& a);
-  ComplexDiagMatrix& fill (const ComplexColumnVector& a);
-  ComplexDiagMatrix& fill (const RowVector& a);
-  ComplexDiagMatrix& fill (const ComplexRowVector& a);
-  ComplexDiagMatrix& fill (const ColumnVector& a, int beg);
-  ComplexDiagMatrix& fill (const ComplexColumnVector& a, int beg);
-  ComplexDiagMatrix& fill (const RowVector& a, int beg);
-  ComplexDiagMatrix& fill (const ComplexRowVector& a, int beg);
-
-  ComplexDiagMatrix hermitian (void) const;  // complex conjugate transpose
-  ComplexDiagMatrix transpose (void) const;
-
-  friend DiagMatrix real (const ComplexDiagMatrix& a);
-  friend DiagMatrix imag (const ComplexDiagMatrix& a);
-  friend ComplexDiagMatrix conj (const ComplexDiagMatrix& a);
-
-// resize is the destructive analog for this one
-
-  ComplexMatrix extract (int r1, int c1, int r2, int c2) const;
+#include "mx-base.h"
 
-// extract row or column i.
-
-  ComplexRowVector row (int i) const;
-  ComplexRowVector row (char *s) const;
-
-  ComplexColumnVector column (int i) const;
-  ComplexColumnVector column (char *s) const;
-
-  ComplexDiagMatrix inverse (int& info) const;
-  ComplexDiagMatrix inverse (void) const;
-
-// diagonal matrix by diagonal matrix -> diagonal matrix operations
-
-  ComplexDiagMatrix& operator += (const DiagMatrix& a);
-  ComplexDiagMatrix& operator -= (const DiagMatrix& a);
-
-  ComplexDiagMatrix& operator += (const ComplexDiagMatrix& a);
-  ComplexDiagMatrix& operator -= (const ComplexDiagMatrix& a);
-
-// diagonal matrix by scalar -> matrix operations
-
-  friend ComplexMatrix operator + (const ComplexDiagMatrix& a, double s);
-  friend ComplexMatrix operator - (const ComplexDiagMatrix& a, double s);
-
-  friend ComplexMatrix operator + (const ComplexDiagMatrix& a,
-				   const Complex& s);
-  friend ComplexMatrix operator - (const ComplexDiagMatrix& a,
-				   const Complex& s);
-
-// diagonal matrix by scalar -> diagonal matrix operations
-
-  friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, double s);
-  friend ComplexDiagMatrix operator / (const ComplexDiagMatrix& a, double s);
-
-// scalar by diagonal matrix -> matrix operations
-
-  friend ComplexMatrix operator + (double s, const ComplexDiagMatrix& a);
-  friend ComplexMatrix operator - (double s, const ComplexDiagMatrix& a);
-
-  friend ComplexMatrix operator + (const Complex& s,
-				   const ComplexDiagMatrix& a);
-  friend ComplexMatrix operator - (const Complex& s,
-				   const ComplexDiagMatrix& a);
-
-// scalar by diagonal matrix -> diagonal matrix operations
-
-  friend ComplexDiagMatrix operator * (double s, const ComplexDiagMatrix& a);
-
-// diagonal matrix by column vector -> column vector operations
-
-  friend ComplexColumnVector operator * (const ComplexDiagMatrix& a,
-					 const ColumnVector& b);
-
-  friend ComplexColumnVector operator * (const ComplexDiagMatrix& a,
-					 const ComplexColumnVector& b);
-
-// diagonal matrix by diagonal matrix -> diagonal matrix operations
-
-  friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a,
-				       const ComplexDiagMatrix& b);
-
-  friend ComplexDiagMatrix operator + (const ComplexDiagMatrix& a,
-				       const DiagMatrix& b);
-  friend ComplexDiagMatrix operator - (const ComplexDiagMatrix& a,
-				       const DiagMatrix& b);
-  friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a,
-				       const DiagMatrix& b);
-
-  friend ComplexDiagMatrix product (const ComplexDiagMatrix& a,
-				    const DiagMatrix& b); 
-
-// diagonal matrix by matrix -> matrix operations
-
-  friend ComplexMatrix operator + (const ComplexDiagMatrix& a,
-				   const Matrix& b); 
-  friend ComplexMatrix operator - (const ComplexDiagMatrix& a,
-				   const Matrix& b);
-  friend ComplexMatrix operator * (const ComplexDiagMatrix& a,
-				   const Matrix& b);
-
-  friend ComplexMatrix operator + (const ComplexDiagMatrix& a,
-				   const ComplexMatrix& b);
-  friend ComplexMatrix operator - (const ComplexDiagMatrix& a,
-				   const ComplexMatrix& b);
-  friend ComplexMatrix operator * (const ComplexDiagMatrix& a,
-				   const ComplexMatrix& b);
-
-// other operations
-
-  ComplexColumnVector diag (void) const;
-  ComplexColumnVector diag (int k) const;
-
-// i/o
-
-  friend ostream& operator << (ostream& os, const ComplexDiagMatrix& a);
-
-#define KLUDGE_DIAG_MATRICES
-#define TYPE Complex
-#define KL_DMAT_TYPE ComplexDiagMatrix
-#include "mx-kludge.h"
-#undef KLUDGE_DIAG_MATRICES
-#undef TYPE
-#undef KL_DMAT_TYPE
-
-private:
-
-  ComplexDiagMatrix (Complex *d, int nr, int nc)
-    : DiagArray<Complex> (d, nr, nc) { }
-};
-
-/*
- * Result of a AEP Balance operation
- */
-
-class AEPBALANCE
-{
-friend class Matrix;
-
-public:
-
-  AEPBALANCE (void) {}
-
-  AEPBALANCE (const Matrix& a, const char *balance_job);
-
-  AEPBALANCE (const AEPBALANCE& a);
-
-  AEPBALANCE& operator = (const AEPBALANCE& a);
-  Matrix balanced_matrix (void) const;
-  Matrix balancing_matrix (void) const;
-  friend ostream& operator << (ostream& os, const AEPBALANCE& a);
-
-private:
-
-  int init (const Matrix& a, const char * balance_job);
-
-  Matrix balanced_mat;
-  Matrix balancing_mat;
-};
-
-inline AEPBALANCE::AEPBALANCE (const Matrix& a,const char * balance_job) 
-{
-  init (a, balance_job); 
-}
-
-inline AEPBALANCE::AEPBALANCE (const AEPBALANCE& a)
-{
-  balanced_mat = a.balanced_mat;
-  balancing_mat = a.balancing_mat;
-}
-
-inline AEPBALANCE&
-AEPBALANCE::operator = (const AEPBALANCE& a)
-{
-  balanced_mat = a.balanced_mat;
-  balancing_mat = a.balancing_mat;
-
-  return *this;
-}
-
-inline Matrix AEPBALANCE::balanced_matrix (void) const
-{
-  return balanced_mat;
-}
-
-inline Matrix AEPBALANCE::balancing_matrix (void) const
-{
-  return balancing_mat;
-}
-
-/*
- * Result of a Complex balancing operation
- */
-
-class ComplexAEPBALANCE
-{
-friend class ComplexMatrix;
-
-public:
-
-  ComplexAEPBALANCE (void) {}
-  ComplexAEPBALANCE (const ComplexMatrix& a, const char *balance_job);
-  ComplexAEPBALANCE (const ComplexAEPBALANCE& a);
-  ComplexAEPBALANCE& operator = (const ComplexAEPBALANCE& a);
-  ComplexMatrix balanced_matrix (void) const;
-  ComplexMatrix balancing_matrix (void) const;
-
-  friend ostream& operator << (ostream& os, const ComplexAEPBALANCE& a);
-
-private:
-
-  int init (const ComplexMatrix& a, const char * balance_job);
-
-  ComplexMatrix balanced_mat;
-  ComplexMatrix balancing_mat;
-};
-
-inline ComplexAEPBALANCE::ComplexAEPBALANCE (const ComplexMatrix& a,
-					     const char * balance_job)
-{
-  init (a, balance_job); 
-}
-
-inline ComplexAEPBALANCE::ComplexAEPBALANCE (const ComplexAEPBALANCE& a)
-{
-  balanced_mat = a.balanced_mat;
-  balancing_mat = a.balancing_mat;
-}
-
-inline ComplexAEPBALANCE&
-ComplexAEPBALANCE::operator = (const ComplexAEPBALANCE& a)
-{
-  balanced_mat = a.balanced_mat;
-  balancing_mat = a.balancing_mat;
-
-  return *this;
-}
-
-inline ComplexMatrix ComplexAEPBALANCE::balanced_matrix (void) const
-{
-  return balanced_mat;
-}
-
-inline ComplexMatrix ComplexAEPBALANCE::balancing_matrix (void) const
-{
-  return balancing_mat;
-}
-
-/*
- * Result of a Determinant calculation.
- */
-
-class DET
-{
-public:
-
-  DET (void) {}
-
-  DET (const DET& a);
-
-  DET& operator = (const DET& a);
-
-  int value_will_overflow (void) const;
-  int value_will_underflow (void) const;
-  double coefficient (void) const;
-  int exponent (void) const;
-  double value (void) const;
-
-  friend ostream&  operator << (ostream& os, const DET& a);
-
-private:
-
-  DET (const double *d);
-
-  double det [2];
-};
-
-inline DET::DET (const DET& a)
-{
-  det[0] = a.det[0];
-  det[1] = a.det[1];
-}
-
-inline DET& DET::operator = (const DET& a)
-{
-  det[0] = a.det[0];
-  det[1] = a.det[1];
-  return *this;
-}
-
-inline DET::DET (const double *d)
-{
-  det[0] = d[0];
-  det[1] = d[1];
-}
-
-/*
- * Result of a Determinant calculation.
- */
-
-class ComplexDET
-{
-public:
-
-  ComplexDET (void) {}
-
-  ComplexDET (const ComplexDET& a);
-
-  ComplexDET& operator = (const ComplexDET& a);
-
-  int value_will_overflow (void) const;
-  int value_will_underflow (void) const;
-  Complex coefficient (void) const;
-  int exponent (void) const;
-  Complex value (void) const;
-
-  friend ostream&  operator << (ostream& os, const ComplexDET& a);
-
-private:
-
-  ComplexDET (const Complex *d);
-
-  Complex det [2];
-};
-
-inline ComplexDET::ComplexDET (const ComplexDET& a)
-{
-  det[0] = a.det[0];
-  det[1] = a.det[1];
-}
-
-inline ComplexDET& ComplexDET::operator = (const ComplexDET& a)
-{
-  det[0] = a.det[0];
-  det[1] = a.det[1];
-  return *this;
-}
-
-inline ComplexDET::ComplexDET (const Complex *d)
-{
-  det[0] = d[0];
-  det[1] = d[1];
-}
-
-/*
- * Result of a GEP Balance operation
- * Note: currenlty only do balancing on real data.  Complex balancing
- * done on magnitudes of complex data.
- */
-
-class GEPBALANCE
-{
-friend class Matrix;
-
-public:
-
-  GEPBALANCE (void) {}
-
-  GEPBALANCE (const Matrix& a, const Matrix &, const char *balance_job);
-
-  GEPBALANCE (const GEPBALANCE& a);
-
-  GEPBALANCE& operator = (const GEPBALANCE& a);
-  Matrix balanced_a_matrix (void) const;
-  Matrix balanced_b_matrix (void) const;
-  Matrix left_balancing_matrix (void) const;
-  Matrix right_balancing_matrix (void) const;
-  friend ostream& operator << (ostream& os, const GEPBALANCE& a);
-
-private:
-
-  int init (const Matrix& a, const Matrix& b, const char * balance_job);
-
-  Matrix balanced_a_mat;
-  Matrix balanced_b_mat;
-  Matrix left_balancing_mat;
-  Matrix right_balancing_mat;
-};
-
-inline GEPBALANCE::GEPBALANCE (const Matrix& a, const Matrix& b, 
-			       const char * balance_job) 
-{
-  init (a, b, balance_job); 
-}
-
-inline GEPBALANCE::GEPBALANCE (const GEPBALANCE& a)
-{
-  balanced_a_mat = a.balanced_a_mat;
-  balanced_b_mat = a.balanced_b_mat;
-  left_balancing_mat = a.left_balancing_mat;
-  right_balancing_mat = a.right_balancing_mat;
-}
-
-inline GEPBALANCE&
-GEPBALANCE::operator = (const GEPBALANCE& a)
-{
-  balanced_a_mat = a.balanced_a_mat;
-  balanced_b_mat = a.balanced_b_mat;
-  left_balancing_mat = a.left_balancing_mat;
-  right_balancing_mat = a.right_balancing_mat;
-
-  return *this;
-}
-
-inline Matrix GEPBALANCE::balanced_a_matrix (void) const 
-{
-  return balanced_a_mat;
-}
-
-inline Matrix GEPBALANCE::balanced_b_matrix (void) const 
-{
-  return balanced_b_mat;
-}
-
-inline Matrix GEPBALANCE::left_balancing_matrix (void) const 
-{
-  return left_balancing_mat;
-}
-
-inline Matrix GEPBALANCE::right_balancing_matrix (void) const 
-{
-  return right_balancing_mat;
-}
-
-/*
- * Result of a Cholesky Factorization
- */
-
-class CHOL
-{
-friend class Matrix;
-
-public:
-
-  CHOL (void) {}
-
-  CHOL (const Matrix& a);
-  CHOL (const Matrix& a, int& info);
-
-  CHOL (const CHOL& a);
-
-  CHOL& operator = (const CHOL& a);
-  Matrix chol_matrix (void) const;
-  friend ostream& operator << (ostream& os, const CHOL& a);
-
-private:
-
-  int init (const Matrix& a);
-
-  Matrix chol_mat;
-};
-
-inline CHOL::CHOL (const Matrix& a)
-{
-  init (a);
-}
-
-inline CHOL::CHOL (const Matrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline CHOL::CHOL (const CHOL& a)
-{
-  chol_mat = a.chol_mat;
-}
-
-inline CHOL&
-CHOL::operator = (const CHOL& a)
-{
-  chol_mat = a.chol_mat;
-
-  return *this;
-}
-
-inline Matrix CHOL::chol_matrix (void) const { return chol_mat; }
-
-/*
- * Result of a Cholesky Factorization
- */
-
-class ComplexCHOL
-{
-friend class ComplexMatrix;
-
-public:
-
-  ComplexCHOL (void) {}
-  ComplexCHOL (const ComplexMatrix& a);
-  ComplexCHOL (const ComplexMatrix& a, int& info);
-  ComplexCHOL (const ComplexCHOL& a);
-  ComplexCHOL& operator = (const ComplexCHOL& a);
-  ComplexMatrix chol_matrix (void) const;
-
-  friend ostream& operator << (ostream& os, const ComplexCHOL& a);
-
-private:
-
-  int init (const ComplexMatrix& a);
-
-  ComplexMatrix chol_mat;
-};
-
-inline ComplexCHOL::ComplexCHOL (const ComplexMatrix& a)
-{
-  init (a);
-}
-
-inline ComplexCHOL::ComplexCHOL (const ComplexMatrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline ComplexCHOL::ComplexCHOL (const ComplexCHOL& a)
-{
-  chol_mat = a.chol_mat;
-}
-
-inline ComplexCHOL&
-ComplexCHOL::operator = (const ComplexCHOL& a)
-{
-  chol_mat = a.chol_mat;
-
-  return *this;
-}
-
-inline ComplexMatrix ComplexCHOL::chol_matrix (void) const
-{
-  return chol_mat;
-}
-
-  
-/*
- * Result of a Hessenberg Decomposition
- */
-
-class HESS
-{
-friend class Matrix;
-
-public:
-
-  HESS (void) {}
-
-  HESS (const Matrix& a);
-  HESS (const Matrix&a, int& info);
-
-  HESS (const HESS& a);
-
-  HESS& operator = (const HESS& a);
-  Matrix hess_matrix (void) const;
-  Matrix unitary_hess_matrix (void) const;
-  friend ostream& operator << (ostream& os, const HESS& a);
-
-private:
-
-  int init (const Matrix& a);
-
-  Matrix hess_mat;
-  Matrix unitary_hess_mat;
-};
-
-inline HESS::HESS (const Matrix& a)
-{
-  init (a);
-}
-
-inline HESS::HESS (const Matrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline HESS::HESS (const HESS& a)
-{
-  hess_mat = a.hess_mat;
-  unitary_hess_mat = a.unitary_hess_mat;
-}
-
-inline HESS&
-HESS::operator = (const HESS& a)
-{
-  hess_mat = a.hess_mat;
-  unitary_hess_mat = a.unitary_hess_mat;
-
-  return *this;
-}
-
-inline Matrix HESS::hess_matrix (void) const
-{
-  return hess_mat;
-}
-
-inline Matrix HESS::unitary_hess_matrix (void) const
-{
-  return unitary_hess_mat;
-}
-
-/*
- * Result of a Hessenberg Decomposition
- */
-
-class ComplexHESS
-{
-friend class ComplexMatrix;
-
-public:
-
-  ComplexHESS (void) {}
-  ComplexHESS (const ComplexMatrix& a);
-  ComplexHESS (const ComplexMatrix& a, int& info);
-  ComplexHESS (const ComplexHESS& a);
-  ComplexHESS& operator = (const ComplexHESS& a);
-  ComplexMatrix hess_matrix (void) const;
-  ComplexMatrix unitary_hess_matrix (void) const;
-
-  friend ostream& operator << (ostream& os, const ComplexHESS& a);
-
-private:
-
-  int init (const ComplexMatrix& a);
-
-  ComplexMatrix hess_mat;
-  ComplexMatrix unitary_hess_mat;
-};
-
-inline ComplexHESS::ComplexHESS (const ComplexMatrix& a)
-{
-  init (a);
-}
-
-inline ComplexHESS::ComplexHESS (const ComplexMatrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline ComplexHESS::ComplexHESS (const ComplexHESS& a)
-{
-  hess_mat = a.hess_mat;
-  unitary_hess_mat = a.unitary_hess_mat;
-}
-
-inline ComplexHESS&
-ComplexHESS::operator = (const ComplexHESS& a)
-{
-  hess_mat = a.hess_mat;
-  unitary_hess_mat = a.unitary_hess_mat;
-
-  return *this;
-}
-
-inline ComplexMatrix ComplexHESS::hess_matrix (void) const
-{
-  return hess_mat;
-}
-
-inline ComplexMatrix ComplexHESS::unitary_hess_matrix (void) const
-{
-  return unitary_hess_mat;
-}
-
-/*
- * Result of a Schur Decomposition
- */
-
-class SCHUR
-{
-friend class Matrix;
-
-public:
-
-  SCHUR (void) {}
-
-  SCHUR (const Matrix& a, const char *ord);
-  SCHUR (const Matrix& a, const char *ord, int& info);
-
-  SCHUR (const SCHUR& a, const char *ord);
-
-  SCHUR& operator = (const SCHUR& a);
-
-  Matrix schur_matrix (void) const;
-  Matrix unitary_matrix (void) const;
-
-  friend ostream& operator << (ostream& os, const SCHUR& a);
-
-private:
-
-  int init (const Matrix& a, const char *ord);
-
-  Matrix schur_mat;
-  Matrix unitary_mat;
-};
-
-inline SCHUR::SCHUR (const Matrix& a, const char *ord)
-{
-  init (a, ord);
-}
-
-inline SCHUR::SCHUR (const Matrix& a, const char *ord, int& info) 
-{
-  info = init (a, ord);
-}
-
-inline SCHUR::SCHUR (const SCHUR& a, const char *ord)
-{
-  schur_mat = a.schur_mat;
-  unitary_mat = a.unitary_mat;
-}
-
-inline SCHUR&
-SCHUR::operator = (const SCHUR& a)
-{
-  schur_mat = a.schur_mat;
-  unitary_mat = a.unitary_mat;
-  
-  return *this;
-}
-
-inline Matrix SCHUR::schur_matrix (void) const
-{
-  return schur_mat;
-}
-
-inline Matrix SCHUR::unitary_matrix (void) const
-{
-  return unitary_mat;
-}
-
-/*
- * Result of a Schur Decomposition
- */
-
-class ComplexSCHUR
-{
-friend class ComplexMatrix;
-
-public:
-
-  ComplexSCHUR (void) {}
-
-  ComplexSCHUR (const ComplexMatrix& a, const char *ord);
-  ComplexSCHUR (const ComplexMatrix& a, const char *ord, int& info);
-
-  ComplexSCHUR (const ComplexSCHUR& a, const char *ord);
-
-  ComplexSCHUR& operator = (const ComplexSCHUR& a);
-
-  ComplexMatrix schur_matrix (void) const;
-  ComplexMatrix unitary_matrix (void) const;
-
-  friend ostream& operator << (ostream& os, const ComplexSCHUR& a);
-
-private:
-
-  int init (const ComplexMatrix& a, const char *ord);
-
-  ComplexMatrix schur_mat;
-  ComplexMatrix unitary_mat;
-};
-
-inline ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& a, const char *ord) 
-{
-  init (a,ord);
-}
-
-inline ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& a, const char *ord,
-				   int& info)
-{
-  info = init (a,ord);
-}
-
-inline ComplexSCHUR::ComplexSCHUR (const ComplexSCHUR& a, const char *ord)
-{
-  schur_mat = a.schur_mat;
-  unitary_mat = a.unitary_mat;
-}
-
-inline ComplexSCHUR&
-ComplexSCHUR::operator = (const ComplexSCHUR& a)
-{
-  schur_mat = a.schur_mat;
-  unitary_mat = a.unitary_mat;
-
-  return *this;
-}
-
-inline ComplexMatrix ComplexSCHUR::schur_matrix (void) const
-{
-  return schur_mat;
-}
-
-inline ComplexMatrix ComplexSCHUR::unitary_matrix (void) const
-{
-  return unitary_mat;
-}
-
-
-/*
- * Result of a Singular Value Decomposition.
- */
-
-class SVD
-{
-friend class Matrix;
-
-public:
-
-  SVD (void) {}
-
-  SVD (const Matrix& a);
-  SVD (const Matrix& a, int& info);
-
-  SVD (const SVD& a);
-
-  SVD& operator = (const SVD& a);
-
-  DiagMatrix singular_values (void) const;
-  Matrix left_singular_matrix (void) const;
-  Matrix right_singular_matrix (void) const;
-
-  friend ostream&  operator << (ostream& os, const SVD& a);
-
-private:
-
-  int init (const Matrix& a);
-
-  DiagMatrix sigma;
-  Matrix left_sm;
-  Matrix right_sm;
-};
-
-inline SVD::SVD (const Matrix& a)
-{
-  init (a);
-}
-
-inline SVD::SVD (const Matrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline SVD::SVD (const SVD& a)
-{
-  sigma = a.sigma;
-  left_sm = a.left_sm;
-  right_sm = a.right_sm;
-}
-
-inline SVD&
-SVD::operator = (const SVD& a)
-{
-  sigma = a.sigma;
-  left_sm = a.left_sm;
-  right_sm = a.right_sm;
-
-  return *this;
-}
-
-inline DiagMatrix SVD::singular_values (void) const
-{
-  return sigma;
-}
-
-inline Matrix SVD::left_singular_matrix (void) const
-{
-  return left_sm;
-}
-
-inline Matrix SVD::right_singular_matrix (void) const
-{
-  return right_sm;
-}
-
-/*
- * Result of a Singular Value Decomposition.
- */
-
-class ComplexSVD
-{
-friend class ComplexMatrix;
-
-public:
-
-  ComplexSVD (void) {}
-
-  ComplexSVD (const ComplexMatrix& a);
-  ComplexSVD (const ComplexMatrix& a, int& info);
-
-  ComplexSVD (const ComplexSVD& a);
-
-  ComplexSVD& operator = (const ComplexSVD& a);
-
-  DiagMatrix singular_values (void) const;
-  ComplexMatrix left_singular_matrix (void) const;
-  ComplexMatrix right_singular_matrix (void) const;
-
-  friend ostream&  operator << (ostream& os, const ComplexSVD& a);
-
-private:
-
-  int init (const ComplexMatrix& a);
-
-  DiagMatrix sigma;
-  ComplexMatrix left_sm;
-  ComplexMatrix right_sm;
-};
-
-inline ComplexSVD::ComplexSVD (const ComplexMatrix& a)
-{
-  init (a);
-}
-
-inline ComplexSVD::ComplexSVD (const ComplexMatrix& a, int& info)
-{
-  info = init (a);
-} 
-
-inline ComplexSVD::ComplexSVD (const ComplexSVD& a)
-{
-  sigma = a.sigma;
-  left_sm = a.left_sm;
-  right_sm = a.right_sm;
-}
-
-inline ComplexSVD&
-ComplexSVD::operator = (const ComplexSVD& a)
-{
-  sigma = a.sigma;
-  left_sm = a.left_sm;
-  right_sm = a.right_sm;
-
-  return *this;
-}
-
-inline DiagMatrix ComplexSVD::singular_values (void) const
-{
-  return sigma;
-}
-
-inline ComplexMatrix ComplexSVD::left_singular_matrix (void) const
-{
-  return left_sm;
-}
-
-inline ComplexMatrix ComplexSVD::right_singular_matrix (void) const
-{
-  return right_sm;
-}
-
-/*
- * Result of an Eigenvalue computation.
- */
-
-class EIG
-{
-friend class Matrix;
-friend class ComplexMatrix;
-
-public:
-
-  EIG (void) {}
-
-  EIG (const Matrix& a);
-  EIG (const Matrix& a, int& info);
-
-  EIG (const ComplexMatrix& a);
-  EIG (const ComplexMatrix& a, int& info);
-
-  EIG (const EIG& a);
-
-  EIG& operator = (const EIG& a);
-
-  ComplexColumnVector eigenvalues (void) const;
-  ComplexMatrix eigenvectors (void) const;
-
-  friend ostream&  operator << (ostream& os, const EIG& a);
-
-private:
-
-  int init (const Matrix& a);
-  int init (const ComplexMatrix& a);
-
-  ComplexColumnVector lambda;
-  ComplexMatrix v;
-};
-
-inline EIG::EIG (const Matrix& a)
-{
-  init (a);
-}
-
-inline EIG::EIG (const Matrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline EIG::EIG (const ComplexMatrix& a)
-{
-  init (a);
-}
-
-inline EIG::EIG (const ComplexMatrix& a, int& info)
-{
-  info = init (a);
-}
-
-inline EIG::EIG (const EIG& a)
-{
-  lambda = a.lambda;
-  v = a.v;
-}
-
-inline EIG& EIG::operator = (const EIG& a)
-{
-  lambda = a.lambda;
-  v = a.v;
-  return *this;
-}
-
-inline ComplexColumnVector EIG::eigenvalues (void) const
-{
-  return lambda;
-}
-
-inline ComplexMatrix EIG::eigenvectors (void) const
-{
-  return v;
-}
-
-/*
- * Result of an LU decomposition.
- */
-
-class LU
-{
-friend class Matrix;
-
-public:
-
-  LU (void) {}
-
-  LU (const Matrix& a);
-
-  LU (const LU& a);
-
-  LU& operator = (const LU& a);
-
-  Matrix L (void) const;
-  Matrix U (void) const;
-  Matrix P (void) const;
-
-  friend ostream&  operator << (ostream& os, const LU& a);
-
-private:
-
-  Matrix l;
-  Matrix u;
-  Matrix p;
-};
-
-inline LU::LU (const LU& a)
-{
-  l = a.l;
-  u = a.u;
-  p = a.p;
-}
-
-inline LU& LU::operator = (const LU& a)
-{
-  l = a.l;
-  u = a.u;
-  p = a.p;
-  return *this;
-}
-
-inline Matrix LU::L (void) const
-{
-  return l;
-}
-
-inline Matrix LU::U (void) const
-{
-  return u;
-}
-
-inline Matrix LU::P (void) const
-{
-  return p;
-}
-
-class ComplexLU
-{
-friend class ComplexMatrix;
-
-public:
-
-  ComplexLU (void) {}
-
-  ComplexLU (const ComplexMatrix& a);
-
-  ComplexLU (const ComplexLU& a);
-
-  ComplexLU& operator = (const ComplexLU& a);
-
-  ComplexMatrix L (void) const;
-  ComplexMatrix U (void) const;
-  Matrix P (void) const;
-
-  friend ostream&  operator << (ostream& os, const ComplexLU& a);
-
-private:
-
-  ComplexMatrix l;
-  ComplexMatrix u;
-  Matrix p;
-};
-
-inline ComplexLU::ComplexLU (const ComplexLU& a)
-{
-  l = a.l;
-  u = a.u;
-  p = a.p;
-}
-
-inline ComplexLU& ComplexLU::operator = (const ComplexLU& a)
-{
-  l = a.l;
-  u = a.u;
-  p = a.p;
-  return *this;
-}
-
-inline ComplexMatrix ComplexLU::L (void) const
-{
-  return l;
-}
-
-inline ComplexMatrix ComplexLU::U (void) const
-{
-  return u;
-}
-
-inline Matrix ComplexLU::P (void) const
-{
-  return p;
-}
-
-/*
- * Result of a QR decomposition.
- */
-
-class QR
-{
-public:
-
-  QR (void) {}
-
-  QR (const Matrix& A);
-
-  QR (const QR& a);
-
-  QR& operator = (const QR& a);
-
-  Matrix Q (void) const;
-  Matrix R (void) const;
-
-  friend ostream&  operator << (ostream& os, const QR& a);
-
-private:
-
-  Matrix q;
-  Matrix r;
-};
-
-inline QR::QR (const QR& a)
-{
-  q = a.q;
-  r = a.r;
-}
-
-inline QR& QR::operator = (const QR& a)
-{
-  q = a.q;
-  r = a.r;
-  return *this;
-}
-
-inline Matrix QR::Q (void) const
-{
-  return q;
-}
-
-inline Matrix QR::R (void) const
-{
-  return r;
-}
-
-class ComplexQR
-{
-public:
-
-  ComplexQR (void) {}
-
-  ComplexQR (const ComplexMatrix& A);
-
-  ComplexQR (const ComplexQR& a);
-
-  ComplexQR& operator = (const ComplexQR& a);
-
-  ComplexMatrix Q (void) const;
-  ComplexMatrix R (void) const;
-
-  friend ostream&  operator << (ostream& os, const ComplexQR& a);
-
-private:
-
-  ComplexMatrix q;
-  ComplexMatrix r;
-};
-
-inline ComplexQR::ComplexQR (const ComplexQR& a)
-{
-  q = a.q;
-  r = a.r;
-}
-
-inline ComplexQR& ComplexQR::operator = (const ComplexQR& a)
-{
-  q = a.q;
-  r = a.r;
-  return *this;
-}
-
-inline ComplexMatrix ComplexQR::Q (void) const
-{
-  return q;
-}
-
-inline ComplexMatrix ComplexQR::R (void) const
-{
-  return r;
-}
-
-} // extern "C++"
+#include "mx-ext.h"
 
 #endif
 
diff --git a/liboctave/NLFunc.cc b/liboctave/NLFunc.cc
--- a/liboctave/NLFunc.cc
+++ b/liboctave/NLFunc.cc
@@ -1,7 +1,7 @@
 // NLFunc.cc                                             -*- C++ -*-
 /*
 
-Copyright (C) 1992, 1993 John W. Eaton
+Copyright (C) 1992, 1993, 1994 John W. Eaton
 
 This file is part of Octave.
 
@@ -25,18 +25,22 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include "NLFunc.h"
 
 NLFunc::NLFunc (void)
 {
-  fun = NULL;
-  jac = NULL;
+  fun = 0;
+  jac = 0;
 }
 
 NLFunc::NLFunc (const nonlinear_fcn f)
 {
   fun = f;
-  jac = NULL;
+  jac = 0;
 }
 
 NLFunc::NLFunc (const nonlinear_fcn f, const jacobian_fcn j)
diff --git a/liboctave/NLFunc.h b/liboctave/NLFunc.h
--- a/liboctave/NLFunc.h
+++ b/liboctave/NLFunc.h
@@ -24,7 +24,12 @@
 #if !defined (octave_NLFunc_h)
 #define octave_NLFunc_h 1
 
-#include "Matrix.h"
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+class ColumnVector;
+class Matrix;
 
 extern "C++" {
 
diff --git a/liboctave/NPSOL.cc b/liboctave/NPSOL.cc
--- a/liboctave/NPSOL.cc
+++ b/liboctave/NPSOL.cc
@@ -25,6 +25,10 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include <math.h>
 #include <float.h>
 #include <string.h>
@@ -32,6 +36,8 @@
 
 #ifndef NPSOL_MISSING
 
+#include "dMatrix.h"
+
 #include "NPSOL.h"
 #include "f77-uscore.h"
 #include "sun-utils.h"
diff --git a/liboctave/NPSOL.h b/liboctave/NPSOL.h
--- a/liboctave/NPSOL.h
+++ b/liboctave/NPSOL.h
@@ -24,9 +24,13 @@
 #if !defined (octave_NPSOL_h)
 #define octave_NPSOL_h 1
 
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
 #ifndef NPSOL_MISSING
 
-#include "Matrix.h"
+#include "dColVector.h"
 #include "NLP.h"
 
 extern "C++" {
diff --git a/liboctave/ODE.h b/liboctave/ODE.h
--- a/liboctave/ODE.h
+++ b/liboctave/ODE.h
@@ -24,9 +24,14 @@
 #if !defined (octave_ODE_h)
 #define octave_ODE_h 1
 
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
 class ostream;
 
-#include "Matrix.h"
+#include "dMatrix.h"
+#include "dColVector.h"
 #include "ODEFunc.h"
 
 extern "C++" {
diff --git a/liboctave/Objective.cc b/liboctave/Objective.cc
--- a/liboctave/Objective.cc
+++ b/liboctave/Objective.cc
@@ -1,7 +1,7 @@
 // Objective.cc                                             -*- C++ -*-
 /*
 
-Copyright (C) 1992, 1993 John W. Eaton
+Copyright (C) 1992, 1993, 1994 John W. Eaton
 
 This file is part of Octave.
 
@@ -25,18 +25,22 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include "Objective.h"
 
 Objective::Objective (void)
 {
-  phi = (objective_fcn) NULL;
-  grad = (gradient_fcn) NULL;
+  phi = 0;
+  grad = 0;
 }
 
 Objective::Objective (const objective_fcn obj)
 {
   phi = obj;
-  grad = (gradient_fcn) NULL;
+  grad = 0;
 }
 
 Objective::Objective (const objective_fcn obj, const gradient_fcn g)
diff --git a/liboctave/Objective.h b/liboctave/Objective.h
--- a/liboctave/Objective.h
+++ b/liboctave/Objective.h
@@ -24,7 +24,11 @@
 #if !defined (octave_Objective_h)
 #define octave_Objective_h 1
 
-#include "Matrix.h"
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include "dColVector.h"
 
 extern "C++" {
 
diff --git a/liboctave/QLD.cc b/liboctave/QLD.cc
--- a/liboctave/QLD.cc
+++ b/liboctave/QLD.cc
@@ -1,7 +1,7 @@
 // QLD.cc                                                -*- C++ -*-
 /*
 
-Copyright (C) 1992, 1993 John W. Eaton
+Copyright (C) 1992, 1993, 1994 John W. Eaton
 
 This file is part of Octave.
 
@@ -25,8 +25,16 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include <math.h>
+#include <iostream.h>
 
+#include "dMatrix.h"
+#include "dColVector.h"
+#include "dRowVector.h"
 #include "QLD.h"
 #include "f77-uscore.h"
 
@@ -78,15 +86,15 @@
   cout << H;
   cout << c;
 
-  double *pc = (double *) NULL;
+  double *pc = 0;
   if (c.capacity () > 0)
     pc = c.fortran_vec ();
 
-  double *pa = (double *) NULL;
+  double *pa = 0;
   if (A1.rows () > 0 && A1.columns () > 0)
     pa = A1.fortran_vec ();
 
-  double *pb = (double *) NULL;
+  double *pb = 0;
   if (b1.capacity () > 0)
     pb = b1.fortran_vec ();
 
diff --git a/liboctave/QLD.h b/liboctave/QLD.h
--- a/liboctave/QLD.h
+++ b/liboctave/QLD.h
@@ -24,7 +24,13 @@
 #if !defined (octave_QLD_h)
 #define octave_QLD_h 1
 
-#include "Matrix.h"
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+class Matrix;
+class ColumnVector;
+
 #include "QP.h"
 
 extern "C++" {
diff --git a/liboctave/Quad.cc b/liboctave/Quad.cc
--- a/liboctave/Quad.cc
+++ b/liboctave/Quad.cc
@@ -25,6 +25,10 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include <math.h>
 #include <float.h>
 
diff --git a/liboctave/Quad.h b/liboctave/Quad.h
--- a/liboctave/Quad.h
+++ b/liboctave/Quad.h
@@ -24,7 +24,11 @@
 #if !defined (octave_Quad_h)
 #define octave_Quad_h 1
 
-#include "Matrix.h"
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include "dColVector.h"
 
 extern "C++" {
 
diff --git a/liboctave/Range.cc b/liboctave/Range.cc
--- a/liboctave/Range.cc
+++ b/liboctave/Range.cc
@@ -25,6 +25,10 @@
 #include "config.h"
 #endif
 
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
 #include <iostream.h>
 #include <limits.h>
 
diff --git a/liboctave/Range.h b/liboctave/Range.h
--- a/liboctave/Range.h
+++ b/liboctave/Range.h
@@ -24,6 +24,10 @@
 #if !defined (octave_Range_h)
 #define octave_Range_h 1
 
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
 extern "C++" {
 
 class istream;
diff --git a/liboctave/mx-inlines.cc b/liboctave/mx-inlines.cc
--- a/liboctave/mx-inlines.cc
+++ b/liboctave/mx-inlines.cc
@@ -1,7 +1,7 @@
 // Helper functions for matrix classes.                 -*- C++ -*-
 /*
 
-Copyright (C) 1992, 1993 John W. Eaton
+Copyright (C) 1992, 1993, 1994 John W. Eaton
 
 This file is part of Octave.
 
@@ -21,6 +21,8 @@
 
 */
 
+#include <Complex.h>
+
 // But first, some helper functions...
 
 static inline double *