Mercurial > hg > octave-max

changeset 1988:7b56630a1e05

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-03-02 00:33:22 by jwe] Initial revision
author jwe
date Sat, 02 Mar 1996 00:33:22 +0000
parents 6df7b42db205
children a4b0826e240c
files liboctave/Array2-idx.h liboctave/Array2.cc liboctave/Array2.h liboctave/Array3-idx.h liboctave/Array3.cc liboctave/Array3.h liboctave/DiagArray2.cc liboctave/DiagArray2.h liboctave/MArray-defs.h liboctave/MArray2.cc liboctave/MArray2.h liboctave/MDiagArray2.cc liboctave/MDiagArray2.h
diffstat 13 files changed, 2330 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
new file mode 100644
--- /dev/null
+++ b/liboctave/Array2-idx.h
@@ -0,0 +1,580 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#include "Array-flags.h"
+#include "idx-vector.h"
+#include "lo-error.h"
+
+template <class T>
+Array2<T>
+Array2<T>::value (void)
+{
+  Array2<T> retval;
+
+  int n_idx = index_count ();
+
+  int nr = d1;
+  int nc = d2;
+
+  if (n_idx == 2)
+    {
+      idx_vector *tmp = get_idx ();
+      idx_vector idx_i = tmp[0];
+      idx_vector idx_j = tmp[1];
+
+      int n = idx_i.freeze (nr, "row", liboctave_pzo_flag);
+      int m = idx_j.freeze (nc, "column", liboctave_pzo_flag);
+
+      if (idx_i && idx_j)
+	{
+	  if (n == 0)
+	    {
+	      if (m == 0 || idx_j.is_colon_equiv (nc, 1))
+		retval.resize (0, 0);
+	      else
+		(*current_liboctave_error_handler)
+		  ("invalid row index = 0");
+	    }
+	  else if (m == 0)
+	    {
+	      if (n == 0 || idx_i.is_colon_equiv (nr, 1))
+		retval.resize (0, 0);
+	      else
+		(*current_liboctave_error_handler)
+		  ("invalid column index = 0");
+	    }
+	  else if (idx_i.is_colon_equiv (nr) && idx_j.is_colon_equiv (nc))
+	    {
+	      retval = *this;
+	    }
+	  else
+	    {
+	      retval.resize (n, m);
+
+	      for (int j = 0; j < m; j++)
+		{
+		  int jj = idx_j.elem (j);
+		  for (int i = 0; i < n; i++)
+		    {
+		      int ii = idx_i.elem (i);
+		      retval.elem (i, j) = elem (ii, jj);
+		    }
+		}
+	    }
+	}
+      // idx_vector::freeze() printed an error message for us.
+    }
+  else if (n_idx == 1)
+    {
+      if (nr == 1 && nc == 1)
+	{
+	  Array<T> tmp = Array<T>::value ();
+
+	  int len = tmp.length ();
+
+	  if (len == 0)
+	    retval = Array2<T> (0, 0);
+	  else
+	    {
+	      if (liboctave_pcv_flag)
+		retval = Array2<T> (tmp, len, 1);
+	      else
+		retval = Array2<T> (tmp, 1, len);
+	    }
+	}
+      else if (nr == 1 || nc == 1)
+	{
+	  int result_is_column_vector = (nc == 1);
+
+	  if (liboctave_dfi_flag)
+	    {
+	      idx_vector *tmp = get_idx ();
+	      idx_vector idx = tmp[0];
+
+	      if (idx.is_colon ())
+		result_is_column_vector = 1;
+	    }
+
+	  Array<T> tmp = Array<T>::value ();
+
+	  int len = tmp.length ();
+
+	  if (len == 0)
+	    retval = Array2<T> (0, 0);
+	  else
+	    {
+	      if (result_is_column_vector)
+		retval = Array2<T> (tmp, len, 1);
+	      else
+		retval = Array2<T> (tmp, 1, len);
+	    }
+	}
+      else if (liboctave_dfi_flag)
+	{
+	  // This code is only for indexing matrices.  The vector
+	  // cases are handled above.
+
+	  idx_vector *tmp = get_idx ();
+	  idx_vector idx = tmp[0];
+
+	  idx.freeze (nr * nc, "matrix", liboctave_pzo_flag);
+
+	  if (idx)
+	    {
+	      int result_nr = idx.orig_rows ();
+	      int result_nc = idx.orig_columns ();
+
+	      if (idx.is_colon ())
+		{
+		  result_nr = nr * nc;
+		  result_nc = 1;
+		}
+	      else if (idx.one_zero_only ())
+		{
+		  result_nr = idx.ones_count ();
+		  result_nc = (result_nr > 0 ? 1 : 0);
+		}
+
+	      retval.resize (result_nr, result_nc);
+
+	      int k = 0;
+	      for (int j = 0; j < result_nc; j++)
+		{
+		  for (int i = 0; i < result_nr; i++)
+		    {
+		      int ii = idx.elem (k++);
+		      int fr = ii % nr;
+		      int fc = (ii - fr) / nr;
+		      retval.elem (i, j) = elem (fr, fc);
+		    }
+		}
+	    }
+	  // idx_vector::freeze() printed an error message for us.
+	}
+      else
+	(*current_liboctave_error_handler)
+	  ("single index only valid for row or column vector");
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("invalid number of indices for matrix expression");
+
+  clear_index ();
+
+  return retval;
+}
+
+template <class T>
+void
+Array2<T>::maybe_delete_elements (idx_vector& idx_i, idx_vector& idx_j)
+{
+  int nr = d1;
+  int nc = d2;
+
+  if (nr == 0 && nc == 0)
+    return;
+
+  if (idx_i.is_colon_equiv (nr, 1))
+    {
+      if (idx_j.is_colon_equiv (nc, 1))
+	resize (0, 0);
+      else
+	{
+	  int num_to_delete = idx_j.length (nc);
+
+	  if (num_to_delete != 0)
+	    {
+	      if (nr == 1 && num_to_delete == nc)
+		resize (0, 0);
+	      else
+		{
+		  int new_nc = nc - num_to_delete;
+		  if (new_nc > 0)
+		    {
+		      T *new_data = new T [nr * new_nc];
+
+		      int jj = 0;
+		      int idx = 0;
+		      for (int j = 0; j < nc; j++)
+			{
+			  if (j == idx_j.elem (idx))
+			    idx++;
+			  else
+			    {
+			      for (int i = 0; i < nr; i++)
+				new_data[nr*jj+i] = elem (i, j);
+			      jj++;
+			    }
+			}
+
+		      if (--rep->count <= 0)
+			delete rep;
+
+		      rep = new ArrayRep (new_data, nr * new_nc);
+
+		      d2 = new_nc;
+
+		      set_max_indices (2);
+		    }
+		  else
+		    (*current_liboctave_error_handler)
+		      ("A(idx) = []: index out of range");
+		}
+	    }
+	}
+    }
+  else if (idx_j.is_colon_equiv (nr, 1))
+    {
+      if (idx_i.is_colon_equiv (nc, 1))
+	resize (0, 0);
+      else
+	{
+	  int num_to_delete = idx_i.length (nr);
+
+	  if (num_to_delete != 0)
+	    {
+	      if (nc == 1 && num_to_delete == nr)
+		resize (0, 0);
+	      else 
+		{
+		  int new_nr = nr - num_to_delete;
+		  if (new_nr > 0)
+		    {
+		      T *new_data = new T [new_nr * nc];
+
+		      int ii = 0;
+		      int idx = 0;
+		      for (int i = 0; i < nr; i++)
+			{
+			  if (i == idx_i.elem (idx))
+			    idx++;
+			  else
+			    {
+			      for (int j = 0; j < nc; j++)
+				new_data[new_nr*j+ii] = elem (i, j);
+			      ii++;
+			    }
+			}
+
+		      if (--rep->count <= 0)
+			delete rep;
+
+		      rep = new ArrayRep (new_data, new_nr * nc);
+
+		      d1 = new_nr;
+
+		      set_max_indices (2);
+		    }
+		  else
+		    (*current_liboctave_error_handler)
+		      ("A(idx) = []: index out of range");
+		}
+	    }
+	}
+    }
+}
+
+template <class LT, class RT>
+int
+assign (Array2<LT>& lhs, const Array2<RT>& rhs)
+{
+  int retval = 1;
+
+  int n_idx = lhs.index_count ();
+
+  int lhs_nr = lhs.rows ();
+  int lhs_nc = lhs.cols ();
+
+  int rhs_nr = rhs.rows ();
+  int rhs_nc = rhs.cols ();
+
+  if (n_idx == 2)
+    {
+      idx_vector *tmp = lhs.get_idx ();
+
+      idx_vector idx_i = tmp[0];
+      idx_vector idx_j = tmp[1];
+
+      int n = idx_i.freeze (lhs_nr, "row", liboctave_pzo_flag,
+			    liboctave_rre_flag);
+
+      int m = idx_j.freeze (lhs_nc, "column", liboctave_pzo_flag,
+			    liboctave_rre_flag);
+
+      int idx_i_is_colon = idx_i.is_colon ();
+      int idx_j_is_colon = idx_j.is_colon ();
+
+      if (idx_i_is_colon)
+	n = rhs_nr;
+
+      if (idx_j_is_colon)
+	m = rhs_nc;
+
+      if (idx_i && idx_j)
+	{
+	  if (rhs_nr == 0 && rhs_nc == 0)
+	    {
+	      lhs.maybe_delete_elements (idx_i, idx_j);
+	    }
+	  else
+	    {
+	      if (liboctave_rre_flag)
+		{
+		  int max_row_idx = idx_i_is_colon ? rhs_nr : idx_i.max () + 1;
+		  int max_col_idx = idx_j_is_colon ? rhs_nc : idx_j.max () + 1;
+
+		  int new_nr = max_row_idx > lhs_nr ? max_row_idx : lhs_nr;
+		  int new_nc = max_col_idx > lhs_nc ? max_col_idx : lhs_nc;
+
+		  lhs.resize (new_nr, new_nc, 0.0);
+		}
+
+	      if (n == rhs_nr && m == rhs_nc)
+		{
+		  for (int j = 0; j < m; j++)
+		    {
+		      int jj = idx_j.elem (j);
+		      for (int i = 0; i < n; i++)
+			{
+			  int ii = idx_i.elem (i);
+			  lhs.elem (ii, jj) = rhs.elem (i, j);
+			}
+		    }
+		}
+	      else if (rhs_nr == 1 && rhs_nc == 1)
+		{
+		  RT scalar = rhs.elem (0, 0);
+
+		  for (int j = 0; j < m; j++)
+		    {
+		      int jj = idx_j.elem (j);
+		      for (int i = 0; i < n; i++)
+			{
+			  int ii = idx_i.elem (i);
+			  lhs.elem (ii, jj) = scalar;
+			}
+		    }
+		}
+	      else
+		{
+		  (*current_liboctave_error_handler)
+    ("A(I, J) = X: X must be a scalar or the number of elements in I must");
+		  (*current_liboctave_error_handler)
+    ("match the number of rows in X and the number of elements in J must");
+		  (*current_liboctave_error_handler)
+    ("match the number of columns in X");
+
+		  retval = 0;
+		}
+	    }
+	}
+      // idx_vector::freeze() printed an error message for us.
+    }
+  else if (n_idx == 1)
+    {
+      if (lhs_nr == 0 || lhs_nc == 0
+	  || (lhs_nr == 1 && lhs_nc == 1))
+	{
+	  idx_vector *tmp = lhs.get_idx ();
+
+	  idx_vector idx = tmp[0];
+
+	  int lhs_len = lhs.length ();
+
+	  int n = idx.freeze (lhs_len, 0, liboctave_pzo_flag,
+			      liboctave_rre_flag);
+
+	  if (idx)
+	    {
+	      if (rhs_nr == 0 && rhs_nc == 0)
+		{
+		  if (n != 0 && (lhs_nr != 0 || lhs_nc != 0))
+		    {
+		      idx_vector tmp (':');
+		      lhs.maybe_delete_elements (idx, tmp);
+		    }
+		}
+	      else
+		{
+		  if (assign ((Array<LT>&) lhs, (Array<RT>&) rhs))
+		    {
+		      int len = lhs.length ();
+
+		      if (len > 0)
+			{
+			  int idx_nr = idx.orig_rows ();
+			  int idx_nc = idx.orig_columns ();
+
+			  if (liboctave_dfi_flag
+			      || (idx_nr == 1 && idx_nc == 1))
+			    {
+			      if (liboctave_pcv_flag)
+				{
+				  lhs.d1 = lhs.length ();
+				  lhs.d2 = 1;
+				}
+			      else
+				{
+				  lhs.d1 = 1;
+				  lhs.d2 = lhs.length ();
+				}
+			    }
+			  else if (idx_nr == 1 && rhs_nr == 1)
+			    {
+			      lhs.d1 = 1;
+			      lhs.d2 = lhs.length ();
+			    }
+			  else if (idx_nc == 1 && rhs_nc == 1)
+			    {
+			      lhs.d1 = lhs.length ();
+			      lhs.d2 = 1;
+			    }
+			  else
+			    (*current_liboctave_error_handler)
+      ("A(I) = X: X must be a scalar or a matrix with the same size as I");
+			}
+		      else
+			{
+			  lhs.d1 = 0;
+			  lhs.d2 = 0;
+			}
+		    }
+		  else
+		    retval = 0;
+		}
+	    }
+	  // idx_vector::freeze() printed an error message for us.
+	}
+      else if (lhs_nr == 1)
+	{
+	  idx_vector *tmp = lhs.get_idx ();
+
+	  idx_vector idx = tmp[0];
+
+	  idx.freeze (lhs_nc, "vector", liboctave_pzo_flag,
+		      liboctave_rre_flag);
+
+	  if (idx)
+	    {
+	      if (rhs_nr == 0 && rhs_nc == 0)
+		{
+		  idx_vector tmp (':');
+		  lhs.maybe_delete_elements (tmp, idx);
+		}
+	      else
+		{
+		  if (assign ((Array<LT>&) lhs, (Array<RT>&) rhs))
+		    lhs.d2 = lhs.length ();
+		  else
+		    retval = 0;
+		}
+	    }
+	  // idx_vector::freeze() printed an error message for us.
+	}
+      else if (lhs_nc == 1)
+	{
+	  idx_vector *tmp = lhs.get_idx ();
+
+	  idx_vector idx = tmp[0];
+
+	  idx.freeze (lhs_nr, "vector", liboctave_pzo_flag,
+		      liboctave_rre_flag);
+
+	  if (idx)
+	    {
+	      if (rhs_nr == 0 && rhs_nc == 0)
+		{
+		  idx_vector tmp (':');
+		  lhs.maybe_delete_elements (idx, tmp);
+		}
+	      else
+		{
+		  if (assign ((Array<LT>&) lhs, (Array<RT>&) rhs))
+		    lhs.d1 = lhs.length ();
+		  else
+		    retval = 0;
+		}
+	    }
+	  // idx_vector::freeze() printed an error message for us.
+	}
+      else if (liboctave_dfi_flag)
+	{
+	  idx_vector *tmp = lhs.get_idx ();
+	  idx_vector idx = tmp[0];
+
+	  int len = idx.freeze (lhs_nr * lhs_nc, "matrix",
+				liboctave_pzo_flag);
+
+	  if (idx)
+	    {
+	      if (len == rhs_nr * rhs_nc)
+		{
+		  int k = 0;
+		  for (int j = 0; j < rhs_nc; j++)
+		    {
+		      for (int i = 0; i < rhs_nr; i++)
+			{
+			  int ii = idx.elem (k++);
+			  int fr = ii % lhs_nr;
+			  int fc = (ii - fr) / lhs_nr;
+			  lhs.elem (fr, fc) = rhs.elem (i, j);
+			}
+		    }
+		}
+	      else
+		{
+		  (*current_liboctave_error_handler)
+      ("A(I) = X: X must be a scalar or a matrix with the same size as I");
+
+		  retval = 0;
+		}
+	    }
+	  // idx_vector::freeze() printed an error message for us.
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("single index only valid for row or column vector");
+
+	  retval = 0;
+	}
+    }
+  else
+    {
+      (*current_liboctave_error_handler)
+	("invalid number of indices for matrix expression");
+
+      retval = 0;
+    }
+
+  lhs.clear_index ();
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/Array2.cc
@@ -0,0 +1,205 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cassert>
+
+#include <iostream.h>
+
+#include "Array2.h"
+
+#if defined (HEAVYWEIGHT_INDEXING)
+#include "idx-vector.h"
+#include "Array2-idx.h"
+#endif
+
+#include "lo-error.h"
+
+// Two dimensional array class.
+
+template <class T>
+T&
+Array2<T>::checkelem (int i, int j)
+{
+  if (i < 0 || j < 0 || i >= d1 || j >= d2)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      static T foo;
+      return foo;
+    }
+  return Array<T>::elem (d1*j+i);
+}
+
+template <class T>
+T
+Array2<T>::elem (int i, int j) const
+{
+  return Array<T>::elem (d1*j+i);
+}
+
+template <class T>
+T
+Array2<T>::checkelem (int i, int j) const
+{
+  if (i < 0 || j < 0 || i >= d1 || j >= d2)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      T foo;
+      static T *bar = &foo;
+      return foo;
+    }
+  return Array<T>::elem (d1*j+i);
+}
+
+template <class T>
+T
+Array2<T>::operator () (int i, int j) const
+{
+  if (i < 0 || j < 0 || i >= d1 || j >= d2)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      T foo;
+      static T *bar = &foo;
+      return foo;
+    }
+  return Array<T>::elem (d1*j+i);
+}
+
+template <class T>
+void
+Array2<T>::resize (int r, int c)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler) ("can't resize to negative dimension");
+      return;
+    }
+
+  if (r == dim1 () && c == dim2 ())
+    return;
+
+  ArrayRep *old_rep = rep;
+  const T *old_data = data ();
+
+  int old_d1 = dim1 ();
+  int old_d2 = dim2 ();
+  int old_len = length ();
+
+  rep = new ArrayRep (r*c);
+
+  d1 = r;
+  d2 = c;
+
+  if (old_data && old_len > 0)
+    {
+      int min_r = old_d1 < r ? old_d1 : r;
+      int min_c = old_d2 < c ? old_d2 : c;
+
+      for (int j = 0; j < min_c; j++)
+	for (int i = 0; i < min_r; i++)
+	  xelem (i, j) = old_data[old_d1*j+i];
+    }
+
+  if (--old_rep->count <= 0)
+    delete old_rep;
+}
+
+template <class T>
+void
+Array2<T>::resize (int r, int c, const T& val)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler) ("can't resize to negative dimension");
+      return;
+    }
+
+  if (r == dim1 () && c == dim2 ())
+    return;
+
+  ArrayRep *old_rep = rep;
+  const T *old_data = data ();
+  int old_d1 = dim1 ();
+  int old_d2 = dim2 ();
+  int old_len = length ();
+
+  rep = new ArrayRep (r*c);
+
+  d1 = r;
+  d2 = c;
+
+  int min_r = old_d1 < r ? old_d1 : r;
+  int min_c = old_d2 < c ? old_d2 : c;
+
+  if (old_data && old_len > 0)
+    {
+      for (int j = 0; j < min_c; j++)
+	for (int i = 0; i < min_r; i++)
+	  xelem (i, j) = old_data[old_d1*j+i];
+    }
+
+  for (int j = 0; j < min_c; j++)
+    for (int i = min_r; i < r; i++)
+      xelem (i, j) = val;
+
+  for (int j = min_c; j < c; j++)
+    for (int i = 0; i < r; i++)
+      xelem (i, j) = val;
+
+  if (--old_rep->count <= 0)
+    delete old_rep;
+}
+
+template <class T>
+Array2<T>&
+Array2<T>::insert (const Array2<T>& a, int r, int c)
+{
+  int a_rows = a.rows ();
+  int a_cols = a.cols ();
+
+  if (r < 0 || r + a_rows > rows () || c < 0 || c + a_cols > 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;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/Array2.h
@@ -0,0 +1,155 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_Array2_h)
+#define octave_Array2_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#define HEAVYWEIGHT_INDEXING 1
+
+#include <cassert>
+#include <cstdlib>
+
+#include "Array.h"
+#include "lo-error.h"
+
+class idx_vector;
+
+// Two dimensional array class.
+
+template <class T>
+class Array2 : public Array<T>
+{
+protected:
+
+  Array2 (T *d, int n, int m) : Array<T> (d, n*m)
+    {
+      d1 = n;
+      d2 = m;
+      set_max_indices (2);
+    }
+
+public:
+
+  // These really need to be protected (and they will be in the
+  // future, so don't depend on them being here!), but they can't be
+  // until template friends work correctly in g++.
+
+  int d1;
+  int d2;
+
+  Array2 (void) : Array<T> ()
+    {
+      d1 = 0;
+      d2 = 0;
+      set_max_indices (2);
+    }
+
+  Array2 (int n, int m) : Array<T> (n*m)
+    {
+      d1 = n;
+      d2 = m;
+      set_max_indices (2);
+    }
+
+  Array2 (int n, int m, const T& val) : Array<T> (n*m, val)
+    {
+      d1 = n;
+      d2 = m;
+      set_max_indices (2);
+    }
+
+  Array2 (const Array2<T>& a) : Array<T> (a)
+    {
+      d1 = a.d1;
+      d2 = a.d2;
+      set_max_indices (2);
+    }
+
+  Array2 (const Array<T>& a, int n, int m) : Array<T> (a)
+    {
+      d1 = n;
+      d2 = m;
+      set_max_indices (2);
+    }
+
+  ~Array2 (void) { }
+
+  Array2<T>& operator = (const Array2<T>& a)
+    {
+      if (this != &a && rep != a.rep)
+	{
+	  Array<T>::operator = (a);
+	  d1 = a.d1;
+	  d2 = a.d2;
+	}
+
+      return *this;
+    }
+
+  int dim1 (void) const { return d1; }
+  int dim2 (void) const { return d2; }
+
+  int rows (void) const { return d1; }
+  int cols (void) const { return d2; }
+  int columns (void) const { return d2; }
+
+  T& elem (int i, int j) { return Array<T>::elem (d1*j+i); }
+  T& checkelem (int i, int j);
+  T& operator () (int i, int j) { return checkelem (i, j); }
+
+  T elem (int i, int j) const;
+  T checkelem (int i, int j) const;
+  T operator () (int i, int j) const;
+
+  // No checking.
+
+  T& xelem (int i, int j) { return Array<T>::xelem (d1*j+i); }
+  T xelem (int i, int j) const { return Array<T>::xelem (d1*j+i); }
+
+  void resize (int n, int m);
+  void resize (int n, int m, const T& val);
+
+  Array2<T>& insert (const Array2<T>& a, int r, int c);
+
+#ifdef HEAVYWEIGHT_INDEXING
+  void maybe_delete_elements (idx_vector& i, idx_vector& j);
+
+  Array2<T> value (void);
+#endif
+};
+
+template <class LT, class RT>
+int assign (Array2<LT>& lhs, const Array2<RT>& rhs);
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/Array3-idx.h
@@ -0,0 +1,57 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#include "Array-flags.h"
+#include "idx-vector.h"
+#include "lo-error.h"
+
+template <class T>
+void
+Array3<T>::maybe_delete_elements (idx_vector&, idx_vector&, idx_vector&)
+{
+  assert (0);
+}
+
+template <class T>
+Array3<T>
+Array3<T>::value (void)
+{
+  Array3<T> retval;
+  assert (0);
+  return retval;
+}
+
+template <class LT, class RT>
+int
+assign (Array3<LT>&, const Array3<RT>&)
+{
+  assert (0);
+  return 0;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/Array3.cc
@@ -0,0 +1,114 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cassert>
+
+#include <iostream.h>
+
+#include "Array3.h"
+
+#if defined (HEAVYWEIGHT_INDEXING)
+#include "idx-vector.h"
+#include "Array3-idx.h"
+#endif
+
+#include "lo-error.h"
+
+// Three dimensional array class.
+
+template <class T>
+T&
+Array3<T>::checkelem (int i, int j, int k)
+{
+  if (i < 0 || j < 0 || k < 0 || i >= d1 || j >= d2 || k >= d3)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      static T foo;
+      return foo;
+    }
+  return Array2<T>::elem (i, d1*k+j);
+}
+
+template <class T>
+T
+Array3<T>::elem (int i, int j, int k) const
+{
+  return Array2<T>::elem (i, d2*k+j);
+}
+
+template <class T>
+T
+Array3<T>::checkelem (int i, int j, int k) const
+{
+  if (i < 0 || j < 0 || k < 0 || i >= d1 || j >= d2 || k >= d3)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      T foo;
+      static T *bar = &foo;
+      return foo;
+    }
+  return Array2<T>::elem (i, d1*k+j);
+}
+
+template <class T>
+T
+Array3<T>::operator () (int i, int j, int k) const
+{
+  if (i < 0 || j < 0 || k < 0 || i >= d1 || j >= d2 || k >= d3)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      T foo;
+      static T *bar = &foo;
+      return foo;
+    }
+  return Array2<T>::elem (i, d2*k+j);
+}
+
+template <class T>
+void
+Array3<T>::resize (int n, int m, int k)
+{
+  assert (0); // XXX FIXME XXX
+}
+
+template <class T>
+void
+Array3<T>::resize (int n, int m, int k, const T& val)
+{
+  assert (0); // XXX FIXME XXX
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/Array3.h
@@ -0,0 +1,136 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_Array3_h)
+#define octave_Array3_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#define HEAVYWEIGHT_INDEXING 1
+
+#include <cassert>
+#include <cstdlib>
+
+#include "Array2.h"
+#include "lo-error.h"
+
+class idx_vector;
+
+// Three dimensional array class.
+
+template <class T>
+class Array3 : public Array2<T>
+{
+protected:
+
+  int d3;
+
+  Array3 (T *d, int n, int m, int k) : Array2<T> (d, n, m*k)
+    {
+      d2 = m;
+      d3 = k;
+      set_max_indices (3);
+    }
+
+public:
+
+  Array3 (void) : Array2<T> ()
+    {
+      d2 = 0;
+      d3 = 0;
+      set_max_indices (3);
+    }
+
+  Array3 (int n, int m, int k) : Array2<T> (n, m*k)
+    {
+      d2 = m;
+      d3 = k;
+      set_max_indices (3);
+    }
+
+  Array3 (int n, int m, int k, const T& val) : Array2<T> (n, m*k, val)
+    {
+      d2 = m;
+      d3 = k;
+      set_max_indices (3);
+    }
+
+  Array3 (const Array3<T>& a) : Array2<T> (a)
+    {
+      d2 = a.d2;
+      d3 = a.d3;
+      set_max_indices (3);
+    }
+
+  ~Array3 (void) { }
+
+  Array3<T>& operator = (const Array3<T>& a)
+    {
+      if (this != &a && rep != a.rep)
+	{
+	  Array<T>::operator = (a);
+	  d1 = a.d1;
+	  d2 = a.d2;
+	  d3 = a.d3;
+	}
+
+      return *this;
+    }
+
+  int dim3 (void) const { return d3; }
+
+  T& elem (int i, int j, int k) { return Array2<T>::elem (i, d2*k+j); }
+  T& checkelem (int i, int j, int k);
+  T& operator () (int i, int j, int k) { return checkelem (i, j, k); }
+
+  // No checking.
+
+  T& xelem (int i, int j, int k) { return Array2<T>::xelem (i, d2*k+j); }
+
+  T elem (int i, int j, int k) const;
+  T checkelem (int i, int j, int k) const;
+  T operator () (int i, int j, int k) const;
+
+  void resize (int n, int m, int k);
+  void resize (int n, int m, int k, const T& val);
+
+#ifdef HEAVYWEIGHT_INDEXING
+  void maybe_delete_elements (idx_vector& i, idx_vector& j, idx_vector& k);
+
+  Array3<T> value (void);
+#endif
+};
+
+template <class LT, class RT>
+int assign (Array3<LT>& lhs, const Array3<RT>& rhs);
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/DiagArray2.cc
@@ -0,0 +1,208 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cassert>
+
+#include <iostream.h>
+
+#include "DiagArray2.h"
+
+#include "lo-error.h"
+
+// A two-dimensional array with diagonal elements only.
+
+#if 0
+template <class T>
+T&
+DiagArray2<T>::elem (int r, int c)
+{
+  static T foo (0);
+  return (r == c) ? Array<T>::xelem (r) : foo;
+}
+
+template <class T>
+T&
+DiagArray2<T>::checkelem (int r, int c)
+{
+  static T foo (0);
+  if (r < 0 || c < 0 || r >= nr || c >= nc)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      return foo;
+    }
+  return (r == c) ? Array<T>::xelem (r) : foo;
+}
+
+template <class T>
+T&
+DiagArray2<T>::operator () (int r, int c)
+{
+  static T foo (0);
+  if (r < 0 || c < 0 || r >= nr || c >= nc)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      return foo;
+    }
+  return (r == c) ? Array<T>::xelem (r) : foo;
+}
+#endif
+
+template <class T>
+T
+DiagArray2<T>::elem (int r, int c) const
+{
+  return (r == c) ? Array<T>::xelem (r) : T (0);
+}
+
+template <class T>
+T
+DiagArray2<T>::checkelem (int r, int c) const
+{
+  if (r < 0 || c < 0 || r >= nr || c >= nc)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      T foo;
+      static T *bar = &foo;
+      return foo;
+    }
+  return (r == c) ? Array<T>::xelem (r) : T (0);
+}
+
+template <class T>
+T
+DiagArray2<T>::operator () (int r, int c) const
+{
+  if (r < 0 || c < 0 || r >= nr || c >= nc)
+    {
+      (*current_liboctave_error_handler) ("range error");
+      T foo;
+      static T *bar = &foo;
+      return foo;
+    }
+  return (r == c) ? Array<T>::xelem (r) : T (0);
+}
+
+template <class T>
+T&
+DiagArray2<T>::xelem (int r, int c)
+{
+  static T foo (0);
+  return (r == c) ? Array<T>::xelem (r) : foo;
+}
+
+template <class T>
+T
+DiagArray2<T>::xelem (int r, int c) const
+{
+  return (r == c) ? Array<T>::xelem (r) : T (0);
+}
+
+template <class T>
+void
+DiagArray2<T>::resize (int r, int c)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler) ("can't resize to negative dimensions");
+      return;
+    }
+
+  if (r == dim1 () && c == dim2 ())
+    return;
+
+  ArrayRep *old_rep = rep;
+  const T *old_data = data ();
+  int old_len = length ();
+
+  int new_len = r < c ? r : c;
+
+  rep = new ArrayRep (new_len);
+
+  nr = r;
+  nc = c;
+
+  if (old_data && old_len > 0)
+    {
+      int min_len = old_len < new_len ? old_len : new_len;
+
+      for (int i = 0; i < min_len; i++)
+	xelem (i, i) = old_data[i];
+    }
+
+  if (--old_rep->count <= 0)
+    delete old_rep;
+}
+
+template <class T>
+void
+DiagArray2<T>::resize (int r, int c, const T& val)
+{
+  if (r < 0 || c < 0)
+    {
+      (*current_liboctave_error_handler) ("can't resize to negative dimensions");
+      return;
+    }
+
+  if (r == dim1 () && c == dim2 ())
+    return;
+
+  ArrayRep *old_rep = rep;
+  const T *old_data = data ();
+  int old_len = length ();
+
+  int new_len = r < c ? r : c;
+
+  rep = new ArrayRep (new_len);
+
+  nr = r;
+  nc = c;
+
+  int min_len = old_len < new_len ? old_len : new_len;
+
+  if (old_data && old_len > 0)
+    {
+      for (int i = 0; i < min_len; i++)
+	xelem (i, i) = old_data[i];
+    }
+
+  for (int i = min_len; i < new_len; i++)
+    xelem (i, i) = val;
+
+  if (--old_rep->count <= 0)
+    delete old_rep;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/DiagArray2.h
@@ -0,0 +1,250 @@
+// Template array classes                              -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_DiagArray2_h)
+#define octave_DiagArray2_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#define HEAVYWEIGHT_INDEXING 1
+
+#include <cassert>
+#include <cstdlib>
+
+#include "Array2.h"
+#include "lo-error.h"
+
+class idx_vector;
+
+// A two-dimensional array with diagonal elements only.
+//
+// Idea and example code for Proxy class and functions from:
+//
+// From: kanze@us-es.sel.de (James Kanze)
+// Subject: Re: How to overload [] to do READ/WRITE differently ?
+// Message-ID: <KANZE.93Nov29151407@slsvhdt.us-es.sel.de>
+// Sender: news@us-es.sel.de
+// Date: 29 Nov 1993 14:14:07 GMT
+// --
+// James Kanze                             email: kanze@us-es.sel.de
+// GABI Software, Sarl., 8 rue du Faisan, F-67000 Strasbourg, France
+
+template <class T>
+class DiagArray2 : public Array<T>
+{
+private:
+
+  T get (int i) { return Array<T>::xelem (i); }
+
+  void set (const T& val, int i) { Array<T>::xelem (i) = val; }
+
+  class Proxy
+  {
+  public:
+
+    Proxy (DiagArray2<T> *ref, int r, int c)
+      : i (r), j (c), object (ref) { } 
+
+    const Proxy& operator = (const T& val) const
+      {
+	if (i == j)
+	  {
+	    if (object)
+	      object->set (val, i);
+	  }
+	else
+	  (*current_liboctave_error_handler)
+	    ("invalid assignment to off-diagonal in diagonal array");
+
+	return *this;
+      }
+
+    operator T () const
+      {
+	if (object && i == j)
+	  return object->get (i);
+	else
+	  {
+	    static T foo (0);
+	    return foo;
+	  }
+      }
+
+  private:
+
+    // XXX FIXME XXX -- this is declared private to keep the user from
+    // taking the address of a Proxy.  Maybe it should be implemented
+    // by means of a companion function in the DiagArray2 class.
+
+    T *operator& () const { assert (0); return (T *) 0; }
+
+    int i;
+    int j;
+
+    DiagArray2<T> *object;
+
+  };
+
+friend class Proxy;
+
+protected:
+
+  int nr;
+  int nc;
+
+  DiagArray2 (T *d, int r, int c) : Array<T> (d, r < c ? r : c)
+    {
+      nr = r;
+      nc = c;
+      set_max_indices (2);
+    }
+
+public:
+
+  DiagArray2 (void) : Array<T> ()
+    {
+      nr = 0;
+      nc = 0;
+      set_max_indices (2);
+    }
+
+  DiagArray2 (int r, int c) : Array<T> (r < c ? r : c)
+    {
+      nr = r;
+      nc = c;
+      set_max_indices (2);
+    }
+
+  DiagArray2 (int r, int c, const T& val) : Array<T> (r < c ? r : c, val)
+    {
+      nr = r;
+      nc = c;
+      set_max_indices (2);
+    }
+
+  DiagArray2 (const Array<T>& a) : Array<T> (a)
+    {
+      nr = nc = a.length ();
+      set_max_indices (2);
+    }
+
+  DiagArray2 (const DiagArray2<T>& a) : Array<T> (a)
+    {
+      nr = a.nr;
+      nc = a.nc;
+      set_max_indices (2);
+    }
+
+  ~DiagArray2 (void) { }
+
+  DiagArray2<T>& operator = (const DiagArray2<T>& a)
+    {
+      if (this != &a)
+	{
+	  Array<T>::operator = (a);
+	  nr = a.nr;
+	  nc = a.nc;
+	}
+
+      return *this;
+    }
+
+#if 0
+  operator Array2<T> () const
+    {
+      Array2<T> retval (nr, nc,  T (0));
+
+      int len = nr < nc ? nr : nc;
+
+      for (int i = 0; i < len; i++)
+	retval.xelem (i, i) = xelem (i, i);
+
+      return retval;
+    }
+#endif
+
+  int dim1 (void) const { return nr; }
+  int dim2 (void) const { return nc; }
+
+  int rows (void) const { return nr; }
+  int cols (void) const { return nc; }
+  int columns (void) const { return nc; }
+
+#if 1
+  Proxy elem (int r, int c)
+    {
+      return Proxy (this, r, c);
+    }
+
+  Proxy checkelem (int r, int c)
+    {
+      if (r < 0 || c < 0 || r >= nr || c >= nc)
+	{
+	  (*current_liboctave_error_handler) ("range error");
+	  return Proxy (0, r, c);
+	}
+      else
+	return Proxy (this, r, c);
+    }
+
+  Proxy operator () (int r, int c)
+    {
+      if (r < 0 || c < 0 || r >= nr || c >= nc)
+	{
+	  (*current_liboctave_error_handler) ("range error");
+	  return Proxy (0, r, c);
+	}
+      else
+	return Proxy (this, r, c);
+  }
+#else
+  T& elem (int r, int c);
+  T& checkelem (int r, int c);
+  T& operator () (int r, int c);
+#endif
+
+  T elem (int r, int c) const;
+  T checkelem (int r, int c) const;
+  T operator () (int r, int c) const;
+
+  // No checking.
+
+  T& xelem (int r, int c);
+  T xelem (int r, int c) const;
+
+  void resize (int n, int m);
+  void resize (int n, int m, const T& val);
+
+  void maybe_delete_elements (idx_vector& i, idx_vector& j);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/MArray-defs.h
@@ -0,0 +1,67 @@
+// Nothing like a little CPP abuse to brighten everyone's day.  Would
+// have been nice to do this with template functions but as of 2.5.x,
+// g++ seems to fail to resolve them properly.
+
+#define DO_VS_OP(OP) \
+  int l = a.length (); \
+  T *result = 0; \
+  if (l > 0) \
+    { \
+      result = new T [l]; \
+      const T *x = a.data (); \
+      for (int i = 0; i < l; i++) \
+	result[i] = x[i] OP s; \
+    }
+
+#define DO_SV_OP(OP) \
+  int l = a.length (); \
+  T *result = 0; \
+  if (l > 0) \
+    { \
+      result = new T [l]; \
+      const T *x = a.data (); \
+      for (int i = 0; i < l; i++) \
+	result[i] = s OP x[i]; \
+    }
+
+#define DO_VV_OP(OP) \
+  T *result = 0; \
+  if (l > 0) \
+    { \
+      result = new T [l]; \
+      const T *x = a.data (); \
+      const T *y = b.data (); \
+      for (int i = 0; i < l; i++) \
+	result[i] = x[i] OP y[i]; \
+    }
+
+#define NEG_V \
+  int l = a.length (); \
+  T *result = 0; \
+  if (l > 0) \
+    { \
+      result = new T [l]; \
+      const T *x = a.data (); \
+      for (int i = 0; i < l; i++) \
+	result[i] = -x[i]; \
+    }
+
+#define DO_VS_OP2(OP) \
+  int l = a.length (); \
+  if (l > 0) \
+    { \
+      T *tmp = a.fortran_vec (); \
+      for (int i = 0; i < l; i++) \
+	tmp[i] OP s; \
+    }
+
+#define DO_VV_OP2(OP) \
+  do \
+    { \
+      T *a_tmp = a.fortran_vec (); \
+      const T *b_tmp = b.data (); \
+      for (int i = 0; i < l; i++) \
+	a_tmp[i] += b_tmp[i]; \
+    } \
+  while (0)
+
new file mode 100644
--- /dev/null
+++ b/liboctave/MArray2.cc
@@ -0,0 +1,177 @@
+// MArray.cc                                             -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "MArray2.h"
+#include "lo-error.h"
+
+#include "MArray-defs.h"
+
+// Two dimensional array with math ops.
+
+// Element by element MArray2 by scalar ops.
+
+template <class T>
+MArray2<T>&
+operator += (MArray2<T>& a, const T& s)
+{
+  DO_VS_OP2 (+=)
+  return a;
+}
+
+template <class T>
+MArray2<T>&
+operator -= (MArray2<T>& a, const T& s)
+{
+  DO_VS_OP2 (-=)
+  return a;
+}
+
+// Element by element MArray2 by MArray2 ops.
+
+template <class T>
+MArray2<T>&
+operator += (MArray2<T>& a, const MArray2<T>& b)
+{
+  int r = a.rows ();
+  int c = a.cols ();
+  if (r != b.rows () || c != b.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant += array operation attempted");
+    }
+  else
+    {
+      if (r > 0 && c > 0)
+	{
+	  int l = a.length ();
+	  DO_VV_OP2 (+=);
+	}
+    }
+  return a;
+}
+
+template <class T>
+MArray2<T>&
+operator -= (MArray2<T>& a, const MArray2<T>& b)
+{
+  int r = a.rows ();
+  int c = a.cols ();
+  if (r != b.rows () || c != b.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant -= array operation attempted");
+    }
+  else
+    {
+      if (r > 0 && c > 0)
+	{
+	  int l = a.length ();
+	  DO_VV_OP2 (-=);
+	}
+    }
+  return a;
+}
+
+// Element by element MArray2 by scalar ops.
+
+#define MARRAY_A2S_OP(OP) \
+  template <class T> \
+  MArray2<T> \
+  operator OP (const MArray2<T>& a, const T& s) \
+  { \
+    DO_VS_OP (OP); \
+    return MArray2<T> (result, a.rows (), a.cols ()); \
+  }
+
+MARRAY_A2S_OP (+)
+MARRAY_A2S_OP (-)
+MARRAY_A2S_OP (*)
+MARRAY_A2S_OP (/)
+
+// Element by element scalar by MArray2 ops.
+
+#define MARRAY_SA2_OP(OP) \
+  template <class T> \
+  MArray2<T> \
+  operator OP (const T& s, const MArray2<T>& a) \
+  { \
+    DO_SV_OP (OP); \
+    return MArray2<T> (result, a.rows (), a.cols ()); \
+  }
+
+MARRAY_SA2_OP (+)
+MARRAY_SA2_OP (-)
+MARRAY_SA2_OP (*)
+MARRAY_SA2_OP (/)
+
+// Element by element MArray2 by MArray2 ops.
+
+#define MARRAY_A2A2_OP(FCN, OP, OP_STR) \
+  template <class T> \
+  MArray2<T> \
+  FCN (const MArray2<T>& a, const MArray2<T>& b) \
+  { \
+    int r = a.rows (); \
+    int c = a.cols (); \
+    if (r != b.rows () || c != b.cols ()) \
+      { \
+	(*current_liboctave_error_handler) \
+	  ("nonconformant array " OP_STR " attempted"); \
+	return MArray2<T> (); \
+      } \
+    if (r == 0 || c == 0) \
+      return MArray2<T> (); \
+    int l = a.length (); \
+    DO_VV_OP (OP); \
+    return MArray2<T> (result, r, c); \
+  }
+
+MARRAY_A2A2_OP (operator +, +, "addition")
+MARRAY_A2A2_OP (operator -, -, "subtraction")
+MARRAY_A2A2_OP (product,    *, "product")
+MARRAY_A2A2_OP (quotient,   /, "quotient")
+
+// Unary MArray2 ops.
+
+template <class T>
+MArray2<T>
+operator - (const MArray2<T>& a)
+{
+  NEG_V;
+  return MArray2<T> (result, a.rows (), a.cols ());
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/MArray2.h
@@ -0,0 +1,115 @@
+// Template array classes with like-type math ops          -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#if !defined (octave_MArray2_h)
+#define octave_MArray2_h 1
+
+#include "Array2.h"
+
+// Two dimensional array with math ops.
+
+template <class T>
+class MArray2 : public Array2<T>
+{
+protected:
+
+  MArray2 (T *d, int n, int m) : Array2<T> (d, n, m) { }
+
+public:
+
+  MArray2 (void) : Array2<T> () { }
+  MArray2 (int n, int m) : Array2<T> (n, m) { }
+  MArray2 (int n, int m, const T& val) : Array2<T> (n, m, val) { }
+  MArray2 (const Array2<T>& a) : Array2<T> (a) { }
+  MArray2 (const MArray2<T>& a) : Array2<T> (a) { }
+
+  ~MArray2 (void) { }
+
+  MArray2<T>& operator = (const MArray2<T>& a)
+    {
+      Array2<T>::operator = (a);
+      return *this;
+    }
+
+  // element by element MArray2 by scalar ops
+
+  friend MArray2<T>& operator += (MArray2<T>& a, const T& s);
+  friend MArray2<T>& operator -= (MArray2<T>& a, const T& s);
+
+  // element by element MArray2 by MArray2 ops
+
+  friend MArray2<T>& operator += (MArray2<T>& a, const MArray2<T>& b);
+  friend MArray2<T>& operator -= (MArray2<T>& a, const MArray2<T>& b);
+
+  // element by element MArray2 by scalar ops
+
+  friend MArray2<T> operator + (const MArray2<T>& a, const T& s);
+  friend MArray2<T> operator - (const MArray2<T>& a, const T& s);
+  friend MArray2<T> operator * (const MArray2<T>& a, const T& s);
+  friend MArray2<T> operator / (const MArray2<T>& a, const T& s);
+
+  // element by element scalar by MArray2 ops
+
+  friend MArray2<T> operator + (const T& s, const MArray2<T>& a);
+  friend MArray2<T> operator - (const T& s, const MArray2<T>& a);
+  friend MArray2<T> operator * (const T& s, const MArray2<T>& a);
+  friend MArray2<T> operator / (const T& s, const MArray2<T>& a);
+
+  // element by element MArray2 by MArray2 ops
+
+  friend MArray2<T> operator + (const MArray2<T>& a, const MArray2<T>& b);
+  friend MArray2<T> operator - (const MArray2<T>& a, const MArray2<T>& b);
+
+  friend MArray2<T> product (const MArray2<T>& a, const MArray2<T>& b);
+  friend MArray2<T> quotient (const MArray2<T>& a, const MArray2<T>& b);
+
+  friend MArray2<T> operator - (const MArray2<T>& a);
+};
+
+#define INSTANTIATE_MARRAY2_FRIENDS(T) \
+  template MArray2<T> operator + (const MArray2<T>& a, const T& s); \
+  template MArray2<T> operator - (const MArray2<T>& a, const T& s); \
+  template MArray2<T> operator * (const MArray2<T>& a, const T& s); \
+  template MArray2<T> operator / (const MArray2<T>& a, const T& s); \
+  template MArray2<T> operator + (const T& s, const MArray2<T>& a); \
+  template MArray2<T> operator - (const T& s, const MArray2<T>& a); \
+  template MArray2<T> operator * (const T& s, const MArray2<T>& a); \
+  template MArray2<T> operator / (const T& s, const MArray2<T>& a); \
+  template MArray2<T> operator + (const MArray2<T>& a, const MArray2<T>& b); \
+  template MArray2<T> operator - (const MArray2<T>& a, const MArray2<T>& b); \
+  template MArray2<T> product (const MArray2<T>& a, const MArray2<T>& b); \
+  template MArray2<T> quotient (const MArray2<T>& a, const MArray2<T>& b); \
+  template MArray2<T> operator - (const MArray2<T>& a);
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/MDiagArray2.cc
@@ -0,0 +1,146 @@
+// MArray.cc                                             -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "MDiagArray2.h"
+#include "lo-error.h"
+
+#include "MArray-defs.h"
+
+// Two dimensional diagonal array with math ops.
+
+// Element by element MDiagArray2 by MDiagArray2 ops.
+
+template <class T>
+MDiagArray2<T>&
+operator += (MDiagArray2<T>& a, const MDiagArray2<T>& b)
+{
+  int r = a.rows ();
+  int c = a.cols ();
+  if (r != b.rows () || c != b.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant array operator += attempted");
+      return MDiagArray2<T> ();
+    }
+  else
+    {
+      int l = a.length ();
+      DO_VV_OP2 (+=);
+    }
+  return a;
+}
+
+template <class T>
+MDiagArray2<T>&
+operator -= (MDiagArray2<T>& a, const MDiagArray2<T>& b)
+{
+  int r = a.rows ();
+  int c = a.cols ();
+  if (r != b.rows () || c != b.cols ())
+    {
+      (*current_liboctave_error_handler)
+	("nonconformant array operator -= attempted");
+      return MDiagArray2<T> ();
+    }
+  else
+    {
+      int l = a.length ();
+      DO_VV_OP2 (-=);
+    }
+  return a;
+}
+
+// Element by element MDiagArray2 by scalar ops.
+
+#define MARRAY_DAS_OP(OP) \
+  template <class T> \
+  MDiagArray2<T> \
+  operator OP (const MDiagArray2<T>& a, const T& s) \
+  { \
+    DO_VS_OP (OP); \
+    return MDiagArray2<T> (result, a.rows (), a.cols ()); \
+  }
+
+MARRAY_DAS_OP (*)
+MARRAY_DAS_OP (/)
+
+// Element by element scalar by MDiagArray2 ops.
+
+template <class T>
+MDiagArray2<T>
+operator * (const T& s, const MDiagArray2<T>& a)
+{
+  DO_SV_OP (*);
+  return MDiagArray2<T> (result, a.rows (), a.cols ());
+}
+
+// Element by element MDiagArray2 by MDiagArray2 ops.
+
+#define MARRAY_DADA_OP(FCN, OP, OP_STR) \
+  template <class T> \
+  MDiagArray2<T> \
+  FCN (const MDiagArray2<T>& a, const MDiagArray2<T>& b) \
+  { \
+    int r = a.rows (); \
+    int c = a.cols (); \
+    if (r != b.rows () || c != b.cols ()) \
+      { \
+	(*current_liboctave_error_handler) \
+	  ("nonconformant diagonal array " OP_STR " attempted"); \
+	return MDiagArray2<T> (); \
+      } \
+    if (c == 0 || r == 0) \
+      return MDiagArray2<T> (); \
+    int l = a.length (); \
+    DO_VV_OP (OP); \
+    return MDiagArray2<T> (result, r, c); \
+  }
+
+MARRAY_DADA_OP (operator +, +, "addition")
+MARRAY_DADA_OP (operator -, -, "subtraction")
+MARRAY_DADA_OP (product,    *, "product")
+
+// Unary MDiagArray2 ops.
+
+template <class T>
+MDiagArray2<T>
+operator - (const MDiagArray2<T>& a)
+{
+  NEG_V;
+  return MDiagArray2<T> (result, a.rows (), a.cols ());
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/MDiagArray2.h
@@ -0,0 +1,120 @@
+// Template array classes with like-type math ops          -*- C++ -*-
+/*
+
+Copyright (C) 1996 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, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#if !defined (octave_MDiagArray2_h)
+#define octave_MDiagArray2_h 1
+
+#include "DiagArray2.h"
+#include "MArray2.h"
+
+// Two dimensional diagonal array with math ops.
+
+template <class T>
+class MDiagArray2 : public DiagArray2<T>
+{
+protected:
+
+  MDiagArray2 (T *d, int r, int c) : DiagArray2<T> (d, r, c) { }
+
+public:
+  
+  MDiagArray2 (void) : DiagArray2<T> () { }
+  MDiagArray2 (int r, int c) : DiagArray2<T> (r, c) { }
+  MDiagArray2 (int r, int c, const T& val) : DiagArray2<T> (r, c, val) { }
+  MDiagArray2 (const Array<T>& a) : DiagArray2<T> (a) { }
+  MDiagArray2 (const MArray<T>& a) : DiagArray2<T> (a) { }
+  MDiagArray2 (const DiagArray2<T>& a) : DiagArray2<T> (a) { }
+  MDiagArray2 (const MDiagArray2<T>& a) : DiagArray2<T> (a) { }
+
+  ~MDiagArray2 (void) { }
+
+  MDiagArray2<T>& operator = (const MDiagArray2<T>& a)
+    {
+      DiagArray2<T>::operator = (a);
+      return *this;
+    }
+
+  operator MArray2<T> () const
+    {
+      MArray2<T> retval (nr, nc,  T (0));
+
+      int len = nr < nc ? nr : nc;
+
+      for (int i = 0; i < len; i++)
+	retval.xelem (i, i) = xelem (i, i);
+
+      return retval;
+    }
+
+  // element by element MDiagArray2 by MDiagArray2 ops
+
+  friend MDiagArray2<T>&
+  operator += (MDiagArray2<T>& a, const MDiagArray2<T>& b);
+
+  friend MDiagArray2<T>&
+  operator -= (MDiagArray2<T>& a, const MDiagArray2<T>& b);
+
+  // element by element MDiagArray2 by scalar ops
+
+  friend MDiagArray2<T> operator * (const MDiagArray2<T>& a, const T& s);
+  friend MDiagArray2<T> operator / (const MDiagArray2<T>& a, const T& s);
+
+  // element by element scalar by MDiagArray2 ops
+
+  friend MDiagArray2<T> operator * (const T& s, const MDiagArray2<T>& a);
+
+  // element by element MDiagArray2 by MDiagArray2 ops
+
+  friend MDiagArray2<T>
+  operator + (const MDiagArray2<T>& a, const MDiagArray2<T>& b); 
+
+  friend MDiagArray2<T>
+  operator - (const MDiagArray2<T>& a, const MDiagArray2<T>& b);
+
+  friend MDiagArray2<T>
+  product (const MDiagArray2<T>& a, const MDiagArray2<T>& b);
+
+  friend MDiagArray2<T> operator - (const MDiagArray2<T>& a);
+};
+
+#define INSTANTIATE_MDIAGARRAY_FRIENDS(T) \
+  template MDiagArray2<T> operator * (const MDiagArray2<T>& a, const T& s); \
+  template MDiagArray2<T> operator / (const MDiagArray2<T>& a, const T& s); \
+  template MDiagArray2<T> operator * (const T& s, const MDiagArray2<T>& a); \
+  template MDiagArray2<T> operator + (const MDiagArray2<T>& a, const MDiagArray2<T>& b); \
+  template MDiagArray2<T> operator - (const MDiagArray2<T>& a, const MDiagArray2<T>& b); \
+  template MDiagArray2<T> product (const MDiagArray2<T>& a, const MDiagArray2<T>& b); \
+  template MDiagArray2<T> operator - (const MDiagArray2<T>& a);
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; page-delimiter: "^/\\*" ***
+;;; End: ***
+*/