# HG changeset patch
# User jwe
# Date 1080926814 0
# Node ID 9f7ef92b50b05e5bf1da7d3a37450f3331eb541e
# Parent  7b4e761009645d403d5cd6f3cab794cb037bfffd
[project @ 2004-04-02 17:26:53 by jwe]

diff --git a/liboctave/CNDArray.cc b/liboctave/CNDArray.cc
--- a/liboctave/CNDArray.cc
+++ b/liboctave/CNDArray.cc
@@ -683,6 +683,192 @@
   return ::cat_ra(*this, ra_arg, dim, iidx, move);
 }
 
+static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
+
+ComplexNDArray
+ComplexNDArray::max (int dim) const
+{
+  ArrayN<int> dummy_idx;
+  return max (dummy_idx, dim);
+}
+
+ComplexNDArray
+ComplexNDArray::max (ArrayN<int>& idx_arg, int dim) const
+{
+  dim_vector dv = dims ();
+  dim_vector dr = dims ();
+
+  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+    return ComplexNDArray ();
+  
+  dr(dim) = 1;
+
+  ComplexNDArray result (dr);
+  idx_arg.resize (dr);
+
+  int x_stride = 1;
+  int x_len = dv(dim);
+  for (int i = 0; i < dim; i++)
+    x_stride *= dv(i);
+
+  for (int i = 0; i < dr.numel (); i++)
+    {
+      int x_offset;
+      if (x_stride == 1)
+	x_offset = i * x_len;
+      else
+	{
+	  int x_offset2 = 0;
+	  x_offset = i;
+	  while (x_offset >= x_stride)
+	    {
+	      x_offset -= x_stride;
+	      x_offset2++;
+	    }
+	  x_offset += x_offset2 * x_stride * x_len;
+	}
+
+      int idx_j;
+
+      Complex tmp_max;
+
+      double abs_max = octave_NaN;
+
+      for (idx_j = 0; idx_j < x_len; idx_j++)
+	{
+	  tmp_max = elem (idx_j * x_stride + x_offset);
+	  
+	  if (! octave_is_NaN_or_NA (tmp_max))
+	    {
+	      abs_max = ::abs(tmp_max);
+	      break;
+	    }
+	}
+
+      for (int j = idx_j+1; j < x_len; j++)
+	{
+	  Complex tmp = elem (j * x_stride + x_offset);
+
+	  if (octave_is_NaN_or_NA (tmp))
+	    continue;
+
+	  double abs_tmp = ::abs (tmp);
+
+	  if (abs_tmp > abs_max)
+	    {
+	      idx_j = j;
+	      tmp_max = tmp;
+	      abs_max = abs_tmp;
+	    }
+	}
+
+      if (octave_is_NaN_or_NA (tmp_max))
+	{
+	  result.elem (i) = Complex_NaN_result;
+	  idx_arg.elem (i) = 0;
+	}
+      else
+	{
+	  result.elem (i) = tmp_max;
+	  idx_arg.elem (i) = idx_j;
+	}
+    }
+
+  return result;
+}
+
+ComplexNDArray
+ComplexNDArray::min (int dim) const
+{
+  ArrayN<int> dummy_idx;
+  return min (dummy_idx, dim);
+}
+
+ComplexNDArray
+ComplexNDArray::min (ArrayN<int>& idx_arg, int dim) const
+{
+  dim_vector dv = dims ();
+  dim_vector dr = dims ();
+
+  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+    return ComplexNDArray ();
+  
+  dr(dim) = 1;
+
+  ComplexNDArray result (dr);
+  idx_arg.resize (dr);
+
+  int x_stride = 1;
+  int x_len = dv(dim);
+  for (int i = 0; i < dim; i++)
+    x_stride *= dv(i);
+
+  for (int i = 0; i < dr.numel (); i++)
+    {
+      int x_offset;
+      if (x_stride == 1)
+	x_offset = i * x_len;
+      else
+	{
+	  int x_offset2 = 0;
+	  x_offset = i;
+	  while (x_offset >= x_stride)
+	    {
+	      x_offset -= x_stride;
+	      x_offset2++;
+	    }
+	  x_offset += x_offset2 * x_stride * x_len;
+	}
+
+      int idx_j;
+
+      Complex tmp_min;
+
+      double abs_min = octave_NaN;
+
+      for (idx_j = 0; idx_j < x_len; idx_j++)
+	{
+	  tmp_min = elem (idx_j * x_stride + x_offset);
+	  
+	  if (! octave_is_NaN_or_NA (tmp_min))
+	    {
+	      abs_min = ::abs(tmp_min);
+	      break;
+	    }
+	}
+
+      for (int j = idx_j+1; j < x_len; j++)
+	{
+	  Complex tmp = elem (j * x_stride + x_offset);
+
+	  if (octave_is_NaN_or_NA (tmp))
+	    continue;
+
+	  double abs_tmp = ::abs (tmp);
+
+	  if (abs_tmp < abs_min)
+	    {
+	      idx_j = j;
+	      tmp_min = tmp;
+	      abs_min = abs_tmp;
+	    }
+	}
+
+      if (octave_is_NaN_or_NA (tmp_min))
+	{
+	  result.elem (i) = Complex_NaN_result;
+	  idx_arg.elem (i) = 0;
+	}
+      else
+	{
+	  result.elem (i) = tmp_min;
+	  idx_arg.elem (i) = idx_j;
+	}
+    }
+
+  return result;
+}
+
 NDArray
 ComplexNDArray::abs (void) const
 {
@@ -836,6 +1022,141 @@
   return is;
 }
 
+// XXX FIXME XXX -- it would be nice to share code among the min/max
+// functions below.
+
+#define EMPTY_RETURN_CHECK(T) \
+  if (nel == 0)	\
+    return T (dv);
+
+ComplexNDArray
+min (const Complex& c, const ComplexNDArray& m)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (ComplexNDArray);
+
+  ComplexNDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmin (c, m (i));
+    }
+
+  return result;
+}
+
+ComplexNDArray
+min (const ComplexNDArray& m, const Complex& c)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (ComplexNDArray);
+
+  ComplexNDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmin (c, m (i));
+    }
+
+  return result;
+}
+
+ComplexNDArray
+min (const ComplexNDArray& a, const ComplexNDArray& b)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  if (dv != b.dims ())
+    {
+      (*current_liboctave_error_handler)
+	("two-arg min expecting args of same size");
+      return ComplexNDArray ();
+    }
+
+  EMPTY_RETURN_CHECK (ComplexNDArray);
+
+  ComplexNDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmin (a (i), b (i));
+    }
+
+  return result;
+}
+
+ComplexNDArray
+max (const Complex& c, const ComplexNDArray& m)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (ComplexNDArray);
+
+  ComplexNDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmax (c, m (i));
+    }
+
+  return result;
+}
+
+ComplexNDArray
+max (const ComplexNDArray& m, const Complex& c)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (ComplexNDArray);
+
+  ComplexNDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmax (c, m (i));
+    }
+
+  return result;
+}
+
+ComplexNDArray
+max (const ComplexNDArray& a, const ComplexNDArray& b)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  if (dv != b.dims ())
+    {
+      (*current_liboctave_error_handler)
+	("two-arg max expecting args of same size");
+      return ComplexNDArray ();
+    }
+
+  EMPTY_RETURN_CHECK (ComplexNDArray);
+
+  ComplexNDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmax (a (i), b (i));
+    }
+
+  return result;
+}
+
 NDS_CMP_OPS(ComplexNDArray, real, Complex, real)
 NDS_BOOL_OPS(ComplexNDArray, Complex, 0.0)
 
diff --git a/liboctave/CNDArray.h b/liboctave/CNDArray.h
--- a/liboctave/CNDArray.h
+++ b/liboctave/CNDArray.h
@@ -89,6 +89,11 @@
   ComplexNDArray sumsq (int dim = -1) const;
   int cat (const ComplexNDArray& ra_arg, int dim, int iidx, int move);
 
+  ComplexNDArray max (int dim = 0) const;
+  ComplexNDArray max (ArrayN<int>& index, int dim = 0) const;
+  ComplexNDArray min (int dim = 0) const;
+  ComplexNDArray min (ArrayN<int>& index, int dim = 0) const;
+
   ComplexNDArray& insert (const NDArray& a, int r, int c);
   ComplexNDArray& insert (const ComplexNDArray& a, int r, int c);
   
@@ -130,6 +135,14 @@
     : MArrayN<Complex> (d, dv) { }
 };
 
+extern ComplexNDArray min (const Complex& c, const ComplexNDArray& m);
+extern ComplexNDArray min (const ComplexNDArray& m, const Complex& c);
+extern ComplexNDArray min (const ComplexNDArray& a, const ComplexNDArray& b);
+
+extern ComplexNDArray max (const Complex& c, const ComplexNDArray& m);
+extern ComplexNDArray max (const ComplexNDArray& m, const Complex& c);
+extern ComplexNDArray max (const ComplexNDArray& a, const ComplexNDArray& b);
+
 NDS_CMP_OP_DECLS (ComplexNDArray, Complex)
 NDS_BOOL_OP_DECLS (ComplexNDArray, Complex)
 
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,21 @@
+2004-04-02  David Bateman  <dbateman@free.fr>
+
+	* lo-specfun.cc (besselj, bessely, besseli, besselk, besselh1, 
+	besselh2, airy, biry, betainc, gammainc, do_bessel):
+	New NDArray versions.
+	(SN_BESSEL, NS_BESSEL, NN_BESSEL): New macros.
+	* lo-specfun.h (besselj, bessely, besseli, besselk, besselh1, 
+	besselh2, airy, biry, betainc, gammainc): Provide decls.
+	
+	* dNDArray.cc (NDArray::min, NDArray::max, min, max):
+	New functions.
+	* dNDArray.h (NDArray::min, NDArray::max, min, max): Provide decls.
+
+	* CNDArray.cc (ComplexNDArray::min, ComplexNDArray::max, min, max):
+	New functions.
+	* CNDArray.h (ComplexNDArray::min, ComplexNDArray::max, min, max): 
+	Provide decls.
+
 2004-03-17  David Hoover  <jazzdaq@yahoo.com>
 
 	* DASPK.cc (DASPK::do_integrate): Always add n*n elements to the
diff --git a/liboctave/dNDArray.cc b/liboctave/dNDArray.cc
--- a/liboctave/dNDArray.cc
+++ b/liboctave/dNDArray.cc
@@ -656,6 +656,156 @@
   MX_ND_REAL_OP_REDUCTION (+= elem (iter_idx), 0);
 }
 
+NDArray
+NDArray::max (int dim) const
+{
+  ArrayN<int> dummy_idx;
+  return max (dummy_idx, dim);
+}
+
+NDArray
+NDArray::max (ArrayN<int>& idx_arg, int dim) const
+{
+  dim_vector dv = dims ();
+  dim_vector dr = dims ();
+
+  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+    return NDArray ();
+  
+  dr(dim) = 1;
+
+  NDArray result (dr);
+  idx_arg.resize (dr);
+
+  int x_stride = 1;
+  int x_len = dv(dim);
+  for (int i = 0; i < dim; i++)
+    x_stride *= dv(i);
+
+  for (int i = 0; i < dr.numel (); i++)
+    {
+      int x_offset;
+      if (x_stride == 1)
+	x_offset = i * x_len;
+      else
+	{
+	  int x_offset2 = 0;
+	  x_offset = i;
+	  while (x_offset >= x_stride)
+	    {
+	      x_offset -= x_stride;
+	      x_offset2++;
+	    }
+	  x_offset += x_offset2 * x_stride * x_len;
+	}
+
+      int idx_j;
+
+      double tmp_max = octave_NaN;
+
+      for (idx_j = 0; idx_j < x_len; idx_j++)
+	{
+	  tmp_max = elem (idx_j * x_stride + x_offset);
+	  
+	  if (! octave_is_NaN_or_NA (tmp_max))
+	    break;
+	}
+
+      for (int j = idx_j+1; j < x_len; j++)
+	{
+	  double tmp = elem (j * x_stride + x_offset);
+
+	  if (octave_is_NaN_or_NA (tmp))
+	    continue;
+	  else if (tmp > tmp_max)
+	    {
+	      idx_j = j;
+	      tmp_max = tmp;
+	    }
+	}
+
+      result.elem (i) = tmp_max;
+      idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j;
+    }
+
+  return result;
+}
+
+NDArray
+NDArray::min (int dim) const
+{
+  ArrayN<int> dummy_idx;
+  return min (dummy_idx, dim);
+}
+
+NDArray
+NDArray::min (ArrayN<int>& idx_arg, int dim) const
+{
+  dim_vector dv = dims ();
+  dim_vector dr = dims ();
+
+  if (dv.numel () == 0 || dim > dv.length () || dim < 0)
+    return NDArray ();
+  
+  dr(dim) = 1;
+
+  NDArray result (dr);
+  idx_arg.resize (dr);
+
+  int x_stride = 1;
+  int x_len = dv(dim);
+  for (int i = 0; i < dim; i++)
+    x_stride *= dv(i);
+
+  for (int i = 0; i < dr.numel (); i++)
+    {
+      int x_offset;
+      if (x_stride == 1)
+	x_offset = i * x_len;
+      else
+	{
+	  int x_offset2 = 0;
+	  x_offset = i;
+	  while (x_offset >= x_stride)
+	    {
+	      x_offset -= x_stride;
+	      x_offset2++;
+	    }
+	  x_offset += x_offset2 * x_stride * x_len;
+	}
+
+      int idx_j;
+
+      double tmp_min = octave_NaN;
+
+      for (idx_j = 0; idx_j < x_len; idx_j++)
+	{
+	  tmp_min = elem (idx_j * x_stride + x_offset);
+	  
+	  if (! octave_is_NaN_or_NA (tmp_min))
+	    break;
+	}
+
+      for (int j = idx_j+1; j < x_len; j++)
+	{
+	  double tmp = elem (j * x_stride + x_offset);
+
+	  if (octave_is_NaN_or_NA (tmp))
+	    continue;
+	  else if (tmp < tmp_min)
+	    {
+	      idx_j = j;
+	      tmp_min = tmp;
+	    }
+	}
+
+      result.elem (i) = tmp_min;
+      idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j;
+    }
+
+  return result;
+}
+
 int
 NDArray::cat (const NDArray& ra_arg, int dim, int iidx, int move)
 {
@@ -776,6 +926,141 @@
   return is;
 }
 
+// XXX FIXME XXX -- it would be nice to share code among the min/max
+// functions below.
+
+#define EMPTY_RETURN_CHECK(T) \
+  if (nel == 0)	\
+    return T (dv);
+
+NDArray
+min (double d, const NDArray& m)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (NDArray);
+
+  NDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmin (d, m (i));
+    }
+
+  return result;
+}
+
+NDArray
+min (const NDArray& m, double d)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (NDArray);
+
+  NDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmin (d, m (i));
+    }
+
+  return result;
+}
+
+NDArray
+min (const NDArray& a, const NDArray& b)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  if (dv != b.dims ())
+    {
+      (*current_liboctave_error_handler)
+	("two-arg min expecting args of same size");
+      return NDArray ();
+    }
+
+  EMPTY_RETURN_CHECK (NDArray);
+
+  NDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmin (a (i), b (i));
+    }
+
+  return result;
+}
+
+NDArray
+max (double d, const NDArray& m)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (NDArray);
+
+  NDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmax (d, m (i));
+    }
+
+  return result;
+}
+
+NDArray
+max (const NDArray& m, double d)
+{
+  dim_vector dv = m.dims ();
+  int nel = dv.numel ();
+
+  EMPTY_RETURN_CHECK (NDArray);
+
+  NDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmax (d, m (i));
+    }
+
+  return result;
+}
+
+NDArray
+max (const NDArray& a, const NDArray& b)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  if (dv != b.dims ())
+    {
+      (*current_liboctave_error_handler)
+	("two-arg max expecting args of same size");
+      return NDArray ();
+    }
+
+  EMPTY_RETURN_CHECK (NDArray);
+
+  NDArray result (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      result (i) = xmax (a (i), b (i));
+    }
+
+  return result;
+}
+
 NDS_CMP_OPS(NDArray, , double, )
 NDS_BOOL_OPS(NDArray, double, 0.0)
 
diff --git a/liboctave/dNDArray.h b/liboctave/dNDArray.h
--- a/liboctave/dNDArray.h
+++ b/liboctave/dNDArray.h
@@ -87,7 +87,12 @@
   NDArray sum (int dim = -1) const;  
   NDArray sumsq (int dim = -1) const;
   int cat (const NDArray& ra_arg, int dim, int iidx, int move);
-     
+
+  NDArray max (int dim = 0) const;
+  NDArray max (ArrayN<int>& index, int dim = 0) const;
+  NDArray min (int dim = 0) const;
+  NDArray min (ArrayN<int>& index, int dim = 0) const;
+  
   NDArray abs (void) const;
 
   ComplexNDArray fourier (int dim = 1) const;
@@ -125,6 +130,14 @@
   NDArray (double *d, const dim_vector& dv) : MArrayN<double> (d, dv) { }
 };
 
+extern NDArray min (double d, const NDArray& m);
+extern NDArray min (const NDArray& m, double d);
+extern NDArray min (const NDArray& a, const NDArray& b);
+
+extern NDArray max (double d, const NDArray& m);
+extern NDArray max (const NDArray& m, double d);
+extern NDArray max (const NDArray& a, const NDArray& b);
+
 NDS_CMP_OP_DECLS (NDArray, double)
 NDS_BOOL_OP_DECLS (NDArray, double)
 
diff --git a/liboctave/lo-specfun.cc b/liboctave/lo-specfun.cc
--- a/liboctave/lo-specfun.cc
+++ b/liboctave/lo-specfun.cc
@@ -29,6 +29,8 @@
 #include "CMatrix.h"
 #include "dRowVector.h"
 #include "dMatrix.h"
+#include "dNDArray.h"
+#include "CNDArray.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
@@ -607,6 +609,62 @@
   return retval;
 }
 
+static inline ComplexNDArray
+do_bessel (fptr f, const char *, double alpha, const ComplexNDArray& x,
+	   bool scaled, ArrayN<int>& ierr)
+{
+  dim_vector dv = x.dims ();
+  int nel = dv.numel ();
+  ComplexNDArray retval (dv);
+
+  ierr.resize (dv);
+
+  for (int i = 0; i < nel; i++)
+      retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
+
+  return retval;
+}
+
+static inline ComplexNDArray
+do_bessel (fptr f, const char *, const NDArray& alpha, const Complex& x,
+	   bool scaled, ArrayN<int>& ierr)
+{
+  dim_vector dv = alpha.dims ();
+  int nel = dv.numel ();
+  ComplexNDArray retval (dv);
+
+  ierr.resize (dv);
+
+  for (int i = 0; i < nel; i++)
+    retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
+
+  return retval;
+}
+
+static inline ComplexNDArray
+do_bessel (fptr f, const char *fn, const NDArray& alpha,
+	   const ComplexNDArray& x, bool scaled, ArrayN<int>& ierr)
+{
+  dim_vector dv = x.dims ();
+  ComplexNDArray retval;
+
+  if (dv == alpha.dims ())
+    {
+      int nel = dv.numel ();
+
+      retval.resize (dv);
+      ierr.resize (dv);
+
+      for (int i = 0; i < nel; i++)
+	retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("%s: the sizes of alpha and x must conform", fn);
+
+  return retval;
+}
+
 static inline ComplexMatrix
 do_bessel (fptr f, const char *, const RowVector& alpha,
 	   const ComplexColumnVector& x, bool scaled, Array2<int>& ierr)
@@ -656,6 +714,30 @@
     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
   }
 
+#define SN_BESSEL(name, fcn) \
+  ComplexNDArray \
+  name (double alpha, const ComplexNDArray& x, bool scaled, \
+	ArrayN<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define NS_BESSEL(name, fcn) \
+  ComplexNDArray \
+  name (const NDArray& alpha, const Complex& x, bool scaled, \
+	ArrayN<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
+#define NN_BESSEL(name, fcn) \
+  ComplexNDArray \
+  name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \
+	ArrayN<int>& ierr) \
+  { \
+    return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
+  }
+
 #define RC_BESSEL(name, fcn) \
   ComplexMatrix \
   name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
@@ -669,6 +751,9 @@
   SM_BESSEL (name, fcn) \
   MS_BESSEL (name, fcn) \
   MM_BESSEL (name, fcn) \
+  SN_BESSEL (name, fcn) \
+  NS_BESSEL (name, fcn) \
+  NN_BESSEL (name, fcn) \
   RC_BESSEL (name, fcn)
 
 ALL_BESSEL (besselj, zbesj)
@@ -778,6 +863,36 @@
   return retval;
 }
 
+ComplexNDArray
+airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
+{
+  dim_vector dv = z.dims ();
+  int nel = dv.numel ();
+  ComplexNDArray retval (dv);
+
+  ierr.resize (dv);
+
+  for (int i = 0; i < nel; i++)
+    retval (i) = airy (z(i), deriv, scaled, ierr(i));
+
+  return retval;
+}
+
+ComplexNDArray
+biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
+{
+  dim_vector dv = z.dims ();
+  int nel = dv.numel ();
+  ComplexNDArray retval (dv);
+
+  ierr.resize (dv);
+
+  for (int i = 0; i < nel; i++)
+    retval (i) = biry (z(i), deriv, scaled, ierr(i));
+
+  return retval;
+}
+
 static void
 gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3,
 			     int c3)
@@ -787,6 +902,21 @@
      r1, c1, r2, c2, r3, c3);
 }
 
+static dim_vector null_dims (0);
+
+static void
+gripe_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2,
+			     const dim_vector& d3)
+{
+  std::string d1_str = d1.str ();
+  std::string d2_str = d2.str ();
+  std::string d3_str = d3.str ();
+
+  (*current_liboctave_error_handler)
+  ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
+   d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
+}
+
 double
 betainc (double x, double a, double b)
 {
@@ -850,6 +980,56 @@
   return retval;
 }
 
+NDArray
+betainc (double x, double a, const NDArray& b)
+{
+  dim_vector dv = b.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval (dv);
+
+  for (int i = 0; i < nel; i++)
+    retval (i) = betainc (x, a, b(i));
+
+  return retval;
+}
+
+NDArray
+betainc (double x, const NDArray& a, double b)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval (dv);
+
+  for (int i = 0; i < nel; i++)
+    retval (i) = betainc (x, a(i), b);
+
+  return retval;
+}
+
+NDArray
+betainc (double x, const NDArray& a, const NDArray& b)
+{
+  NDArray retval;
+  dim_vector dv = a.dims ();
+
+  if (dv == b.dims ())
+    {
+      int nel = dv.numel ();
+
+      retval.resize (dv);
+
+      for (int i = 0; i < nel; i++)
+	retval (i) = betainc (x, a(i), b(i));
+    }
+  else
+    gripe_betainc_nonconformant (dim_vector (0), dv, b.dims ());
+  
+  return retval;
+}
+
+
 Matrix
 betainc (const Matrix& x, double a, double b)
 {
@@ -943,6 +1123,83 @@
   return retval;
 }
 
+NDArray
+betainc (const NDArray& x, double a, double b)
+{
+  dim_vector dv = x.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval (dv);
+
+  for (int i = 0; i < nel; i++)
+    retval (i) = betainc (x(i), a, b);
+
+  return retval;
+}
+
+NDArray
+betainc (const NDArray& x, double a, const NDArray& b)
+{
+  NDArray retval;
+  dim_vector dv = x.dims ();
+
+  if (dv == b.dims ())
+    {
+      int nel = dv.numel ();
+
+      retval.resize (dv);
+
+      for (int i = 0; i < nel; i++)
+	retval (i) = betainc (x(i), a, b(i));
+    }
+  else
+    gripe_betainc_nonconformant (dv, dim_vector (0), b.dims ());
+  
+  return retval;
+}
+
+NDArray
+betainc (const NDArray& x, const NDArray& a, double b)
+{
+  NDArray retval;
+  dim_vector dv = x.dims ();
+
+  if (dv == a.dims ())
+    {
+      int nel = dv.numel ();
+
+      retval.resize (dv);
+
+      for (int i = 0; i < nel; i++)
+	retval (i) = betainc (x(i), a(i), b);
+    }
+  else
+    gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0));
+  
+  return retval;
+}
+
+NDArray
+betainc (const NDArray& x, const NDArray& a, const NDArray& b)
+{
+  NDArray retval;
+  dim_vector dv = x.dims ();
+
+  if (dv == a.dims () && dv == b.dims ())
+    {
+      int nel = dv.numel ();
+
+      retval.resize (dv);
+
+      for (int i = 0; i < nel; i++)
+	retval (i) = betainc (x(i), a(i), b(i));
+    }
+  else
+    gripe_betainc_nonconformant (dv, a.dims (), b.dims ());
+
+  return retval;
+}
+
 // XXX FIXME XXX -- there is still room for improvement here...
 
 double
@@ -1058,6 +1315,98 @@
   return retval;
 }
 
+NDArray
+gammainc (double x, const NDArray& a)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval;
+  NDArray result (dv);
+
+  bool err;
+
+  for (int i = 0; i < nel; i++)
+    {
+      result (i) = gammainc (x, a(i), err);
+
+      if (err)
+	goto done;
+    }
+
+  retval = result;
+
+ done:
+
+  return retval;
+}
+
+NDArray
+gammainc (const NDArray& x, double a)
+{
+  dim_vector dv = x.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval;
+  NDArray result (dv);
+
+  bool err;
+
+  for (int i = 0; i < nel; i++)
+    {
+      result (i) = gammainc (x(i), a, err);
+
+      if (err)
+	goto done;
+    }
+
+  retval = result;
+
+ done:
+
+  return retval;
+}
+
+NDArray
+gammainc (const NDArray& x, const NDArray& a)
+{
+  dim_vector dv = x.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval;
+  NDArray result;
+
+  if (dv == a.dims ())
+    {
+      result.resize (dv);
+
+      bool err;
+
+      for (int i = 0; i < nel; i++)
+	{
+	  result (i) = gammainc (x(i), a(i), err);
+	  
+	  if (err)
+	    goto done;
+	}
+
+      retval = result;
+    }
+  else
+    {
+      std::string x_str = dv.str ();
+      std::string a_str = a.dims ().str ();
+
+      (*current_liboctave_error_handler)
+	("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
+	 x_str.c_str (), a_str. c_str ());
+    }
+
+ done:
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/lo-specfun.h b/liboctave/lo-specfun.h
--- a/liboctave/lo-specfun.h
+++ b/liboctave/lo-specfun.h
@@ -24,10 +24,13 @@
 #define octave_liboctave_specfun_h 1
 
 #include "oct-cmplx.h"
+#include "ArrayN.h"
 
 template <class T> class Array2;
 class Matrix;
 class ComplexMatrix;
+class NDArray;
+class ComplexNDArray;
 class RowVector;
 class ComplexColumnVector;
 class Range;
@@ -145,6 +148,78 @@
 besselh2 (const Matrix& alpha, const ComplexMatrix& x, bool scaled,
 	  Array2<int>& ierr);
 
+extern ComplexNDArray
+besselj (double alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+bessely (double alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besseli (double alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselk (double alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselh1 (double alpha, const ComplexNDArray& x, bool scaled,
+	  ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselh2 (double alpha, const ComplexNDArray& x, bool scaled,
+	  ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselj (const NDArray& alpha, const Complex& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+bessely (const NDArray& alpha, const Complex& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besseli (const NDArray& alpha, const Complex& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselk (const NDArray& alpha, const Complex& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselh1 (const NDArray& alpha, const Complex& x, bool scaled,
+	  ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselh2 (const NDArray& alpha, const Complex& x, bool scaled,
+	  ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselj (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+bessely (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besseli (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselk (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+	 ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselh1 (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+	  ArrayN<int>& ierr);
+
+extern ComplexNDArray
+besselh2 (const NDArray& alpha, const ComplexNDArray& x, bool scaled,
+	  ArrayN<int>& ierr);
+
 extern ComplexMatrix
 besselj (const RowVector& alpha, const ComplexColumnVector& x, bool scaled,
 	 Array2<int>& ierr);
@@ -178,21 +253,40 @@
 extern ComplexMatrix
 biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr);
 
+extern ComplexNDArray
+airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr);
+
+extern ComplexNDArray
+biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr);
+
 extern double betainc (double x, double a, double b);
 extern Matrix betainc (double x, double a, const Matrix& b);
 extern Matrix betainc (double x, const Matrix& a, double b);
 extern Matrix betainc (double x, const Matrix& a, const Matrix& b);
 
+extern NDArray betainc (double x, double a, const NDArray& b);
+extern NDArray betainc (double x, const NDArray& a, double b);
+extern NDArray betainc (double x, const NDArray& a, const NDArray& b);
+
 extern Matrix betainc (const Matrix& x, double a, double b);
 extern Matrix betainc (const Matrix& x, double a, const Matrix& b);
 extern Matrix betainc (const Matrix& x, const Matrix& a, double b);
 extern Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b);
 
+extern NDArray betainc (const NDArray& x, double a, double b);
+extern NDArray betainc (const NDArray& x, double a, const NDArray& b);
+extern NDArray betainc (const NDArray& x, const NDArray& a, double b);
+extern NDArray betainc (const NDArray& x, const NDArray& a, const NDArray& b);
+
 extern double gammainc (double x, double a, bool& err);
 extern Matrix gammainc (double x, const Matrix& a);
 extern Matrix gammainc (const Matrix& x, double a);
 extern Matrix gammainc (const Matrix& x, const Matrix& a);
 
+extern NDArray gammainc (double x, const NDArray& a);
+extern NDArray gammainc (const NDArray& x, double a);
+extern NDArray gammainc (const NDArray& x, const NDArray& a);
+
 inline double gammainc (double x, double a)
 {
   bool err;
diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,25 @@
+2004-04-02  David Bateman  <dbateman@free.fr>
+
+	* statistics/base/std.m: Allow optional args for type and dim.
+	* statistics/base/center.m statistics/base/meansq.m
+	statistics/base/moment.m statistics/base/range.m: Update for NDArrays.
+	* signal/fftshift.m: Fix dimensioning error.
+	
+	* statistics/base/std.m: Use repmat not ones(nr,1)*mean to allow
+	NDArrays.
+	
+	* general/mod.m, general/mod.m: Allow NDArrays with one scalar arg.
+
+	* signal/fftshift.m: Update for NDArrays, allow optional dim arg.
+	
+	* specfun/erfinv.m, general/repmat.m: Update for NDArrays.
+	
+	* control/base/bode.m, control/base/lqg.m, control/system/ss2sys.m,
+	control/system/cellidx.m, control/system/dmr2d.m control/system/ss.m,
+	control/system/sysprune.m: Doc update for usage of cell arrays.
+
+	* control/system/sysidx.m: Use cellidx and not listidx.
+
 2004-03-12  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* plot/__pltopt1__.m: Always add title clause to plot command with
diff --git a/scripts/control/base/bode.m b/scripts/control/base/bode.m
--- a/scripts/control/base/bode.m
+++ b/scripts/control/base/bode.m
@@ -67,7 +67,7 @@
 ##
 ## @strong{Example}
 ## @example
-## bode(sys,[],"y_3",list("u_1","u_4");
+## bode(sys,[],"y_3", @{"u_1","u_4"@});
 ## @end example
 ## @end table
 ## @strong{Outputs}
diff --git a/scripts/control/base/lqg.m b/scripts/control/base/lqg.m
--- a/scripts/control/base/lqg.m
+++ b/scripts/control/base/lqg.m
@@ -40,7 +40,7 @@
 ## @itemx  r
 ## state, control weighting respectively.  Control ARE is
 ## @item  in_idx
-## names or indices of controlled inputs (see @code{sysidx}, @code{listidx})
+## names or indices of controlled inputs (see @code{sysidx}, @code{cellidx})
 ##
 ## default: last dim(R) inputs are assumed to be controlled inputs, all
 ## others are assumed to be noise inputs.
diff --git a/scripts/control/system/cellidx.m b/scripts/control/system/cellidx.m
--- a/scripts/control/system/cellidx.m
+++ b/scripts/control/system/cellidx.m
@@ -17,7 +17,7 @@
 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{idxvec}, @var{errmsg}] =} listidx (@var{listvar}, @var{strlist})
+## @deftypefn {Function File} {[@var{idxvec}, @var{errmsg}] =} cellidx (@var{listvar}, @var{strlist})
 ## Return indices of string entries in @var{listvar} that match strings
 ## in @var{strlist}.
 ##
@@ -29,7 +29,7 @@
 ##
 ## If @var{strlist} contains a string not in @var{listvar}, then
 ## an error message is returned in @var{errmsg}.  If only one output
-## argument is requested, then @var{listidx} prints @var{errmsg} to the
+## argument is requested, then @var{cellidx} prints @var{errmsg} to the
 ## screen and exits with an error.
 ## @end deftypefn
 
diff --git a/scripts/control/system/dmr2d.m b/scripts/control/system/dmr2d.m
--- a/scripts/control/system/dmr2d.m
+++ b/scripts/control/system/dmr2d.m
@@ -29,7 +29,7 @@
 ## @code{dmr2d} exits with an error if @var{sys} is not discrete
 ## @item   idx
 ## indices or names of states with sampling time 
-## @code{sysgettsam(@var{sys})} (may be empty); see @code{listidx}
+## @code{sysgettsam(@var{sys})} (may be empty); see @code{cellidx}
 ## @item   sprefix
 ## list of string prefixes of states with sampling time
 ## @code{sysgettsam(@var{sys})} (may be empty)
diff --git a/scripts/control/system/ss.m b/scripts/control/system/ss.m
--- a/scripts/control/system/ss.m
+++ b/scripts/control/system/ss.m
@@ -46,18 +46,18 @@
 ## see below for system partitioning
 ##
 ## @item  stname
-## list of strings of state signal names
+## cell array of strings of state signal names
 ##
 ## default (@var{stname}=[] on input): @code{x_n} for continuous states,
 ##                     @code{xd_n} for discrete states
 ##
 ## @item inname
-## list of strings of input signal names
+## cell array of strings of input signal names
 ##
 ## default (@var{inname} = [] on input): @code{u_n}
 ##
 ## @item outname
-## list of strings of input signal names
+## cell array of strings of input signal names
 ##
 ## default (@var{outname} = [] on input): @code{y_n}
 ##
@@ -140,7 +140,7 @@
 ## octave:1> a = [1 2 3; 4 5 6; 7 8 10];
 ## octave:2> b = [0 0 ; 0 1 ; 1 0];
 ## octave:3> c = eye(3);
-## octave:4> sys = ss(a,b,c,[],0,3,0,list("volts","amps","joules"));
+## octave:4> sys = ss(a,b,c,[],0,3,0,@{"volts","amps","joules"@});
 ## octave:5> sysout(sys);
 ## Input(s)
 ##         1: u_1
diff --git a/scripts/control/system/ss2sys.m b/scripts/control/system/ss2sys.m
--- a/scripts/control/system/ss2sys.m
+++ b/scripts/control/system/ss2sys.m
@@ -46,18 +46,18 @@
 ## see below for system partitioning
 ##
 ## @item  stname
-## list of strings of state signal names
+## cell array of strings of state signal names
 ##
 ## default (@var{stname}=[] on input): @code{x_n} for continuous states,
 ##                     @code{xd_n} for discrete states
 ##
 ## @item inname
-## list of strings of input signal names
+## cell array of strings of input signal names
 ##
 ## default (@var{inname} = [] on input): @code{u_n}
 ##
 ## @item outname
-## list of strings of input signal names
+## cell array of strings of input signal names
 ##
 ## default (@var{outname} = [] on input): @code{y_n}
 ##
@@ -140,7 +140,7 @@
 ## octave:1> a = [1 2 3; 4 5 6; 7 8 10];
 ## octave:2> b = [0 0 ; 0 1 ; 1 0];
 ## octave:3> c = eye(3);
-## octave:4> sys = ss(a,b,c,[],0,3,0,list("volts","amps","joules"));
+## octave:4> sys = ss(a,b,c,[],0,3,0,@{"volts","amps","joules"@});
 ## octave:5> sysout(sys);
 ## Input(s)
 ##         1: u_1
diff --git a/scripts/control/system/sysidx.m b/scripts/control/system/sysidx.m
--- a/scripts/control/system/sysidx.m
+++ b/scripts/control/system/sysidx.m
@@ -38,13 +38,13 @@
   endif
 
   ## extract correct set of signal names values
-  [idxvec, msg] = listidx (list ("in", "out", "st", "yd"), sigtype);
+  [idxvec, msg] = cellidx ({"in", "out", "st", "yd"}, sigtype);
   if (msg)
     error ("invalid sigtype = %s", sigtype);
   endif
 
   syssiglist = sysgetsignals (sys, sigtype);
-  [idxvec, msg] = listidx (syssiglist, signamelist);
+  [idxvec, msg] = cellidx (syssiglist, signamelist);
   if (length (msg))
     error ("sysidx (sigtype = %s): %s", sigtype,
 	   strrep (msg, "strlist", "signamelist"));
diff --git a/scripts/control/system/sysprune.m b/scripts/control/system/sysprune.m
--- a/scripts/control/system/sysprune.m
+++ b/scripts/control/system/sysprune.m
@@ -33,7 +33,7 @@
 ##
 ## @example
 ## retsys = sysprune(Asys,[1:3,4],"u_1");
-## retsys = sysprune(Asys,list("tx","ty","tz"), 4);
+## retsys = sysprune(Asys,@{"tx","ty","tz"@}, 4);
 ## @end example
 ##
 ## @end table
diff --git a/scripts/general/mod.m b/scripts/general/mod.m
--- a/scripts/general/mod.m
+++ b/scripts/general/mod.m
@@ -44,7 +44,8 @@
     usage ("r = mod (x, y)");
   endif
 
-  if (any (size (x) != size (y)) && ! (isscalar (x) || isscalar (y)))
+  if (((ndims (x) != ndims (y)) || any (size (x) != size (y))) &&
+	 ! (isscalar (x) || isscalar (y)))
     error ("mod: argument sizes must agree");
   endif
 
diff --git a/scripts/general/rem.m b/scripts/general/rem.m
--- a/scripts/general/rem.m
+++ b/scripts/general/rem.m
@@ -39,7 +39,8 @@
     usage ("rem (x, y)");
   endif
 
-  if (any (size (x) != size (y)) && ! (isscalar (x) || isscalar (y)))
+  if (((ndims (x) != ndims (y)) || any (size (x) != size (y))) &&
+	 ! (isscalar (x) || isscalar (y)))
     error ("rem: argument sizes must agree");
   endif
 
diff --git a/scripts/general/repmat.m b/scripts/general/repmat.m
--- a/scripts/general/repmat.m
+++ b/scripts/general/repmat.m
@@ -34,21 +34,47 @@
     usage ("repmat (a, m, n)");
   endif
 
-  if (nargin == 2)
+  if (nargin == 3)
+    if (!isscalar (m) && !isscalar (n))
+      error ("repmat: with 3 arguments m and n must be scalar");
+    endif
+    idx = [m, n];
+  else 
     if (isscalar (m))
+      idx = [m, m];
       n = m;
-    elseif (isvector (m) && length (m) == 2)
-      n = m(2);
-      m = m(1);
+    elseif (isvector (m) && length (m) > 1)
+      # Ensure that we have a row vector
+      idx = m(:).';
     else
-      error ("repmat: only builds 2D matrices")
+      error ("repmat: invalid dimensional argument");
     endif
   endif
 
-  if (isstr (a))
-    x = setstr (kron (ones (m, n), toascii (a)));
+  if (numel (a) == 1)
+    if (isstr (a))
+      x = setstr (toascii (a) * ones (idx));
+    else
+      x = a * ones(idx);
+    endif
+  elseif (ndims (a) == 2 && length (idx) < 3)
+    if (isstr (a))
+      x = setstr (kron (ones (idx), toascii (a)));
+    else
+      x = kron (ones (idx), a);
+    endif
   else
-    x = kron (ones (m, n), a);
+    aidx = size(a);
+    if (length(aidx) > length(idx))
+      idx = [idx, ones(1,length(aidx)-length(idx))];
+    elseif (length(aidx) < length(idx))
+      aidx = [aidx, ones(1,length(idx)-length(aidx))];
+    endif
+    cidx = cell ();
+    for i=1:length(aidx)
+      cidx{i} = kron (ones (1, idx(i)), 1:aidx(i));
+    endfor
+    x = a (cidx{:});
   endif
 
 endfunction
diff --git a/scripts/signal/fftshift.m b/scripts/signal/fftshift.m
--- a/scripts/signal/fftshift.m
+++ b/scripts/signal/fftshift.m
@@ -19,6 +19,7 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} fftshift (@var{v})
+## @deftypefnx {Function File} {} fftshift (@var{v}, @var{dim})
 ## Perform a shift of the vector @var{v}, for use with the @code{fft}
 ## and @code{ifft} functions, in order the move the frequency 0 to the
 ## center of the vector or matrix.
@@ -31,32 +32,55 @@
 ## f = ((1:N) - ceil(N/2)) / N / Dt
 ## @end example
 ##
-## If @var{v} is a matrix, the same holds for rows and columns.
+## If @var{v} is a matrix, the same holds for rows and columns. If 
+## @var{v} is an array, then the same holds along each dimension.
+##
+## The optional @var{dim} argument can be used to limit the dimension
+## along which the permutation occurs.
 ## @end deftypefn
 
 ## Author: Vincent Cautaerts <vincent@comf5.comm.eng.osaka-u.ac.jp>
 ## Created: July 1997
 ## Adapted-By: jwe
 
-function retval = fftshift (V)
+function retval = fftshift (V, dim)
 
   retval = 0;
 
-  if (nargin != 1)
-    usage ("usage: fftshift (X)");
+  if (nargin != 1 && nargin != 2)
+    usage ("usage: fftshift (X, dim)");
   endif
 
-  if (isvector (V))
-    x = length (V);
-    xx = ceil (x/2);
-    retval = V([xx+1:x, 1:xx]);
-  elseif (ismatrix (V))
-    [x, y] = size (V);
-    xx = ceil (x/2);
-    yy = ceil (y/2);
-    retval = V([xx+1:x, 1:xx], [yy+1:y, 1:yy]);
+  if (nargin == 2)
+    if (!isscalar (dim))
+      error ("fftshift: dimension must be an integer scalar");
+    endif
+    nd = ndims (V);
+    sz = size (V);
+    sz2 = ceil (sz(dim) / 2);
+    idx = cell ();
+    for i=1:nd
+      idx {i} = 1:sz(i);
+    endfor
+    idx {dim} = [sz2+1:sz(dim), 1:sz2];
+    retval = V (idx{:});
   else
-    error ("fftshift: expecting vector or matrix argument");
+    if (isvector (V))
+      x = length (V);
+      xx = ceil (x/2);
+      retval = V([xx+1:x, 1:xx]);
+    elseif (ismatrix (V))
+      nd = ndims (V);
+      sz = size (V);
+      sz2 = ceil (sz ./ 2);
+      idx = cell ();
+      for i=1:nd
+        idx{i} = [sz2(i)+1:sz(i), 1:sz2(i)];
+      endfor
+      retval = V (idx{:});
+    else
+      error ("fftshift: expecting vector or matrix argument");
+    endif
   endif
 
 endfunction
diff --git a/scripts/specfun/erfinv.m b/scripts/specfun/erfinv.m
--- a/scripts/specfun/erfinv.m
+++ b/scripts/specfun/erfinv.m
@@ -38,9 +38,11 @@
 
   iterations = 0;
 
-  [m, n] = size (x);
-  x = reshape (x, m * n, 1);
-  y = zeros (m * n, 1);
+  sz = size (x);
+  nel = numel (x);
+
+  x = reshape (x, nel, 1);
+  y = zeros (nel, 1);
 
   i = find ((x < -1) | (x > 1) | isnan(x));
   if any (i)
@@ -69,6 +71,6 @@
     y(i) = z_new;
   endif
 
-  y = reshape (y, m, n);
+  y = reshape (y, sz);
 
 endfunction
diff --git a/scripts/statistics/base/center.m b/scripts/statistics/base/center.m
--- a/scripts/statistics/base/center.m
+++ b/scripts/statistics/base/center.m
@@ -19,27 +19,40 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} center (@var{x})
+## @deftypefnx {Function File} {} center (@var{x}, @var{dim})
 ## If @var{x} is a vector, subtract its mean.
 ## If @var{x} is a matrix, do the above for each column.
+## If the optional argument @var{dim} is given, perform the above
+## operation along this dimension
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Description: Center by subtracting means
 
-function retval = center (x)
+function retval = center (x, varargin)
 
-  if (nargin != 1)
+  if (nargin != 1 && nargin != 2)
     usage ("center (x)");
   endif
 
   if (isvector (x))
-    retval = x - mean (x);
+    retval = x - mean (x, varargin{:});
   elseif (ismatrix (x))
-    retval = x - ones (rows (x), 1) * mean (x);
+    if nargin < 2
+      dim = min (find (size (x) > 1));
+      if isempty (dim), 
+	dim=1; 
+      endif;
+    else
+      dim = varargin {1};
+    endif
+    sz = ones (1, ndims (x));
+    sz (dim) = size (x, dim);
+    retval = x - repmat (mean (x, dim), sz);
   elseif (isempty (x))
     retval = x;
   else
     error ("center: x must be a vector or a matrix");
   endif
 
-endfunction
\ No newline at end of file
+endfunction
diff --git a/scripts/statistics/base/meansq.m b/scripts/statistics/base/meansq.m
--- a/scripts/statistics/base/meansq.m
+++ b/scripts/statistics/base/meansq.m
@@ -19,20 +19,22 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} meansq (@var{x})
+## @deftypefnx {Function File} {} meansq (@var{x}, @var{dim})
 ## For vector arguments, return the mean square of the values.
 ## For matrix arguments, return a row vector contaning the mean square
-## of each column.
+## of each column. With the optional @var{dim} argument, returns the
+## mean squared of the values along this dimension
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Description: Compute mean square
 
-function y = meansq (x)
+function y = meansq (x, varargin)
 
-  if (nargin != 1)
+  if (nargin != 1 && nargin != 2)
     usage ("meansq (x)");
   endif
 
-  y = mean (x.^2);
+  y = mean (x.^2, varargin{:});
 
 endfunction
diff --git a/scripts/statistics/base/moment.m b/scripts/statistics/base/moment.m
--- a/scripts/statistics/base/moment.m
+++ b/scripts/statistics/base/moment.m
@@ -18,7 +18,7 @@
 ## 02111-1307, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} moment (@var{x}, @var{p}, @var{opt})
+## @deftypefn {Function File} {} moment (@var{x}, @var{p}, @var{opt},@var{dim})
 ## If @var{x} is a vector, compute the @var{p}-th moment of @var{x}.
 ##
 ## If @var{x} is a matrix, return the row vector containing the
@@ -34,6 +34,9 @@
 ##
 ## @noindent
 ## computes the third central absolute moment of @var{x}.
+##
+## If the optional argument @var{dim} is supplied, work along dimension
+## @var{dim}.
 ## @end deftypefn
 
 ## Can easily be made to work for continuous distributions (using quad)
@@ -42,35 +45,73 @@
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Description: Compute moments
 
-function m = moment (x, p, opt)
+function m = moment (x, p, opt1, opt2)
 
-  if ((nargin < 2) || (nargin > 3))
-    usage ("moment (x, p, type");
+  if ((nargin < 2) || (nargin > 4))
+    usage ("moment (x, p, type, dim)");
   endif
 
-  [nr, nc] = size (x);
-  if (nr == 0 || nc == 0)
-    error ("moment: x must not be empty");
-  elseif (nr == 1)
-    x  = reshape (x, nc, 1);
-    nr = nc;
+  need_dim = 0;
+
+  if (nargin == 1)
+    opt = "";
+    need_dim = 1;
+  elseif (nargin == 2)
+    if (isstr (opt1))
+      opt = opt1;
+      need_dim = 1;
+    else
+      dim = opt1;
+      opt = "";
+    endif
+  elseif (nargin == 3)
+    if (isstr (opt1))
+      opt = opt1;
+      dim = opt2;
+    elseif (isstr (opt2))
+      opt = opt2;
+      dim = opt1;
+    else
+      usage ("moment: expecting opt to be a string");
+    endif
+  else
+    usage ("moment (x, p, dim, opt) or moment (x, p, dim, opt)");
   endif
 
-  if (nargin == 3)
-    tmp = warn_str_to_num;
-    unwind_protect
-      warn_str_to_num = 0;
-      if any (opt == "c")
-	x = x - ones (nr, 1) * sum (x) / nr;
-      endif
-      if any (opt == "a")
-	x = abs (x);
-      endif
-    unwind_protect_cleanup
-      warn_str_to_num = tmp;
-    end_unwind_protect
+  sz = size(x);
+  n = sz (dim);
+
+  if (need_dim)
+    t = find (size (x) != 1);
+    if (isempty (t))
+      dim = 1;
+    else
+      dim = t(1);
+    endif
+  endif
+
+  sz = size (x);
+  n = sz (dim);
+
+  if (numels (x) < 1)
+    error ("moment: x must not be empty");
   endif
 
-  m = sum(x .^ p) / nr;
+  tmp = warn_str_to_num;
+  unwind_protect
+    warn_str_to_num = 0;
+    if any (opt == "c")
+      rng = ones(1, length (sz));
+      rng (dim) = sz (dim);
+      x = x - repmat (sum (x, dim), rng) / n;
+    endif
+    if any (opt == "a")
+      x = abs (x);
+    endif
+  unwind_protect_cleanup
+    warn_str_to_num = tmp;
+  end_unwind_protect
+
+  m = sum(x .^ p, dim) / n;
 
 endfunction
diff --git a/scripts/statistics/base/std.m b/scripts/statistics/base/std.m
--- a/scripts/statistics/base/std.m
+++ b/scripts/statistics/base/std.m
@@ -19,6 +19,8 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} std (@var{x})
+## @deftypefnx {Function File} {} std (@var{x}, @var{opt})
+## @deftypefnx {Function File} {} std (@var{x}, @var{opt}, @var{dim})
 ## If @var{x} is a vector, compute the standard deviation of the elements
 ## of @var{x}.
 ## @iftex
@@ -38,26 +40,50 @@
 ## @end ifinfo
 ## If @var{x} is a matrix, compute the standard deviation for
 ## each column and return them in a row vector.
+##
+## The argument @var{opt} determines the type of normalization to use. Valid values
+## are
+##
+## @table @asis 
+## @item 0:
+##   normalizes with N-1, provides the square root of best unbiased estimator of 
+##   the variance [default]
+## @item 1:
+##   normalizes with N, this provides the square root of the second moment around 
+##   the mean
+## @end table
+##
+## The third argument @var{dim} determines the dimension along which the standard
+## deviation is calculated.
 ## @end deftypefn
 ## @seealso{mean and median}
 
 ## Author: jwe
 
-function retval = std (a)
+function retval = std (a, opt, dim)
 
-  if (nargin != 1)
-    usage ("std (a)");
+  if (nargin < 1 || nargin > 3)
+    usage ("std (a, opt, dim)");
+  endif
+  if nargin < 3
+    dim = min(find(size(a)>1));
+    if isempty(dim), dim=1; endif;
+  endif
+  if ((nargin < 2) || isempty(opt))
+    opt = 0;
   endif
 
-  nr = rows (a);
-  nc = columns (a);
-  if (nc == 1 && nr == 1)
-    retval = 0;
-  elseif (nc == 1 || nr == 1)
-    n = length (a);
-    retval = sqrt (sumsq (a - mean (a)) / (n - 1));
-  elseif (nr > 1 && nc > 0)
-    retval = sqrt (sumsq (a - ones (nr, 1) * mean (a)) / (nr - 1));
+  sz = size(a);
+  if (sz (dim) == 1)
+    retval = zeros(sz);
+  elseif (numel (a) > 0)
+    rng = ones (1, length (sz));
+    rng (dim) = sz (dim);
+    if (opt == 0)
+      retval = sqrt (sumsq (a - repmat(mean (a, dim), rng), dim) / (sz(dim) - 1));
+    else
+      retval = sqrt (sumsq (a - repmat(mean (a, dim), rng), dim) / sz(dim));
+    endif
   else
     error ("std: invalid matrix argument");
   endif
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,26 @@
+2004-04-02  Quentin Spencer  <qspencer@ieee.org>
+
+	* parse.y: Use persistent instead of static in warnings messages.
+
+2004-04-02  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* pt-decl.cc (tree_static_command::do_init): Initialize to empty
+	matrix by default.
+
+2004-04-02  David Bateman  <dbateman@free.fr>
+
+	* ov-re-mat.cc (octave_matrix::convert_to_str_internal):
+	Return charNDArray.
+	* ov-bool-mat.cc (octave_bool_matrix::convert_to_str_internal):
+	Call array_value.
+	
+	* DLD-FUNCTIONS/besselj.cc, DLD-FUNCTIONS/betainc.cc,
+	DLD-FUNCTIONS/gammainc.cc, DLD-FUNCTIONS/minmax.cc,
+	DLD-FUNCTIONS/filter.cc:
+	Convert for NDArray, better Matlab compatibility.
+
+	* load-save.cc (Fload): Better handling of non existent files.
+	
 2004-03-19  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* ov-list.cc (octave_list::subsref): Correctly create return value.
diff --git a/src/DLD-FUNCTIONS/besselj.cc b/src/DLD-FUNCTIONS/besselj.cc
--- a/src/DLD-FUNCTIONS/besselj.cc
+++ b/src/DLD-FUNCTIONS/besselj.cc
@@ -97,6 +97,24 @@
   return retval;
 }
 
+static inline NDArray
+int_arrayN_to_array (const ArrayN<int>& a)
+{
+  dim_vector dv = a.dims ();
+  int nel = dv.numel ();
+
+  NDArray retval (dv);
+
+  for (int i = 0; i < nel; i++)
+    {
+      OCTAVE_QUIT;
+      
+      retval(i) = (double) (a(i));
+    }
+
+  return retval;
+}
+
 static void
 gripe_bessel_arg (const char *fn, const char *arg)
 {
@@ -145,17 +163,17 @@
 		}
 	      else
 		{
-		  ComplexMatrix x = x_arg.complex_matrix_value ();
+		  ComplexNDArray x = x_arg.complex_array_value ();
 
 		  if (! error_state)
 		    {
-		      Array2<int> ierr;
+		      ArrayN<int> ierr;
 		      octave_value result;
 
 		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
 
 		      if (nargout > 1)
-			retval(1) = int_array2_to_matrix (ierr);
+			retval(1) = int_arrayN_to_array (ierr);
 
 		      retval(0) = result;
 		    }
@@ -168,21 +186,28 @@
 	}
       else
 	{
-	  Matrix alpha = args(0).matrix_value ();
+	  dim_vector dv0 = args(0).dims ();
+	  dim_vector dv1 = args(1).dims ();
+	  
+	  bool args0_is_row_vector = (dv0 (1) == dv0.numel ());
+	  bool args1_is_col_vector = (dv1 (0) == dv1.numel ());
 
-	  if (! error_state)
+	  if (args0_is_row_vector && args1_is_col_vector)
 	    {
-	      if (x_arg.is_scalar_type ())
+	      RowVector ralpha = args(0).row_vector_value ();
+
+	      if (! error_state)
 		{
-		  Complex x = x_arg.complex_value ();
+		  ComplexColumnVector cx = 
+		    x_arg.complex_column_vector_value ();
 
 		  if (! error_state)
 		    {
 		      Array2<int> ierr;
 		      octave_value result;
 
-		      DO_BESSEL (type, alpha, x, scaled, ierr, result);
-
+		      DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
+		      
 		      if (nargout > 1)
 			retval(1) = int_array2_to_matrix (ierr);
 
@@ -192,45 +217,56 @@
 		    gripe_bessel_arg (fn, "second");
 		}
 	      else
-		{
-		  ComplexMatrix x = x_arg.complex_matrix_value ();
-
-		  if (! error_state)
-		    {
-		      if (alpha.rows () == 1 && x.columns () == 1)
-			{
-			  RowVector ralpha = alpha.row (0);
-			  ComplexColumnVector cx = x.column (0);
+		gripe_bessel_arg (fn, "first");
+	    }
+	  else
+	    {
+	      NDArray alpha = args(0).array_value ();
 
-			  Array2<int> ierr;
-			  octave_value result;
-
-			  DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
+	      if (! error_state)
+		{
+		  if (x_arg.is_scalar_type ())
+		    {
+		      Complex x = x_arg.complex_value ();
 
-			  if (nargout > 1)
-			    retval(1) = int_array2_to_matrix (ierr);
-
-			  retval(0) = result;
-			}
-		      else
+		      if (! error_state)
 			{
-			  Array2<int> ierr;
+			  ArrayN<int> ierr;
 			  octave_value result;
 
 			  DO_BESSEL (type, alpha, x, scaled, ierr, result);
 
 			  if (nargout > 1)
-			    retval(1) = int_array2_to_matrix (ierr);
+			    retval(1) = int_arrayN_to_array (ierr);
 
 			  retval(0) = result;
 			}
+		      else
+			gripe_bessel_arg (fn, "second");
 		    }
 		  else
-		    gripe_bessel_arg (fn, "second");
+		    {
+		      ComplexNDArray x = x_arg.complex_array_value ();
+
+		      if (! error_state)
+			{
+			  ArrayN<int> ierr;
+			  octave_value result;
+			  
+			  DO_BESSEL (type, alpha, x, scaled, ierr, result);
+			  
+			  if (nargout > 1)
+			    retval(1) = int_arrayN_to_array (ierr);
+			  
+			  retval(0) = result;
+			}
+		      else
+			gripe_bessel_arg (fn, "second");
+		    }
 		}
+	      else
+		gripe_bessel_arg (fn, "first");
 	    }
-	  else
-	    gripe_bessel_arg (fn, "first");
 	}
     }
   else
@@ -422,7 +458,7 @@
 
       int kind = 0;
 
-      ComplexMatrix z;
+      ComplexNDArray z;
 
       if (nargin > 1)
 	{
@@ -441,11 +477,11 @@
 
       if (! error_state)
 	{
-	  z = args(nargin == 1 ? 0 : 1).complex_matrix_value ();
+	  z = args(nargin == 1 ? 0 : 1).complex_array_value ();
 
 	  if (! error_state)
 	    {
-	      Array2<int> ierr;
+	      ArrayN<int> ierr;
 	      octave_value result;
 
 	      if (kind > 1)
@@ -454,7 +490,7 @@
 		result = airy (z, kind == 1, scale, ierr);
 
 	      if (nargout > 1)
-		retval(1) = int_array2_to_matrix (ierr);
+		retval(1) = int_arrayN_to_array (ierr);
 
 	      retval(0) = result;
 	    }
diff --git a/src/DLD-FUNCTIONS/betainc.cc b/src/DLD-FUNCTIONS/betainc.cc
--- a/src/DLD-FUNCTIONS/betainc.cc
+++ b/src/DLD-FUNCTIONS/betainc.cc
@@ -88,7 +88,7 @@
 		    }
 		  else
 		    {
-		      Matrix b = b_arg.matrix_value ();
+		      NDArray b = b_arg.array_value ();
 
 		      if (! error_state)
 			retval = betainc (x, a, b);
@@ -97,7 +97,7 @@
 	    }
 	  else
 	    {
-	      Matrix a = a_arg.matrix_value ();
+	      NDArray a = a_arg.array_value ();
 
 	      if (! error_state)
 		{
@@ -110,7 +110,7 @@
 		    }
 		  else
 		    {
-		      Matrix b = b_arg.matrix_value ();
+		      NDArray b = b_arg.array_value ();
 
 		      if (! error_state)
 			retval = betainc (x, a, b);
@@ -120,7 +120,7 @@
 	}
       else
 	{
-	  Matrix x = x_arg.matrix_value ();
+	  NDArray x = x_arg.array_value ();
 
 	  if (a_arg.is_scalar_type ())
 	    {
@@ -137,7 +137,7 @@
 		    }
 		  else
 		    {
-		      Matrix b = b_arg.matrix_value ();
+		      NDArray b = b_arg.array_value ();
 
 		      if (! error_state)
 			retval = betainc (x, a, b);
@@ -146,7 +146,7 @@
 	    }
 	  else
 	    {
-	      Matrix a = a_arg.matrix_value ();
+	      NDArray a = a_arg.array_value ();
 
 	      if (! error_state)
 		{
@@ -159,7 +159,7 @@
 		    }
 		  else
 		    {
-		      Matrix b = b_arg.matrix_value ();
+		      NDArray b = b_arg.array_value ();
 
 		      if (! error_state)
 			retval = betainc (x, a, b);
diff --git a/src/DLD-FUNCTIONS/filter.cc b/src/DLD-FUNCTIONS/filter.cc
--- a/src/DLD-FUNCTIONS/filter.cc
+++ b/src/DLD-FUNCTIONS/filter.cc
@@ -39,34 +39,28 @@
 #include "oct-obj.h"
 
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-extern MArray<double>
-filter (MArray<double>&, MArray<double>&, MArray<double>&);
+extern MArrayN<double>
+filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
 
-extern MArray<Complex>
-filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&);
+extern MArrayN<Complex>
+filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
 #endif
 
 template <class T>
-MArray<T>
-filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si)
+MArrayN<T>
+filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, MArrayN<T>& si, 
+	int dim = 0)
 {
-  MArray<T> y;
+  MArrayN<T> y;
 
   int a_len  = a.length ();
   int b_len  = b.length ();
-  int x_len  = x.length ();
-
-  int si_len = si.length ();
 
   int ab_len = a_len > b_len ? a_len : b_len;
 
   b.resize (ab_len, 0.0);
-
-  if (si.length () != ab_len - 1)
-    {
-      error ("filter: si must be a vector of length max (length (a), length (b)) - 1");
-      return y;
-    }
+  if (a_len > 1)
+    a.resize (ab_len, 0.0);
 
   T norm = a (0);
 
@@ -76,95 +70,183 @@
       return y;
     }
 
-  y.resize (x_len, 0.0);
+  dim_vector x_dims = x.dims ();
+  if ((dim < 0) || (dim > x_dims.length ()))
+    {
+      error ("filter: filtering over invalid dimension");
+      return y;
+    }
+
+  int x_len = x_dims (dim);
+
+  dim_vector si_dims = si.dims ();
+  int si_len = si_dims (0);
+
+  if (si_len != ab_len - 1)
+    {
+      error ("filter: first dimension of si must be of length max (length (a), length (b)) - 1");
+      return y;
+    }
+
+  if (si_dims.length () != x_dims.length ())
+    {
+      error ("filter: dimensionality of si and x must agree");
+      return y;
+    }
+
+  int si_dim = 0;
+  for (int i = 0; i < x_dims.length (); i++)
+    {
+      if (i == dim)
+	continue;
+      
+      if (si_dims (++si_dim) != x_dims (i))
+	{
+	  error ("filter: dimensionality of si and x must agree");
+	  return y;
+	}
+    }
 
   if (norm != 1.0)
-    b = b / norm;
+    {
+      a = a / norm;
+      b = b / norm;
+    }
 
-  if (a_len > 1)
-    {
-      a.resize (ab_len, 0.0);
+  if ((a_len <= 1) && (si_len <= 1))
+    return b(0) * x;
+
+  y.resize (x_dims, 0.0);
+
+  int x_stride = 1;
+  for (int i = 0; i < dim; i++)
+    x_stride *= x_dims(i);
 
-      if (norm != 1.0)
-	a = a / norm;
+  int x_num = x_dims.numel () / x_len;
+  for (int num = 0; num < x_num; num++)
+    {
+      int x_offset;
+      if (x_stride == 1)
+	x_offset = num * x_len;
+      else
+	{
+	  int x_offset2 = 0;
+	  x_offset = num;
+	  while (x_offset >= x_stride)
+	    {
+	      x_offset -= x_stride;
+	      x_offset2++;
+	    }
+	  x_offset += x_offset2 * x_stride * x_len;
+	}
+      int si_offset = num * si_len;
 
-      for (int i = 0; i < x_len; i++)
+      if (a_len > 1)
 	{
-	  y (i) = si (0) + b (0) * x (i);
-
-	  if (si_len > 1)
+	  for (int i = 0; i < x_len; i++)
 	    {
-	      for (int j = 0; j < si_len - 1; j++)
+	      int idx = i * x_stride + x_offset; 
+	      y (idx) = si (si_offset) + b (0) * x (idx);
+
+	      if (si_len > 1)
 		{
-		  OCTAVE_QUIT;
+		  for (int j = 0; j < si_len - 1; j++)
+		    {
+		      OCTAVE_QUIT;
 
-		  si (j) = si (j+1) - a (j+1) * y (i)
-		    + b (j+1) * x (i);
+		      si (j + si_offset) =  si (j + 1 + si_offset) - 
+			a (j+1) * y (idx) + b (j+1) * x (idx);
+		    }
+
+		  si (si_len - 1 + si_offset) = b (si_len) * x (idx)
+		    - a (si_len) * y (idx);
 		}
-
-	      si (si_len-1) = b (si_len) * x (i)
-		- a (si_len) * y (i);
+	      else
+		si (si_offset) = b (si_len) * x (idx)
+		  - a (si_len) * y (idx);
 	    }
-	  else
-	    si (0) = b (si_len) * x (i)
-	      - a (si_len) * y (i);
+	}
+      else if (si_len > 0)
+	{
+	  for (int i = 0; i < x_len; i++)
+	    {
+	      int idx = i * x_stride + x_offset; 
+	      y (idx) = si (si_offset) + b (0) * x (idx);
+
+	      if (si_len > 1)
+		{
+		  for (int j = 0; j < si_len - 1; j++)
+		    {
+		      OCTAVE_QUIT;
+
+		      si (j + si_offset) = si (j + 1 + si_offset) + 
+			b (j+1) * x (idx);
+		    }
+
+		  si (si_len - 1 + si_offset) = b (si_len) * x (idx);
+		}
+	      else
+		si (si_offset) = b (1) * x (idx);
+	    }
 	}
     }
-  else if (si_len > 0)
-    {
-      for (int i = 0; i < x_len; i++)
-	{
-	  y (i) = si (0) + b (0) * x (i);
-
-	  if (si_len > 1)
-	    {
-	      for (int j = 0; j < si_len - 1; j++)
-		{
-		  OCTAVE_QUIT;
-
-		  si (j) = si (j+1) + b (j+1) * x (i);
-		}
-
-	      si (si_len-1) = b (si_len) * x (i);
-	    }
-	  else
-	    si (0) = b (1) * x (i);
-	}
-    }
-  else
-    y = b (0) * x;
-
+  
   return y;
 }
 
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-extern MArray<double>
-filter (MArray<double>&, MArray<double>&, MArray<double>&,
-	MArray<double>&);
+extern MArrayN<double>
+filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
+	MArrayN<double>&, int dim);
 
-extern MArray<Complex>
-filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
-	MArray<Complex>&);
+extern MArrayN<Complex>
+filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
+	MArrayN<Complex>&, int dim);
 #endif
 
 template <class T>
-MArray<T>
-filter (MArray<T>& b, MArray<T>& a, MArray<T>& x)
+MArrayN<T>
+filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, int dim = -1)
 {
+  dim_vector x_dims = x.dims ();
+
+  if (dim < 0)
+    {
+      // Find first non-singleton dimension
+      while ((dim < x_dims.length()) && (x_dims (dim) <= 1))
+	dim++;
+  
+      // All dimensions singleton, pick first dimension
+      if (dim == x_dims.length ())
+	dim = 0;
+    }
+  else
+    if (dim < 0 || dim > x_dims.length ())
+      {
+	error ("filter: filtering over invalid dimension");
+	return MArrayN<T> ();
+      }
+
   int a_len = a.length ();
   int b_len = b.length ();
 
   int si_len = (a_len > b_len ? a_len : b_len) - 1;
+  dim_vector si_dims = x.dims ();
+  for (int i = dim; i > 0; i--)
+    si_dims (i) = si_dims (i-1);
+  si_dims (0) = si_len;
+  
+  MArrayN<T> si (si_dims, T (0.0));
 
-  MArray<T> si (si_len, T (0.0));
-
-  return filter (b, a, x, si);
+  return filter (b, a, x, si, dim);
 }
 
 DEFUN_DLD (filter, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
 @deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
+@deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})\n\
+@deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})\n\
 Return the solution to the following linear, time-invariant difference\n\
 equation:\n\
 @iftex\n\
@@ -194,7 +276,8 @@
  $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
 @end tex\n\
 @end iftex\n\
-An equivalent form of this equation is:\n\
+over the first non-singleton dimension of @var{x} or over @var{dim} if\n\
+supplied. An equivalent form of this equation is:\n\
 @iftex\n\
 @tex\n\
 $$\n\
@@ -259,59 +342,91 @@
 
   int nargin  = args.length ();
 
-  if (nargin < 3 || nargin > 4)
+  if (nargin < 3 || nargin > 5)
     {
       print_usage ("filter");
       return retval;
     }
 
-  const char *errmsg = "filter: arguments must be vectors";
+  const char *errmsg = "filter: arguments a and b must be vectors";
+
+  int dim;
+  dim_vector x_dims = args(2).dims ();
 
-  bool x_is_row_vector = (args(2).rows () == 1);
-
-  bool si_is_row_vector = (nargin == 4 && args(3).rows () == 1);
+  if (nargin == 5)
+    {
+      dim = args(4).nint_value() - 1;
+      if (dim < 0 || dim >= x_dims.length ())
+	{
+	  error ("filter: filtering over invalid dimension");
+	  return retval;
+	}
+    }
+  else
+    {
+      // Find first non-singleton dimension
+      dim = 0;
+      while ((dim < x_dims.length()) && (x_dims (dim) <= 1))
+	dim++;
+  
+      // All dimensions singleton, pick first dimension
+      if (dim == x_dims.length ())
+	dim = 0;
+    }
 
   if (args(0).is_complex_type ()
       || args(1).is_complex_type ()
       || args(2).is_complex_type ()
-      || (nargin == 4 && args(3).is_complex_type ()))
+      || (nargin >= 4 && args(3).is_complex_type ()))
     {
       ComplexColumnVector b (args(0).complex_vector_value ());
       ComplexColumnVector a (args(1).complex_vector_value ());
-      ComplexColumnVector x (args(2).complex_vector_value ());
+
+      ComplexNDArray x (args(2).complex_array_value ());
 
       if (! error_state)
 	{
-	  ComplexColumnVector si;
+	  ComplexNDArray si;
 
-	  if (nargin == 3)
+	  if (nargin == 3 || args(3).is_empty ())
 	    {
 	      int a_len = a.length ();
 	      int b_len = b.length ();
 
 	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
 
-	      si.resize (si_len, 0.0);
+	      dim_vector si_dims = x.dims ();
+	      for (int i = dim; i > 0; i--)
+		si_dims (i) = si_dims (i-1);
+	      si_dims (0) = si_len;
+
+	      si.resize (si_dims, 0.0);
 	    }
 	  else
-	    si = ComplexColumnVector (args(3).complex_vector_value ());
+	    {
+	      dim_vector si_dims = args (3).dims ();
+	      bool si_is_vector = true;
+	      for (int i=0; i < si_dims.length (); i++)
+		if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ()))
+		  {
+		    si_is_vector = false;
+		    break;
+		  }
+
+	      if (si_is_vector)
+		si = ComplexNDArray (args(3).complex_vector_value ());
+	      else
+		si = args(3).complex_array_value ();
+	    }
 
 	  if (! error_state)
 	    {
-	      ComplexColumnVector y (filter (b, a, x, si));
+	      ComplexNDArray y (filter (b, a, x, si, dim));
 
 	      if (nargout == 2)
-		{
-		  if (si_is_row_vector)
-		    retval(1) = si.transpose ();
-		  else
-		    retval(1) = si;
-		}
+		retval(1) = si;
 
-	      if (x_is_row_vector)
-		retval(0) = y.transpose ();
-	      else
-		retval(0) = y;
+	      retval(0) = y;
 	    }
 	  else
 	    error (errmsg);
@@ -323,40 +438,52 @@
     {
       ColumnVector b (args(0).vector_value ());
       ColumnVector a (args(1).vector_value ());
-      ColumnVector x (args(2).vector_value ());
+
+      NDArray x (args(2).array_value ());
 
       if (! error_state)
 	{
-	  ColumnVector si;
+	  NDArray si;
 
-	  if (nargin == 3)
+	  if (nargin == 3 || args(3).is_empty ())
 	    {
 	      int a_len = a.length ();
 	      int b_len = b.length ();
 
 	      int si_len = (a_len > b_len ? a_len : b_len) - 1;
 
-	      si.resize (si_len, 0.0);
+	      dim_vector si_dims = x.dims ();
+	      for (int i = dim; i > 0; i--)
+		si_dims (i) = si_dims (i-1);
+	      si_dims (0) = si_len;
+
+	      si.resize (si_dims, 0.0);
 	    }
 	  else
-	    si = ColumnVector (args(3).vector_value ());
+	    {
+	      dim_vector si_dims = args (3).dims ();
+	      bool si_is_vector = true;
+	      for (int i=0; i < si_dims.length (); i++)
+		if ((si_dims (i) != 1) && (si_dims (i) < si_dims.numel ()))
+		  {
+		    si_is_vector = false;
+		    break;
+		  }
+
+	      if (si_is_vector)
+		si = NDArray (args(3).vector_value ());
+	      else
+		si = args(3).array_value ();
+	    }
 
 	  if (! error_state)
 	    {
-	      ColumnVector y (filter (b, a, x, si));
+	      NDArray y (filter (b, a, x, si, dim));
 
 	      if (nargout == 2)
-		{
-		  if (si_is_row_vector)
-		    retval(1) = si.transpose ();
-		  else
-		    retval(1) = si;
-		}
+		retval(1) = si;
 
-	      if (x_is_row_vector)
-		retval(0) = y.transpose ();
-	      else
-		retval(0) = y;
+	      retval(0) = y;
 	    }
 	  else
 	    error (errmsg);
@@ -368,19 +495,19 @@
   return retval;
 }
 
-template MArray<double>
-filter (MArray<double>&, MArray<double>&, MArray<double>&,
-	MArray<double>&);
+template MArrayN<double>
+filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
+	MArrayN<double>&, int dim);
 
-template MArray<double>
-filter (MArray<double>&, MArray<double>&, MArray<double>&);
+template MArrayN<double>
+filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
 
-template MArray<Complex>
-filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
-	MArray<Complex>&);
+template MArrayN<Complex>
+filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
+	MArrayN<Complex>&, int dim);
 
-template MArray<Complex>
-filter (MArray<Complex>&, MArray<Complex>&, MArray <Complex>&);
+template MArrayN<Complex>
+filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
 
 /*
 ;;; Local Variables: ***
diff --git a/src/DLD-FUNCTIONS/gammainc.cc b/src/DLD-FUNCTIONS/gammainc.cc
--- a/src/DLD-FUNCTIONS/gammainc.cc
+++ b/src/DLD-FUNCTIONS/gammainc.cc
@@ -86,7 +86,7 @@
 		}
 	      else
 		{
-		  Matrix a = a_arg.matrix_value ();
+		  NDArray a = a_arg.array_value ();
 
 		  if (! error_state)
 		    retval = gammainc (x, a);
@@ -95,7 +95,7 @@
 	}
       else
 	{
-	  Matrix x = x_arg.matrix_value ();
+	  NDArray x = x_arg.array_value ();
 
 	  if (! error_state)
 	    {
@@ -108,7 +108,7 @@
 		}
 	      else
 		{
-		  Matrix a = a_arg.matrix_value ();
+		  NDArray a = a_arg.array_value ();
 
 		  if (! error_state)
 		    retval = gammainc (x, a);
diff --git a/src/DLD-FUNCTIONS/minmax.cc b/src/DLD-FUNCTIONS/minmax.cc
--- a/src/DLD-FUNCTIONS/minmax.cc
+++ b/src/DLD-FUNCTIONS/minmax.cc
@@ -28,8 +28,8 @@
 
 #include "lo-ieee.h"
 #include "lo-mappers.h"
-#include "dMatrix.h"
-#include "CMatrix.h"
+#include "dNDArray.h"
+#include "CNDArray.h"
 #include "quit.h"
 
 #include "defun-dld.h"
@@ -37,13 +37,15 @@
 #include "gripes.h"
 #include "oct-obj.h"
 
+#include "ov-cx-mat.h"
+
 #define MINMAX_BODY(FCN) \
  \
   octave_value_list retval;  \
  \
   int nargin = args.length (); \
  \
-  if (nargin < 1 || nargin > 2 || nargout > 2) \
+  if (nargin < 1 || nargin > 3 || nargout > 2) \
     { \
       print_usage (#FCN); \
       return retval; \
@@ -51,9 +53,13 @@
  \
   octave_value arg1; \
   octave_value arg2; \
+  octave_value arg3; \
  \
   switch (nargin) \
     { \
+    case 3: \
+      arg3 = args(2); \
+ \
     case 2: \
       arg2 = args(1); \
  \
@@ -66,97 +72,111 @@
       break; \
     } \
  \
-  if (nargin == 1 && (nargout == 1 || nargout == 0)) \
+  int dim; \
+  dim_vector dv = ((const octave_complex_matrix&) arg1) .dims (); \
+  if (error_state) \
+    { \
+      gripe_wrong_type_arg (#FCN, arg1);  \
+      return retval; \
+    } \
+ \
+  if (nargin == 3) \
+    { \
+      dim = arg3.nint_value () - 1;  \
+      if (dim < 0 || dim >= dv.length ()) \
+        { \
+	  error ("%s: invalid dimension", #FCN); \
+	  return retval; \
+	} \
+    } \
+  else \
+    { \
+      dim = 0; \
+      while ((dim < dv.length ()) && (dv (dim) <= 1)) \
+	dim++; \
+      if (dim == dv.length ()) \
+	dim = 0; \
+    } \
+ \
+  bool single_arg = (nargin == 1) || arg2.is_empty();	\
+ \
+  if (single_arg) \
+    { \
+      dv(dim) = 1; \
+      int n_dims = dv.length (); \
+      for (int i = n_dims; i > 1; i--) \
+	{ \
+	  if (dv(i-1) == 1) \
+	    n_dims--; \
+	  else \
+	    break; \
+	} \
+      dv.resize (n_dims); \
+    } \
+ \
+  if (single_arg && (nargout == 1 || nargout == 0)) \
     { \
       if (arg1.is_real_type ()) \
 	{ \
-	  Matrix m = arg1.matrix_value (); \
+	  NDArray m = arg1.array_value (); \
  \
 	  if (! error_state) \
 	    { \
-	      if (m.rows () == 1) \
-		retval(0) = m.row_ ## FCN (); \
-	      else \
-		{ \
-		  if (m.rows () == 0 || m.columns () == 0) \
-		    retval(0) = Matrix (); \
-		  else \
-		    retval(0) = m.column_ ## FCN (); \
-		} \
+	      NDArray n = m. FCN (dim); \
+	      n.resize (dv); \
+	      retval(0) = n; \
 	    } \
 	} \
       else if (arg1.is_complex_type ()) \
 	{ \
-	  ComplexMatrix m = arg1.complex_matrix_value (); \
+	  ComplexNDArray m = arg1.complex_array_value (); \
  \
 	  if (! error_state) \
 	    { \
-	      if (m.rows () == 1) \
-		retval(0) = m.row_ ## FCN (); \
-	      else \
-		{ \
-		  if (m.rows () == 0 || m.columns () == 0) \
-		    retval(0) = Matrix (); \
-		  else \
-		    retval(0) = m.column_ ## FCN (); \
-		} \
+	      ComplexNDArray n = m. FCN (dim); \
+	      n.resize (dv); \
+	      retval(0) = n; \
 	    } \
 	} \
       else \
 	gripe_wrong_type_arg (#FCN, arg1); \
     } \
-  else if (nargin == 1 && nargout == 2) \
+  else if (single_arg && nargout == 2) \
     { \
-      Array<int> index; \
+      ArrayN<int> index; \
  \
       if (arg1.is_real_type ()) \
 	{ \
-	  Matrix m = arg1.matrix_value (); \
+	  NDArray m = arg1.array_value (); \
  \
 	  if (! error_state) \
 	    { \
-	      retval.resize (2); \
- \
-	      if (m.rows () == 1) \
-		retval(0) = m.row_ ## FCN (index); \
-	      else \
-		{ \
-		  if (m.rows () == 0 || m.columns () == 0) \
-		    retval(0) = Matrix (); \
-		  else \
-		    retval(0) = m.column_ ## FCN (index); \
-		} \
+	      NDArray n = m. FCN (index, dim);	\
+	      n.resize (dv); \
+	      retval(0) = n; \
 	    } \
 	} \
       else if (arg1.is_complex_type ()) \
 	{ \
-	  ComplexMatrix m = arg1.complex_matrix_value (); \
+	  ComplexNDArray m = arg1.complex_array_value (); \
  \
 	  if (! error_state) \
 	    { \
-	      retval.resize (2); \
- \
-	      if (m.rows () == 1) \
-		retval(0) = m.row_ ## FCN (index); \
-	      else \
-		{ \
-		  if (m.rows () == 0 || m.columns () == 0) \
-		    retval(0) = Matrix (); \
-		  else \
-		    retval(0) = m.column_ ## FCN (index); \
-		} \
+	      ComplexNDArray n = m. FCN (index, dim);	\
+	      n.resize (dv); \
+	      retval(0) = n; \
 	    } \
 	} \
       else \
 	gripe_wrong_type_arg (#FCN, arg1); \
  \
-      int len = index.length (); \
+      int len = index.numel (); \
  \
       if (len > 0) \
 	{ \
 	  double nan_val = lo_ieee_nan_value (); \
  \
-	  RowVector idx (len); \
+	  NDArray idx (index.dims ()); \
  \
 	  for (int i = 0; i < len; i++) \
 	    { \
@@ -169,9 +189,9 @@
 	  retval(1) = idx; \
 	} \
       else \
-	retval(1) = Matrix (); \
+	retval(1) = NDArray (); \
     } \
-  else if (nargin == 2) \
+  else \
     { \
       int arg1_is_scalar = arg1.is_scalar_type (); \
       int arg2_is_scalar = arg2.is_scalar_type (); \
@@ -184,10 +204,10 @@
 	  if (arg1_is_complex || arg2_is_complex) \
 	    { \
 	      Complex c1 = arg1.complex_value (); \
-	      ComplexMatrix m2 = arg2.complex_matrix_value (); \
+	      ComplexNDArray m2 = arg2.complex_array_value (); \
 	      if (! error_state) \
 		{ \
-		  ComplexMatrix result = FCN (c1, m2); \
+		  ComplexNDArray result = FCN (c1, m2); \
 		  if (! error_state) \
 		    retval(0) = result; \
 		} \
@@ -195,11 +215,11 @@
 	  else \
 	    { \
 	      double d1 = arg1.double_value (); \
-	      Matrix m2 = arg2.matrix_value (); \
+	      NDArray m2 = arg2.array_value (); \
  \
 	      if (! error_state) \
 		{ \
-		  Matrix result = FCN (d1, m2); \
+		  NDArray result = FCN (d1, m2); \
 		  if (! error_state) \
 		    retval(0) = result; \
 		} \
@@ -209,24 +229,24 @@
 	{ \
 	  if (arg1_is_complex || arg2_is_complex) \
 	    { \
-	      ComplexMatrix m1 = arg1.complex_matrix_value (); \
+	      ComplexNDArray m1 = arg1.complex_array_value (); \
  \
 	      if (! error_state) \
 		{ \
 		  Complex c2 = arg2.complex_value (); \
-		  ComplexMatrix result = FCN (m1, c2); \
+		  ComplexNDArray result = FCN (m1, c2); \
 		  if (! error_state) \
 		    retval(0) = result; \
 		} \
 	    } \
 	  else \
 	    { \
-	      Matrix m1 = arg1.matrix_value (); \
+	      NDArray m1 = arg1.array_value (); \
  \
 	      if (! error_state) \
 		{ \
 		  double d2 = arg2.double_value (); \
-		  Matrix result = FCN (m1, d2); \
+		  NDArray result = FCN (m1, d2); \
 		  if (! error_state) \
 		    retval(0) = result; \
 		} \
@@ -236,15 +256,15 @@
 	{ \
 	  if (arg1_is_complex || arg2_is_complex) \
 	    { \
-	      ComplexMatrix m1 = arg1.complex_matrix_value (); \
+	      ComplexNDArray m1 = arg1.complex_array_value (); \
  \
 	      if (! error_state) \
 		{ \
-		  ComplexMatrix m2 = arg2.complex_matrix_value (); \
+		  ComplexNDArray m2 = arg2.complex_array_value (); \
  \
 		  if (! error_state) \
 		    { \
-		      ComplexMatrix result = FCN (m1, m2); \
+		      ComplexNDArray result = FCN (m1, m2); \
 		      if (! error_state) \
 			retval(0) = result; \
 		    } \
@@ -252,15 +272,15 @@
 	    } \
 	  else \
 	    { \
-	      Matrix m1 = arg1.matrix_value (); \
+	      NDArray m1 = arg1.array_value (); \
  \
 	      if (! error_state) \
 		{ \
-		  Matrix m2 = arg2.matrix_value (); \
+		  NDArray m2 = arg2.array_value (); \
  \
 		  if (! error_state) \
 		    { \
-		      Matrix result = FCN (m1, m2); \
+		      NDArray result = FCN (m1, m2); \
 		      if (! error_state) \
 			retval(0) = result; \
 		    } \
@@ -268,21 +288,18 @@
 	    } \
 	} \
     } \
-  else \
-    panic_impossible (); \
  \
   return retval
 
 DEFUN_DLD (min, args, nargout,
   "-*- texinfo -*-\n\
-@deftypefn {Mapping Function} {} min (@var{x}, @var{y})\n\
+@deftypefn {Mapping Function} {} min (@var{x}, @var{y}, @var{dim})\n\
 @deftypefnx {Mapping Function} {[@var{w}, @var{iw}] =} min (@var{x})\n\
 @cindex Utility Functions\n\
 For a vector argument, return the minimum value.  For a matrix\n\
 argument, return the minimum value from each column, as a row\n\
-vector.\n\
-For two matrices (or a matrix and scalar),\n\
-return the pair-wise minimum.\n\
+vector, or over the dimension @var{dim} if defined. For two matrices\n\
+(or a matrix and scalar), return the pair-wise minimum.\n\
 Thus,\n\
 \n\
 @example\n\
@@ -323,14 +340,13 @@
 
 DEFUN_DLD (max, args, nargout,
   "-*- texinfo -*-\n\
-@deftypefn {Mapping Function} {} max (@var{x}, @var{y})\n\
+@deftypefn {Mapping Function} {} max (@var{x}, @var{y}, @var{dim})\n\
 @deftypefnx {Mapping Function} {[@var{w}, @var{iw}] =} max (@var{x})\n\
 @cindex Utility Functions\n\
 For a vector argument, return the maximum value.  For a matrix\n\
 argument, return the maximum value from each column, as a row\n\
-vector.\n\
-For two matrices (or a matrix and scalar),\n\
-return the pair-wise maximum.\n\
+vector, or over the dimension @var{dim} if defined. For two matrices\n\
+(or a matrix and scalar), return the pair-wise maximum.\n\
 Thus,\n\
 \n\
 @example\n\
diff --git a/src/load-save.cc b/src/load-save.cc
--- a/src/load-save.cc
+++ b/src/load-save.cc
@@ -306,6 +306,13 @@
 {
   load_save_format retval = LS_UNKNOWN;
 
+  // If the file doesn't exist do nothing
+  std::ifstream file_exist (fname.c_str ());
+  if (file_exist)
+    file_exist.close ();
+  else
+    return LS_UNKNOWN;
+
 #ifdef HAVE_HDF5
   // check this before we open the file
   if (H5Fis_hdf5 (fname.c_str ()) > 0)
@@ -722,19 +729,26 @@
 	{
 	  i++;
 
-	  hdf5_ifstream hdf5_file (fname.c_str ());
-
-	  if (hdf5_file.file_id >= 0)
+	  // If the file doesn't exist do nothing
+	  std::ifstream file (fname.c_str (), std::ios::in);
+	  if (file)
 	    {
-	      retval = do_load (hdf5_file, orig_fname, force, format,
-				flt_fmt, list_only, swap, verbose,
-				argv, i, argc, nargout);
+	      file.close ();
+	      
+	      hdf5_ifstream hdf5_file (fname.c_str ());
 
-	      hdf5_file.close ();
+	      if (hdf5_file.file_id >= 0)
+		{
+		  retval = do_load (hdf5_file, orig_fname, force, format,
+				    flt_fmt, list_only, swap, verbose,
+				    argv, i, argc, nargout);
+
+		  hdf5_file.close ();
+		}
+	      else
+		error ("load: couldn't open input file `%s'",
+		       orig_fname.c_str ());
 	    }
-	  else
-	    error ("load: couldn't open input file `%s'",
-		   orig_fname.c_str ());
 	}
       else
 #endif /* HAVE_HDF5 */
diff --git a/src/ov-bool-mat.cc b/src/ov-bool-mat.cc
--- a/src/ov-bool-mat.cc
+++ b/src/ov-bool-mat.cc
@@ -142,7 +142,7 @@
 octave_value
 octave_bool_matrix::convert_to_str_internal (bool pad, bool force) const
 {
-  octave_value tmp = octave_value (matrix_value ());
+  octave_value tmp = octave_value (array_value ());
   return tmp.convert_to_str (pad, force);
 }
 
diff --git a/src/ov-re-mat.cc b/src/ov-re-mat.cc
--- a/src/ov-re-mat.cc
+++ b/src/ov-re-mat.cc
@@ -186,66 +186,54 @@
 octave_matrix::convert_to_str_internal (bool, bool) const
 {
   octave_value retval;
+  dim_vector dv = dims ();
+  int nel = dv.numel ();
 
-  int nr = matrix.rows ();
-  int nc = matrix.columns ();
-
-  if (nr == 0 && nc == 0)
+  if (nel == 0)
     {
       char s = '\0';
       retval = octave_value (&s);
     }
   else
     {
-      if (nr == 0 || nc == 0)
-	{
-	  char s = '\0';
-	  retval = octave_value (&s);
-	}
-      else
-	{
-	  charMatrix chm (nr, nc);
+      charNDArray chm (dv);
 	  
-	  bool warned = false;
+      bool warned = false;
+
+      for (int i = 0; i < nel; i++)
+	{
+	  OCTAVE_QUIT;
+
+	  double d = matrix (i);
 
-	  for (int j = 0; j < nc; j++)
+	  if (xisnan (d))
 	    {
-	      for (int i = 0; i < nr; i++)
+	      ::error ("invalid conversion from NaN to character");
+	      return retval;
+	    }
+	  else
+	    {
+	      int ival = NINT (d);
+
+	      if (ival < 0 || ival > UCHAR_MAX)
 		{
-		  OCTAVE_QUIT;
+		  // XXX FIXME XXX -- is there something
+		  // better we could do?
 
-		  double d = matrix (i, j);
+		  ival = 0;
 
-		  if (xisnan (d))
+		  if (! warned)
 		    {
-		      ::error ("invalid conversion from NaN to character");
-		      return retval;
-		    }
-		  else
-		    {
-		      int ival = NINT (d);
-
-		      if (ival < 0 || ival > UCHAR_MAX)
-			{
-			  // XXX FIXME XXX -- is there something
-			  // better we could do?
-
-			  ival = 0;
-
-			  if (! warned)
-			    {
-			      ::warning ("range error for conversion to character value");
-			      warned = true;
-			    }
-			}
-
-		      chm (i, j) = static_cast<char> (ival);
+		      ::warning ("range error for conversion to character value");
+		      warned = true;
 		    }
 		}
-	    }
 
-	  retval = octave_value (chm, 1);
+	      chm (i) = static_cast<char> (ival);
+	    }
 	}
+
+      retval = octave_value (chm, 1);
     }
 
   return retval;
diff --git a/src/parse.y b/src/parse.y
--- a/src/parse.y
+++ b/src/parse.y
@@ -2810,10 +2810,10 @@
       else
 	{
 	  if (reading_script_file)
-	    warning ("ignoring static declaration near line %d of file `%s'",
+	    warning ("ignoring persistent declaration near line %d of file `%s'",
 		     l, curr_fcn_file_full_name.c_str ());
 	  else
-	    warning ("ignoring static declaration near line %d", l);
+	    warning ("ignoring persistent declaration near line %d", l);
 	}
       break;
 
diff --git a/src/pt-decl.cc b/src/pt-decl.cc
--- a/src/pt-decl.cc
+++ b/src/pt-decl.cc
@@ -147,13 +147,18 @@
     {
       id->mark_as_static ();
 
-      tree_expression *expr = elt.expression ();
+      octave_lvalue ult = id->lvalue ();
 
-      if (expr)
+      if (ult.is_defined ())
 	{
-	  octave_value init_val = expr->rvalue ();
+	  tree_expression *expr = elt.expression ();
+
+	  octave_value init_val;
 
-	  octave_lvalue ult = id->lvalue ();
+	  if (expr)
+	    init_val = expr->rvalue ();
+	  else
+	    init_val = Matrix ();
 
 	  ult.assign (octave_value::op_asn_eq, init_val);
 	}
diff --git a/test/octave.test/arith/max-4.m b/test/octave.test/arith/max-4.m
--- a/test/octave.test/arith/max-4.m
+++ b/test/octave.test/arith/max-4.m
@@ -1,1 +1,1 @@
-max (1, 2, 3)
+max (1, 2, 3, 4)
diff --git a/test/octave.test/arith/min-4.m b/test/octave.test/arith/min-4.m
--- a/test/octave.test/arith/min-4.m
+++ b/test/octave.test/arith/min-4.m
@@ -1,1 +1,1 @@
-min (1, 2, 3)
+min (1, 2, 3, 4)
diff --git a/test/octave.test/stats/std-3.m b/test/octave.test/stats/std-3.m
--- a/test/octave.test/stats/std-3.m
+++ b/test/octave.test/stats/std-3.m
@@ -1,1 +1,1 @@
-std (1, 2)
+std (1, 2, 3, 4)