changeset 9042:97aa01a85ea4

Merge documentation cleanup changes to main branch.
author Rik <rdrider0-list@yahoo.com>
date Wed, 25 Mar 2009 18:13:08 -0700
parents 12ca81f1fa99 (diff) 853f96e8008f (current diff)
children cfad8f9a77fa
files NEWS scripts/plot/print.m src/DLD-FUNCTIONS/cellfun.cc src/DLD-FUNCTIONS/find.cc src/parse.y
diffstat 94 files changed, 1801 insertions(+), 517 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,12 @@
+2009-03-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* NEWS: More updates.
+
+2009-03-20  Jaroslav Hajek  <highegg@gmail.com>
+
+	* aclocal.m4 (OCTAVE_CMATH_FUNC): New macro.
+	* configure.in: Use it.
+
 2009-03-09  John W. Eaton  <jwe@octave.org>
 
 	* run-octave.in: Use doc-cache instead of DOC for doc cache file.
--- a/NEWS
+++ b/NEWS
@@ -202,6 +202,39 @@
     The performance of the sum, prod, sumsq, cumsum, cumprod, any, all,
     max and min functions has been significantly improved.
 
+ ** Sorting and searching.
+    
+    The performance of sort has been improved, especially when sorting
+    indices are requested. An efficient built-in issorted implementation
+    was added. sortrows now uses a more efficient algorithm, especially
+    in the homegeneous case. lookup is now a built-in function performing
+    a binary search, optimized for long runs of close elements. Lookup
+    also works with cell arrays of strings.
+
+ ** Range arithmetics
+
+    For some operations on ranges, Octave will attempt to keep the result as a
+    range.  These include negation, adding a scalar, subtracting a scalar, and
+    multiplying by a scalar. Ranges with zero increment are allowed and can be
+    constructed using the built-in function `ones'.
+
+ ** Various performance improvements.
+
+    Performance of a number of other built-in operations and functions was
+    improved, including:
+
+    * logical operations
+    * comparison operators
+    * element-wise power
+    * accumarray
+    * cellfun
+    * isnan
+    * isinf
+    * isfinite
+    * nchoosek
+    * repmat
+    * strcmp
+
  ** 64-bit integer arithmetic.
 
     Arithmetic with 64-bit integers (int64 and uint64 types) is fully
@@ -262,21 +295,22 @@
 
       addtodate          hypot                       reallog
       bicgstab           idivide                     realpow
-      cgs                info                        realsqrt
-      command_line_path  interp1q                    rectint
-      contrast           isdebugmode                 regexptranslate
-      convn              isfloat                     restoredefaultpath
-      cummin             isstrprop                   roundb
-      cummax             log1p                       rundemos
-      datetick           lsqnonneg                   runlength
-      display            matlabroot                  saveobj
-      expm1              namelengthmax               spaugment
-      filemarker         nargoutchk                  strchr
-      fstat              pathdef                     strvcat
-      full               perl                        subspace
-      fzero              prctile                     symvar
-      genvarname         quantile                    treelayout
-      histc              re_read_readline_init_file  validatestring
+      cellslices         info                        realsqrt
+      cgs                interp1q                    rectint
+      command_line_path  isdebugmode                 regexptranslate
+      contrast           isfloat                     restoredefaultpath
+      convn              isstrprop                   roundb
+      cummin             log1p                       rundemos
+      cummax             lsqnonneg                   runlength
+      datetick           matlabroot                  saveobj
+      display            namelengthmax               spaugment
+      expm1              nargoutchk                  strchr
+      filemarker         pathdef                     strvcat
+      fstat              perl                        subspace
+      full               prctile                     symvar
+      fzero              quantile                    treelayout
+      genvarname         re_read_readline_init_file  validatestring
+      histc
 
  ** Changes to strcat.
 
--- a/aclocal.m4
+++ b/aclocal.m4
@@ -1285,6 +1285,41 @@
    AC_SUBST([FT2_LIBS])])
 dnl end of freetype2.m4
 
+dnl Check whether a math mapper function is available in <cmath>.
+dnl Will define HAVE_CMATH_FUNC if there is a double variant and
+dnl HAVE_CMATH_FUNCF if there is a float variant.
+dnl Currently capable of checking for functions with single 
+dnl argument and returning bool/int/real.
+AC_DEFUN([OCTAVE_CMATH_FUNC],[
+AC_MSG_CHECKING([for std::$1 in <cmath>])
+AC_LANG_PUSH(C++)
+AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+#include <cmath>
+void take_func (bool (*func) (double x));
+void take_func (int (*func) (double x));
+void take_func (double (*func) (double x));
+]],
+[[
+take_func(std::$1);
+]])],
+[AC_MSG_RESULT([yes])
+ AC_DEFINE(HAVE_CMATH_[]AS_TR_CPP($1),1,[Define if <cmath> provides $1])],
+[AC_MSG_RESULT([no])])
+AC_MSG_CHECKING([for std::$1 (float variant) in <cmath>])
+AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+#include <cmath>
+void take_func (bool (*func) (float x));
+void take_func (int (*func) (float x));
+void take_func (float (*func) (float x));
+]],
+[[
+take_func(std::$1);
+]])],
+[AC_MSG_RESULT([yes])
+ AC_DEFINE(HAVE_CMATH_[]AS_TR_CPP($1)F,1,[Define if <cmath> provides float variant of $1])],
+[AC_MSG_RESULT([no])])
+])
+
 dnl Check whether fast signed integer arithmetics using bit tricks
 dnl can be used in oct-inttypes.h. Defines HAVE_FAST_INT_OPS if
 dnl the following conditions hold:
--- a/configure.in
+++ b/configure.in
@@ -1776,6 +1776,12 @@
      [Define if your system has a single-arg prototype for gettimeofday.])])
 fi
 
+dnl Maybe <cmath> defines the IEEE functions we need.
+
+OCTAVE_CMATH_FUNC(isnan)
+OCTAVE_CMATH_FUNC(isinf)
+OCTAVE_CMATH_FUNC(isfinite)
+
 dnl Would like to get rid of this cruft, and just have
 dnl
 dnl   AC_CHECK_FUNCS(finite isnan isinf)
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,8 @@
+2009-03-25  John W. Eaton  <jwe@octave.org>
+
+	* interpreter/munge-texi.cc (process_texi_input_file):
+	Copy leading comment with file name info to output.
+
 2009-03-09  John W. Eaton  <jwe@octave.org>
 
 	* interpreter/Makefile.in (DISTFILES): Use doc-cache instead of
--- a/doc/interpreter/munge-texi.cc
+++ b/doc/interpreter/munge-texi.cc
@@ -245,7 +245,23 @@
 		    {
 		      std::string doc_string = help_text[symbol_name];
 
+		      size_t len = doc_string.length ();
+
 		      int j = 0;
+
+		      // If there is a leading comment with the file
+		      // name, copy it to the output.
+		      if (len > 1
+			  && doc_string[j] == '@'
+			  && doc_string[j+1] == 'c')
+			{
+			  j = 2;
+			  while (doc_string[j++] != '\n')
+			    /* find eol */;
+
+			  os << doc_string.substr (0, j);
+			}
+
 		      while (doc_string[j] == ' ')
 			j++;
 
--- a/liboctave/Array-d.cc
+++ b/liboctave/Array-d.cc
@@ -27,7 +27,7 @@
 
 // Instantiate Arrays of double values.
 
-#include "lo-ieee.h"
+#include "lo-mappers.h"
 #include "Array.h"
 #include "Array.cc"
 #include "oct-locbuf.h"
@@ -40,19 +40,19 @@
 inline bool
 sort_isnan<double> (double x)
 {
-  return lo_ieee_isnan (x);
+  return xisnan (x);
 }
 
 static bool
 nan_ascending_compare (double x, double y)
 {
-  return lo_ieee_isnan (y) ? ! lo_ieee_isnan (x) : x < y;
+  return xisnan (y) ? ! xisnan (x) : x < y;
 }
 
 static bool
 nan_descending_compare (double x, double y)
 {
-  return lo_ieee_isnan (x) ? ! lo_ieee_isnan (y) : x > y;
+  return xisnan (x) ? ! xisnan (y) : x > y;
 }
 
 Array<double>::compare_fcn_type
--- a/liboctave/Array-f.cc
+++ b/liboctave/Array-f.cc
@@ -27,7 +27,7 @@
 
 // Instantiate Arrays of float values.
 
-#include "lo-ieee.h"
+#include "lo-mappers.h"
 #include "Array.h"
 #include "Array.cc"
 #include "oct-locbuf.h"
@@ -40,19 +40,19 @@
 inline bool
 sort_isnan<float> (float x)
 {
-  return lo_ieee_isnan (x);
+  return xisnan (x);
 }
 
 static bool
 nan_ascending_compare (float x, float y)
 {
-  return lo_ieee_isnan (y) ? ! lo_ieee_isnan (x) : x < y;
+  return xisnan (y) ? ! xisnan (x) : x < y;
 }
 
 static bool
 nan_descending_compare (float x, float y)
 {
-  return lo_ieee_isnan (x) ? ! lo_ieee_isnan (y) : x > y;
+  return xisnan (x) ? ! xisnan (y) : x > y;
 }
 
 Array<float>::compare_fcn_type
--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -461,6 +461,7 @@
 class rec_permute_helper
 {
   octave_idx_type *dim, *stride;
+  bool use_blk;
   int top;
 
 public:
@@ -472,37 +473,78 @@
       dim = new octave_idx_type [2*n];
       // A hack to avoid double allocation
       stride = dim + n;
-      top = 0;
 
       // Get cumulative dimensions.
       OCTAVE_LOCAL_BUFFER (octave_idx_type, cdim, n+1);
       cdim[0] = 1;
       for (int i = 1; i < n+1; i++) cdim[i] = cdim[i-1] * dv(i-1);
 
-      int k = 0;
-      // Reduce leading identity
-      for (; k < n && perm(k) == k; k++) ;
-      if (k > 0)
+      // Setup the permuted strides.
+      for (int k = 0; k < n; k++)
         {
-          dim[0] = cdim[k];
-          stride[0] = 1;
+          int kk = perm(k);
+          dim[k] = dv(kk);
+          stride[k] = cdim[kk];
         }
-      else
-        top = -1;
-
-      // Setup the permuted strides.
-      for (; k < n; k++)
+
+      // Reduce contiguous runs.
+      top = 0;
+      for (int k = 1; k < n; k++)
         {
-          ++top;
-          int kk = perm(k);
-          dim[top] = dv(kk);
-          stride[top] = cdim[kk];
+          if (stride[k] == stride[top]*dim[top])
+            dim[top] *= dim[k];
+          else
+            {
+              top++;
+              dim[top] = dim[k];
+              stride[top] = stride[k];
+            }
         }
 
+      // Determine whether we can use block transposes.
+      use_blk = top >= 1 && stride[1] == 1 && stride[0] == dim[1];
+
     }
 
   ~rec_permute_helper (void) { delete [] dim; }
 
+  // Helper method for fast blocked transpose.
+  template <class T>
+  T *blk_trans (const T *src, T *dest, octave_idx_type nr, octave_idx_type nc) const
+    {
+      static const octave_idx_type m = 8;
+      OCTAVE_LOCAL_BUFFER (T, blk, m*m);
+      for (octave_idx_type kr = 0; kr < nr; kr += m)
+        for (octave_idx_type kc = 0; kc < nc; kc += m)
+          {
+            octave_idx_type lr = std::min (m, nr - kr);
+            octave_idx_type lc = std::min (m, nc - kc);
+            if (lr == m && lc == m)
+              {
+                const T *ss = src + kc * nr + kr;
+                for (octave_idx_type j = 0; j < m; j++)
+                  for (octave_idx_type i = 0; i < m; i++)
+                    blk[j*m+i] = ss[j*nr + i];
+                T *dd = dest + kr * nc + kc;
+                for (octave_idx_type j = 0; j < m; j++)
+                  for (octave_idx_type i = 0; i < m; i++)
+                    dd[j*nc+i] = blk[i*m+j];
+              }
+            else
+              {
+                const T *ss = src + kc * nr + kr;
+                for (octave_idx_type j = 0; j < lc; j++)
+                  for (octave_idx_type i = 0; i < lr; i++)
+                    blk[j*m+i] = ss[j*nr + i];
+                T *dd = dest + kr * nc + kc;
+                for (octave_idx_type j = 0; j < lr; j++)
+                  for (octave_idx_type i = 0; i < lc; i++)
+                    dd[j*nc+i] = blk[i*m+j];
+              }
+          }
+
+      return dest + nr*nc;
+    }
 private:
 
   // Recursive N-d generalized transpose
@@ -521,8 +563,9 @@
 
               dest += len;
             }
-
         }
+      else if (use_blk && lev == 1)
+        dest = blk_trans (src, dest, dim[1], dim[0]);
       else
         {
           octave_idx_type step = stride[lev], len = dim[lev];
--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -601,6 +601,26 @@
 
     return result;
   }
+
+  // This is non-breakable map, suitable for fast functions. Efficiency
+  // relies on compiler's ability to inline a function pointer. This seems
+  // to be OK with recent GCC.
+  template <class U>
+  Array<U>
+  fastmap (U (*fcn) (typename ref_param<T>::type)) const
+  {
+    octave_idx_type len = length ();
+
+    const T *m = data ();
+
+    Array<U> result (dims ());
+    U *p = result.fortran_vec ();
+
+    std::transform (m, m + len, p, fcn);
+
+    return result;
+  }
+
 };
 
 #endif
--- a/liboctave/CColVector.cc
+++ b/liboctave/CColVector.cc
@@ -524,9 +524,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       double tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/CDiagMatrix.cc
+++ b/liboctave/CDiagMatrix.cc
@@ -30,6 +30,7 @@
 
 #include "Array-util.h"
 #include "lo-error.h"
+#include "lo-ieee.h"
 #include "mx-base.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -3640,9 +3640,7 @@
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       Complex tmp;
       for (octave_idx_type i = 0; i < nr; i++)
--- a/liboctave/CNDArray.cc
+++ b/liboctave/CNDArray.cc
@@ -754,6 +754,24 @@
                   dims ());
 }
 
+boolNDArray
+ComplexNDArray::isnan (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisnan));
+}
+
+boolNDArray
+ComplexNDArray::isinf (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisinf));
+}
+
+boolNDArray
+ComplexNDArray::isfinite (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xfinite));
+}
+
 ComplexNDArray
 conj (const ComplexNDArray& a)
 {
@@ -898,9 +916,7 @@
 {
   octave_idx_type nel = a.nelem ();
 
-  if (nel < 1 )
-    is.clear (std::ios::badbit);
-  else
+  if (nel > 0)
     {
       Complex tmp;
       for (octave_idx_type i = 0; i < nel; i++)
--- a/liboctave/CNDArray.h
+++ b/liboctave/CNDArray.h
@@ -98,6 +98,9 @@
   ComplexNDArray& insert (const ComplexNDArray& a, const Array<octave_idx_type>& ra_idx);
   
   NDArray abs (void) const;
+  boolNDArray isnan (void) const;
+  boolNDArray isinf (void) const;
+  boolNDArray isfinite (void) const;
 
   friend ComplexNDArray conj (const ComplexNDArray& a);
 
--- a/liboctave/CRowVector.cc
+++ b/liboctave/CRowVector.cc
@@ -434,9 +434,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       Complex tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -646,28 +646,28 @@
   octave_idx_type nz = nnz ();
   SparseComplexMatrix retval (nc, nr, nz);
 
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr + 1);
-  for (octave_idx_type i = 0; i < nr; i++)
-    w[i] = 0;
   for (octave_idx_type i = 0; i < nz; i++)
-    w[ridx(i)]++;
+    retval.xcidx (ridx (i) + 1)++;
+  // retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1)
   nz = 0;
-  for (octave_idx_type i = 0; i < nr; i++)
+  for (octave_idx_type i = 1; i <= nr; i++)
     {
-      retval.xcidx(i) = nz;
-      nz += w[i];
-      w[i] = retval.xcidx(i);
+      const octave_idx_type tmp = retval.xcidx (i);
+      retval.xcidx (i) = nz;
+      nz += tmp;
     }
-  retval.xcidx(nr) = nz;
-  w[nr] = nz;
+  // retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1)
 
   for (octave_idx_type j = 0; j < nc; j++)
     for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
       {
-	octave_idx_type q = w [ridx(k)]++;
+	octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
 	retval.xridx (q) = j;
 	retval.xdata (q) = conj (data (k));
       }
+  assert (nnz () == retval.xcidx (nr));
+  // retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1)
+  // and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets
 
   return retval;
 }
@@ -7470,9 +7470,7 @@
   octave_idx_type nc = a.cols ();
   octave_idx_type nz = a.nzmax ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       octave_idx_type itmp, jtmp, jold = 0;
       Complex tmp;
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,73 @@
+2009-03-25  John W. Eaton  <jwe@octave.org>
+
+	* Makefile.in (MATRIX_INC): Add Sparse-diag-op-defs.h and
+	Sparse-perm-op-defs.h to the list.
+
+2009-03-25  Jaroslav Hajek  <highegg@gmail.com>
+	
+	* oct-inttypes.cc (INT_DOUBLE_BINOP_DECL (*, uint64),
+	INT_DOUBLE_BINOP_DECL (*, int64)): x -> y where appropriate.
+
+2009-03-25  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Array.cc (rec_permute_helper::use_blk): New field.
+	(rec_permute_helper::blk_trans): New method.
+	(rec_permute_helper::rec_permute_helper): Use smart reductions,
+	detect possibility of using blocked transpose.
+	(rec_permute_helper::do_permute): Use blocked transpose if possible.
+
+2009-03-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* idx-vector.cc (convert_index(double,...)): Simplify.
+
+2009-03-21  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Array-d.cc: lo_ieee_isnan -> xisnan.
+	* Array-f.cc: Ditto.
+	* oct-inttypes.cc: Ditto.
+	* oct-inttypes.h: Ditto.
+	* CDiagMatrix.cc: Add missing include.
+	* fCDiagMatrix.cc: Ditto.
+
+2009-03-20  Jaroslav Hajek  <highegg@gmail.com>
+	
+	* CColVector.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc,
+	boolSparse.cc, dColVector.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc,
+	dSparse.cc, fCColVector.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc,
+	fColVector.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, intNDArray.cc:
+	Allow empty arrays in stream input operators.
+
+2009-03-20  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Array.h (Array<T>::fastmap): New method.
+	* dNDArray.cc (NDArray::isnan, NDArray::isinf, NDArray::isfinite):
+	New methods.
+	* dNDArray.h: Declare them.
+	* fNDArray.cc (FloatNDArray::isnan, FloatNDArray::isinf,
+	FloatNDArray::isfinite): New methods.
+	* fNDArray.h: Declare them.
+	* CNDArray.cc (ComplexNDArray::isnan, ComplexNDArray::isinf,
+	ComplexNDArray::isfinite): New methods.
+	* CNDArray.h: Declare them.
+	* fCNDArray.cc (FloatComplexNDArray::isnan, FloatComplexNDArray::isinf,
+	FloatComplexNDArray::isfinite): New methods.
+	* fCNDArray.h: Declare them.
+	* lo-mappers.h (xisnan, xisinf, xfinite): If possible, use definitions
+	from <cmath>.
+
+2009-03-18  Jaroslav Hajek <highegg@gmail.com>
+
+	* oct-norm.cc (get_eps): Remove that hack.
+	(higham): Use std::numeric_limits instead.
+	Include OCTAVE_QUIT.
+
+2009-03-16  Jason Riedy  <jason@acm.org>
+
+	* Sparse.cc (transpose): Eliminate the workspace by computing in
+	retval.xcidx.
+	* CSparse.cc (hermitian): Eliminate the workspace by computing in
+	retval.xcidx.
+
 2009-03-14  Jaroslav Hajek  <highegg@gmail.com>
 
 	* mx-op-decl.h (NDS_BOOL_OP_DECLS, SND_BOOL_OP_DECLS, NDND_BOOL_OP_DECLS): Support compound binary ops.
--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -53,9 +53,9 @@
 	dbleGEPBAL.h dbleHESS.h dbleLU.h dbleQR.h dbleQRP.h dbleSCHUR.h \
 	dbleSVD.h boolSparse.h CSparse.h dSparse.h MSparse-defs.h MSparse.h \
 	Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \
-	sparse-base-chol.h SparseCmplxCHOL.h \
-	SparsedbleCHOL.h SparseCmplxQR.h SparseQR.h Sparse-op-defs.h \
-	MatrixType.h PermMatrix.h \
+	sparse-base-chol.h SparseCmplxCHOL.h SparsedbleCHOL.h \
+	SparseCmplxQR.h SparseQR.h Sparse-op-defs.h Sparse-diag-op-defs.h \
+	Sparse-perm-op-defs.h MatrixType.h PermMatrix.h \
 	int8NDArray.h uint8NDArray.h int16NDArray.h uint16NDArray.h \
 	int32NDArray.h uint32NDArray.h int64NDArray.h uint64NDArray.h \
 	intNDArray.h \
--- a/liboctave/Sparse.cc
+++ b/liboctave/Sparse.cc
@@ -1062,28 +1062,28 @@
   octave_idx_type nz = nnz ();
   Sparse<T> retval (nc, nr, nz);
 
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr + 1);
-  for (octave_idx_type i = 0; i < nr; i++)
-    w[i] = 0;
   for (octave_idx_type i = 0; i < nz; i++)
-    w[ridx(i)]++;
+    retval.xcidx (ridx (i) + 1)++;
+  // retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1)
   nz = 0;
-  for (octave_idx_type i = 0; i < nr; i++)
+  for (octave_idx_type i = 1; i <= nr; i++)
     {
-      retval.xcidx(i) = nz;
-      nz += w[i];
-      w[i] = retval.xcidx(i);
+      const octave_idx_type tmp = retval.xcidx (i);
+      retval.xcidx (i) = nz;
+      nz += tmp;
     }
-  retval.xcidx(nr) = nz;
-  w[nr] = nz;
+  // retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1)
 
   for (octave_idx_type j = 0; j < nc; j++)
     for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
       {
-	octave_idx_type q = w [ridx(k)]++;
+	octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
 	retval.xridx (q) = j;
 	retval.xdata (q) = data (k);
       }
+  assert (nnz () == retval.xcidx (nr));
+  // retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1)
+  // and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets
 
   return retval;
 }
--- a/liboctave/boolSparse.cc
+++ b/liboctave/boolSparse.cc
@@ -184,9 +184,7 @@
   octave_idx_type nc = a.cols ();
   octave_idx_type nz = a.nzmax ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       octave_idx_type itmp, jtmp, jold = 0;
       bool tmp;
--- a/liboctave/dColVector.cc
+++ b/liboctave/dColVector.cc
@@ -321,9 +321,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       double tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -3058,9 +3058,7 @@
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       double tmp;
       for (octave_idx_type i = 0; i < nr; i++)
--- a/liboctave/dNDArray.cc
+++ b/liboctave/dNDArray.cc
@@ -874,6 +874,24 @@
                   dims ());
 }
 
+boolNDArray
+NDArray::isnan (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisnan));
+}
+
+boolNDArray
+NDArray::isinf (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisinf));
+}
+
+boolNDArray
+NDArray::isfinite (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xfinite));
+}
+
 Matrix
 NDArray::matrix_value (void) const
 {
@@ -947,9 +965,7 @@
 {
   octave_idx_type nel = a.nelem ();
 
-  if (nel < 1 )
-    is.clear (std::ios::badbit);
-  else
+  if (nel > 0)
     {
       double tmp;
       for (octave_idx_type i = 0; i < nel; i++)
--- a/liboctave/dNDArray.h
+++ b/liboctave/dNDArray.h
@@ -109,6 +109,9 @@
   NDArray& insert (const NDArray& a, const Array<octave_idx_type>& ra_idx);
 
   NDArray abs (void) const;
+  boolNDArray isnan (void) const;
+  boolNDArray isinf (void) const;
+  boolNDArray isfinite (void) const;
 
   ComplexNDArray fourier (int dim = 1) const;
   ComplexNDArray ifourier (int dim = 1) const;
--- a/liboctave/dRowVector.cc
+++ b/liboctave/dRowVector.cc
@@ -291,9 +291,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       double tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -7588,9 +7588,7 @@
   octave_idx_type nc = a.cols ();
   octave_idx_type nz = a.nzmax ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       octave_idx_type itmp, jtmp, jold = 0;
       double tmp;
--- a/liboctave/fCColVector.cc
+++ b/liboctave/fCColVector.cc
@@ -524,9 +524,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       float tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/fCDiagMatrix.cc
+++ b/liboctave/fCDiagMatrix.cc
@@ -30,6 +30,7 @@
 
 #include "Array-util.h"
 #include "lo-error.h"
+#include "lo-ieee.h"
 #include "mx-base.h"
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
--- a/liboctave/fCMatrix.cc
+++ b/liboctave/fCMatrix.cc
@@ -3633,9 +3633,7 @@
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       FloatComplex tmp;
       for (octave_idx_type i = 0; i < nr; i++)
--- a/liboctave/fCNDArray.cc
+++ b/liboctave/fCNDArray.cc
@@ -749,6 +749,24 @@
                        dims ());
 }
 
+boolNDArray
+FloatComplexNDArray::isnan (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisnan));
+}
+
+boolNDArray
+FloatComplexNDArray::isinf (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisinf));
+}
+
+boolNDArray
+FloatComplexNDArray::isfinite (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xfinite));
+}
+
 FloatComplexNDArray
 conj (const FloatComplexNDArray& a)
 {
@@ -893,9 +911,7 @@
 {
   octave_idx_type nel = a.nelem ();
 
-  if (nel < 1 )
-    is.clear (std::ios::badbit);
-  else
+  if (nel > 0)
     {
       FloatComplex tmp;
       for (octave_idx_type i = 0; i < nel; i++)
--- a/liboctave/fCNDArray.h
+++ b/liboctave/fCNDArray.h
@@ -98,6 +98,9 @@
   FloatComplexNDArray& insert (const FloatComplexNDArray& a, const Array<octave_idx_type>& ra_idx);
   
   FloatNDArray abs (void) const;
+  boolNDArray isnan (void) const;
+  boolNDArray isinf (void) const;
+  boolNDArray isfinite (void) const;
 
   friend FloatComplexNDArray conj (const FloatComplexNDArray& a);
 
--- a/liboctave/fCRowVector.cc
+++ b/liboctave/fCRowVector.cc
@@ -434,9 +434,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       FloatComplex tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/fColVector.cc
+++ b/liboctave/fColVector.cc
@@ -321,9 +321,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       float tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -3057,9 +3057,7 @@
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
 
-  if (nr < 1 || nc < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (nr > 0 && nc > 0)
     {
       float tmp;
       for (octave_idx_type i = 0; i < nr; i++)
--- a/liboctave/fNDArray.cc
+++ b/liboctave/fNDArray.cc
@@ -829,6 +829,24 @@
                        dims ());
 }
 
+boolNDArray
+FloatNDArray::isnan (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisnan));
+}
+
+boolNDArray
+FloatNDArray::isinf (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xisinf));
+}
+
+boolNDArray
+FloatNDArray::isfinite (void) const
+{
+  return ArrayN<bool> (fastmap<bool> (xfinite));
+}
+
 FloatMatrix
 FloatNDArray::matrix_value (void) const
 {
@@ -902,9 +920,7 @@
 {
   octave_idx_type nel = a.nelem ();
 
-  if (nel < 1 )
-    is.clear (std::ios::badbit);
-  else
+  if (nel > 0)
     {
       float tmp;
       for (octave_idx_type i = 0; i < nel; i++)
--- a/liboctave/fNDArray.h
+++ b/liboctave/fNDArray.h
@@ -106,6 +106,9 @@
   FloatNDArray& insert (const FloatNDArray& a, const Array<octave_idx_type>& ra_idx);
 
   FloatNDArray abs (void) const;
+  boolNDArray isnan (void) const;
+  boolNDArray isinf (void) const;
+  boolNDArray isfinite (void) const;
 
   FloatComplexNDArray fourier (int dim = 1) const;
   FloatComplexNDArray ifourier (int dim = 1) const;
--- a/liboctave/fRowVector.cc
+++ b/liboctave/fRowVector.cc
@@ -291,9 +291,7 @@
 {
   octave_idx_type len = a.length();
 
-  if (len < 1)
-    is.clear (std::ios::badbit);
-  else
+  if (len > 0)
     {
       float tmp;
       for (octave_idx_type i = 0; i < len; i++)
--- a/liboctave/idx-vector.cc
+++ b/liboctave/idx-vector.cc
@@ -182,18 +182,9 @@
 inline octave_idx_type
 convert_index (double x, bool& conv_error, octave_idx_type& ext)
 {
-  octave_idx_type i;
-  if (xisnan (x) || xisinf (x))
-    {
-      i = 0;
-      conv_error = true;
-    }
-  else
-    {
-      i = static_cast<octave_idx_type> (x);
-      if (static_cast<double> (i) != x)
-        conv_error = true;
-    }
+  octave_idx_type i = static_cast<octave_idx_type> (x);
+  if (static_cast<double> (i) != x)
+    conv_error = true;
 
   return convert_index (i, conv_error, ext);
 }
--- a/liboctave/intNDArray.cc
+++ b/liboctave/intNDArray.cc
@@ -146,9 +146,7 @@
 {
   octave_idx_type nel = a.nelem ();
 
-  if (nel < 1 )
-    is.clear (std::ios::badbit);
-  else
+  if (nel > 0)
     {
       T tmp;
 
--- a/liboctave/lo-mappers.cc
+++ b/liboctave/lo-mappers.cc
@@ -190,23 +190,29 @@
 
 // double -> bool mappers.
 
+#if ! defined(HAVE_CMATH_ISNAN)
 bool
 xisnan (double x)
 {
   return lo_ieee_isnan (x);
 }
+#endif
 
+#if ! defined(HAVE_CMATH_ISFINITE)
 bool
 xfinite (double x)
 {
   return lo_ieee_finite (x);
 }
+#endif
 
+#if ! defined(HAVE_CMATH_ISINF)
 bool
 xisinf (double x)
 {
   return lo_ieee_isinf (x);
 }
+#endif
 
 bool
 octave_is_NA (double x)
@@ -321,28 +327,6 @@
 // complex -> bool mappers.
 
 bool
-xisnan (const Complex& x)
-{
-  return (xisnan (real (x)) || xisnan (imag (x)));
-}
-
-bool
-xfinite (const Complex& x)
-{
-  double rx = real (x);
-  double ix = imag (x);
-
-  return (xfinite (rx) && ! xisnan (rx)
-	  && xfinite (ix) && ! xisnan (ix));
-}
-
-bool
-xisinf (const Complex& x)
-{
-  return (xisinf (real (x)) || xisinf (imag (x)));
-}
-
-bool
 octave_is_NA (const Complex& x)
 {
   return (octave_is_NA (real (x)) || octave_is_NA (imag (x)));
@@ -524,23 +508,29 @@
 
 // float -> bool mappers.
 
+#if ! defined(HAVE_CMATH_ISNANF)
 bool
 xisnan (float x)
 {
   return lo_ieee_isnan (x);
 }
+#endif
 
+#if ! defined(HAVE_CMATH_ISFINITEF)
 bool
 xfinite (float x)
 {
   return lo_ieee_finite (x);
 }
+#endif
 
+#if ! defined(HAVE_CMATH_ISINFF)
 bool
 xisinf (float x)
 {
   return lo_ieee_isinf (x);
 }
+#endif
 
 bool
 octave_is_NA (float x)
@@ -655,28 +645,6 @@
 // complex -> bool mappers.
 
 bool
-xisnan (const FloatComplex& x)
-{
-  return (xisnan (real (x)) || xisnan (imag (x)));
-}
-
-bool
-xfinite (const FloatComplex& x)
-{
-  float rx = real (x);
-  float ix = imag (x);
-
-  return (xfinite (rx) && ! xisnan (rx)
-	  && xfinite (ix) && ! xisnan (ix));
-}
-
-bool
-xisinf (const FloatComplex& x)
-{
-  return (xisinf (real (x)) || xisinf (imag (x)));
-}
-
-bool
 octave_is_NA (const FloatComplex& x)
 {
   return (octave_is_NA (real (x)) || octave_is_NA (imag (x)));
--- a/liboctave/lo-mappers.h
+++ b/liboctave/lo-mappers.h
@@ -25,6 +25,7 @@
 #define octave_liboctave_mappers_h 1
 
 #include "oct-cmplx.h"
+#include "lo-math.h"
 
 // Double Precision 
 extern OCTAVE_API double arg (double x);
@@ -46,9 +47,24 @@
 inline bool xisnan (bool) { return false; }
 inline bool xisnan (char) { return false; }
 
+#if defined (HAVE_CMATH_ISNAN)
+inline bool xisnan (double x)
+{ return std::isnan (x); }
+#else
 extern OCTAVE_API bool xisnan (double x);
+#endif
+#if defined (HAVE_CMATH_ISFINITE)
+inline bool xfinite (double x)
+{ return std::isfinite (x); }
+#else
 extern OCTAVE_API bool xfinite (double x);
+#endif
+#if defined (HAVE_CMATH_ISINF)
+inline bool xisinf (double x)
+{ return std::isinf (x); }
+#else
 extern OCTAVE_API bool xisinf (double x);
+#endif
 
 extern OCTAVE_API bool octave_is_NA (double x);
 extern OCTAVE_API bool octave_is_NaN_or_NA (double x) GCC_ATTR_DEPRECATED;
@@ -70,9 +86,15 @@
 extern OCTAVE_API Complex xroundb (const Complex& x);
 extern OCTAVE_API Complex signum (const Complex& x);
 
-extern OCTAVE_API bool xisnan (const Complex& x);
-extern OCTAVE_API bool xfinite (const Complex& x);
-extern OCTAVE_API bool xisinf (const Complex& x);
+inline bool
+xisnan (const Complex& x)
+{ return (xisnan (real (x)) || xisnan (imag (x))); }
+inline bool
+xfinite (const Complex& x)
+{ return (xfinite (real (x)) && xfinite (imag (x))); }
+inline bool
+xisinf (const Complex& x)
+{ return (xisinf (real (x)) || xisinf (imag (x))); }
 
 extern OCTAVE_API bool octave_is_NA (const Complex& x);
 extern OCTAVE_API bool octave_is_NaN_or_NA (const Complex& x);
@@ -96,9 +118,25 @@
 extern OCTAVE_API FloatComplex xlog2 (const FloatComplex& x, int& exp);
 extern OCTAVE_API float xexp2 (float x);
 
+#if defined (HAVE_CMATH_ISNANF)
+inline bool xisnan (float x)
+{ return std::isnan (x); }
+#else
 extern OCTAVE_API bool xisnan (float x);
+#endif
+#if defined (HAVE_CMATH_ISFINITEF)
+inline bool xfinite (float x)
+{ return std::isfinite (x); }
+#else
 extern OCTAVE_API bool xfinite (float x);
+#endif
+#if defined (HAVE_CMATH_ISINFF)
+inline bool xisinf (float x)
+{ return std::isinf (x); }
+#else
 extern OCTAVE_API bool xisinf (float x);
+#endif
+
 
 extern OCTAVE_API bool octave_is_NA (float x);
 extern OCTAVE_API bool octave_is_NaN_or_NA (float x) GCC_ATTR_DEPRECATED;
@@ -120,9 +158,15 @@
 extern OCTAVE_API FloatComplex xroundb (const FloatComplex& x);
 extern OCTAVE_API FloatComplex signum (const FloatComplex& x);
 
-extern OCTAVE_API bool xisnan (const FloatComplex& x);
-extern OCTAVE_API bool xfinite (const FloatComplex& x);
-extern OCTAVE_API bool xisinf (const FloatComplex& x);
+inline bool
+xisnan (const FloatComplex& x)
+{ return (xisnan (real (x)) || xisnan (imag (x))); }
+inline bool
+xfinite (const FloatComplex& x)
+{ return (xfinite (real (x)) && xfinite (imag (x))); }
+inline bool
+xisinf (const FloatComplex& x)
+{ return (xisinf (real (x)) || xisinf (imag (x))); }
 
 extern OCTAVE_API bool octave_is_NA (const FloatComplex& x);
 extern OCTAVE_API bool octave_is_NaN_or_NA (const FloatComplex& x);
--- a/liboctave/mx-op-defs.h
+++ b/liboctave/mx-op-defs.h
@@ -1120,6 +1120,13 @@
 #define MPM_BIN_OPS(R, M, PM) \
   MPM_MULTIPLY_OP(M, PM);
 
+#define NDND_MAPPER_BODY(R, NAME) \
+  R retval (dims ()); \
+  octave_idx_type n = numel (); \
+  for (octave_idx_type i = 0; i < n; i++) \
+    retval.xelem (i) = NAME (elem (i)); \
+  return retval;
+
 #endif
 
 
--- a/liboctave/oct-inttypes.cc
+++ b/liboctave/oct-inttypes.cc
@@ -407,7 +407,7 @@
     {
       return x / octave_uint64 (static_cast<uint64_t> (2));
     }
-  else if (y < 0 || lo_ieee_isnan (x) || lo_ieee_isinf (x))
+  else if (y < 0 || xisnan (y) || xisinf (y))
     {
       return octave_uint64 (x.value () * y); 
     }
@@ -442,7 +442,7 @@
     {
       return x / octave_int64 (static_cast<uint64_t> (4*y));
     }
-  else if (lo_ieee_isnan (x) || lo_ieee_isinf (x))
+  else if (xisnan (y) || xisinf (y))
     {
       return octave_int64 (x.value () * y); 
     }
--- a/liboctave/oct-inttypes.h
+++ b/liboctave/oct-inttypes.h
@@ -33,11 +33,11 @@
 #include "lo-traits.h"
 #include "lo-math.h"
 #include "oct-types.h"
-#include "lo-ieee.h"
 #include "lo-mappers.h"
 
 #ifdef OCTAVE_INT_USE_LONG_DOUBLE
 inline long double xround (long double x) { return roundl (x); }
+inline long double xisnan (long double x) { return xisnan (static_cast<double> (x)); }
 #endif
 
 // Undefine min/max if needed (this may happen under Windows)
@@ -289,7 +289,7 @@
       // Compute proper thresholds.
       static const S thmin = compute_threshold (static_cast<S> (min_val ()), min_val ());
       static const S thmax = compute_threshold (static_cast<S> (max_val ()), max_val ());
-      if (lo_ieee_isnan (value))
+      if (xisnan (value))
         {
           fnan = true;
           return static_cast<T> (0);
--- a/liboctave/oct-norm.cc
+++ b/liboctave/oct-norm.cc
@@ -413,6 +413,7 @@
   RR lambda = 0, mu = 0;
   for (octave_idx_type k = 0; k < m.columns (); k++)
     {
+      OCTAVE_QUIT;
       VectorT col (m.column (k));
       if (k > 0)
         higham_subp (y, col, 4*k, p, lambda, mu);
@@ -430,6 +431,7 @@
   int iter = 0;
   while (iter < maxiter)
     {
+      OCTAVE_QUIT;
       y = m*x;
       gamma1 = gamma;
       gamma = vector_norm (y, p);
@@ -449,9 +451,6 @@
   return gamma;
 }
 
-// TODO: is there a better way?
-inline float get_eps (float) { return FLT_EPSILON; }
-inline double get_eps (double) { return DBL_EPSILON; }
 // derive column vector and SVD types 
 
 static const char *p_less1_gripe = "xnorm: p must be at least 1";
@@ -477,7 +476,8 @@
   else if (p > 1)
     {
       VectorT x;
-      res = higham (m, p, std::sqrt (get_eps (p)), max_norm_iter, x);
+      const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
+      res = higham (m, p, sqrteps, max_norm_iter, x);
     }
   else
     (*current_liboctave_error_handler) (p_less1_gripe); 
@@ -497,7 +497,8 @@
   else if (p > 1)
     {
       VectorT x;
-      res = higham (m, p, std::sqrt (get_eps (p)), max_norm_iter, x);
+      const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
+      res = higham (m, p, sqrteps, max_norm_iter, x);
     }
   else
     (*current_liboctave_error_handler) (p_less1_gripe); 
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,52 @@
+2009-03-25  Kai Habel  <kai.habel@gmx.de>
+
+	* general/gradient.m: Fix calculation for more than two
+	dimensions.  Change interpretation of vector arguments from
+	spacing to coordinates.  New tests.
+
+2009-03-25  John W. Eaton  <jwe@octave.org>
+
+	* mkdoc: Pass full file name to gethelp.
+	* gethelp.cc (main): Handle second argument.  Write comment with
+	full file name to output.
+
+2009-03-24  Ben Abbott <bpabbott@mac.com>
+
+	* plot/gnuplot_drawnow.m: When printing, pass scalar plot_stream
+	to __gnuplot_draw_figure__, and close all plot streams when done.
+
+2009-03-24  John W. Eaton  <jwe@octave.org>
+
+	* general/isa.m: Handle parent classes.
+
+2009-03-23  Ben Abbott <bpabbott@mac.com>
+
+	* plot/gnuplot_drawnow.m: Check that gnuplot has internal variable
+	"GPVAL_TERMINALS".
+	* plot/__gnuplot_has_feature__.m: Add "variable_GPVAL_TERMINALS".
+
+2009-03-21  Ben Abbott <bpabbott@mac.com>
+
+	* plot/gnuplot_drawnow.m: Verify the gnuplot terminal is supported.
+	* plot/__gnuplot_get_var__.m: Add function to get gnuplot variables.
+	* plot/print.m: Restore the behavior for option -S<num>,<num>.
+
+2009-03-19  Jaroslav Hajek <highegg@gmail.com>
+
+	* optimization/fsolve.m (guarded_eval): Simplify & fix missing
+	semicolon.
+
+2009-03-17  Jaroslav Hajek  <highegg@gmail.com>
+
+	* optimization/__fdjac__.m: Pass in fvec to save one evaluation.
+	* optimization/fsolve.m: Avoid redundant reevaluation when using
+	FD jacobians. Document how it can be done with user jacobians.  Make
+	first iteration special and call outputfcn after it. Skip updates
+	unless two successful iterations have occured.
+	* optimization/__dogleg__.m: Add missing alpha in the zero-gradient
+	case.
+	* optimization/fsolve.m: Remove autodg (not used), simplify.
+
 2009-03-14  Jaroslav Hajek  <highegg@gmail.com>
 
 	* statistics/base/var.m: a -> x.
--- a/scripts/general/gradient.m
+++ b/scripts/general/gradient.m
@@ -17,36 +17,36 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} gradient (@var{m})
-## @deftypefnx {Function File} {[@var{x}, @var{y}, @dots{}] =} gradient (@var{m})
+## @deftypefn {Function File} {@var{dx} =} gradient (@var{m})
+## @deftypefnx {Function File} {[@var{dx}, @var{dy}, @var{dz}, @dots{}] =} gradient (@var{m})
 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{s})
-## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{dx}, @var{dy}, @dots{})
+## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{x}, @var{y}, @var{z}, @dots{})
 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0})
 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{s})
-## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{dx}, @var{dy}, @dots{})
+## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{x}, @var{y}, @dots{})
 ##
 ## Calculate the gradient of sampled data, or of a function.  If @var{m}
 ## is a vector, calculate the one dimensional gradient of @var{m}. If
-## @var{m} is a matrix the gradient is calculated for each row.
+## @var{m} is a matrix the gradient is calculated for each dimension.
 ##
-## @code{[@var{x}, @var{y}] = gradient (@var{m})} calculates the one
-## dimensional gradient for each direction if @var{m} if @var{m} is a
+## @code{[@var{dx}, @var{dy}] = gradient (@var{m})} calculates the one
+## dimensional gradient for @var{x} and @var{y} direction if @var{m} is a
 ## matrix. Additional return arguments can be use for multi-dimensional
 ## matrices.
 ##
-## Spacing values between two points can be provided by the
-## @var{dx}, @var{dy} or @var{h} parameters. If @var{h} is supplied it
-## is assumed to be the spacing in all directions. Otherwise, separate
-## values of the spacing can be supplied by the @var{dx}, etc variables.
-## A scalar value specifies an equidistant spacing, while a vector value
-## can be used to specify a variable spacing. The length must match
-## their respective dimension of @var{m}.
+## A constant spacing between two points can be provided by the
+## @var{s} parameter. If @var{s} is a scalar, it is assumed to be the spacing
+## for all dimensions. 
+## Otherwise, separate values of the spacing can be supplied by
+## the @var{x}, @dots{} arguments. Scalar values specify an equidistant spacing.
+## Vector values for the @var{x}, @dots{} arguments specify the coordinate for that
+## dimension. The length must match their respective dimension of @var{m}.
 ## 
 ## At boundary points a linear extrapolation is applied. Interior points
 ## are calculated with the first approximation of the numerical gradient
 ##
 ## @example
-## y'(i) = 1/(x(i+1)-x(i-1)) *(y(i-1)-y(i+1)).
+## y'(i) = 1/(x(i+1)-x(i-1)) * (y(i-1)-y(i+1)).
 ## @end example
 ## 
 ## If the first argument @var{f} is a function handle, the gradient of the
@@ -83,81 +83,85 @@
 function varargout = matrix_gradient (m, varargin)
   transposed = false;
   if (isvector (m))
-    ## make a column vector.
+    ## make a row vector.
     transposed = (size (m, 2) == 1);
     m = m(:)';
   endif
 
   nd = ndims (m);
   sz = size (m);
+  if (length(sz) > 1)
+    tmp = sz(1); sz(1) = sz(2); sz(2) = tmp;
+  endif
+
   if (nargin > 2 && nargin != nd + 1)
     print_usage ()
   endif
-    
+  
+  ## cell d stores a spacing vector for each dimension
   d = cell (1, nd);
   if (nargin == 1)
+    ## no spacing given - assume 1.0 for all dimensions
     for i = 1:nd
-      d{i} = ones (sz(i), 1);
+      d{i} = ones (sz(i) - 1, 1);
     endfor
   elseif (nargin == 2)
     if (isscalar (varargin{1}))
+      ## single scalar value for all dimensions
       for i = 1:nd
-	d{i} = varargin{1} * ones (sz(i), 1);
+        d{i} = varargin{1} * ones (sz(i) - 1, 1);
       endfor
     else
-      for i = 1:nd
-	d{i} = varargin{1};
-      endfor
+      ## vector for one-dimensional derivative
+      d{1} = diff (varargin{1}(:));
     endif
   else
+    ## have spacing value for each dimension
+    if (length(varargin) != nd)
+      error ("dimensions and number of spacing values do not match.");
+    end
     for i = 1:nd
       if (isscalar (varargin{i}))
-	## Why the hell did Matlab decide to swap these two values?
-	if (i == 1)
-	  d{2} = varargin{1} * ones (sz(2), 1);
-	elseif (i == 2)
-	  d{1} = varargin{2} * ones (sz(1), 1);
-	else
-	  d{i} = varargin{i} * ones (sz(i), 1);
-	endif
+        d{i} = varargin{i} * ones (sz(i) - 1, 1);
       else
-	## Why the hell did Matlab decide to swap these two values?
-	if (i == 1)
-	  d{2} = varargin{1};
-	elseif (i == 2)
-	  d{1} = varargin{2};
-	else
-	  d{i} = varargin{i};
-	endif
+        d{i} = diff (varargin{i}(:));
       endif
     endfor
   endif
 
-  for i = 1:max (2, min (nd, nargout))
-    mr = sz(i);
-    mc = prod ([sz(1:i-1), sz(i+1:nd)]);
+  m = shiftdim (m, 1);
+  for i = 1:min (nd, nargout)
+    mr = rows (m);
+    mc = numel (m) / mr;
     Y = zeros (size (m), class (m));
-
+	
     if (mr > 1)
       ## Top and bottom boundary.
-      Y(1,:) = diff (m(1:2,:)) / d{i}(1);
-      Y(mr,:) = diff (m(mr-1:mr,:)) / d{i}(mr-1);
+      Y(1,:) = diff (m(1:2, :)) / d{i}(1);
+      Y(mr,:) = diff (m(mr-1:mr, :) / d{i}(mr - 1));
     endif
 
     if (mr > 2)
       ## Interior points.
       Y(2:mr-1,:) = ((m(3:mr,:) - m(1:mr-2,:))
-		     ./ kron (d{i}(1:mr-2) + d{i}(2:mr-1), ones (1, mc)));
+          ./ kron (d{i}(1:mr-2) + d{i}(2:mr-1), ones (1, mc)));
     endif
-    varargout{i} = ipermute (Y, [i:nd,1:i-1]);
-    m = permute (m, [2:nd,1]);
+
+    ## turn multi-dimensional matrix in a way, that gradient
+    ## along x-direction is calculated first then y, z, ...
+
+    if (i == 1)
+      varargout{i} = shiftdim (Y, nd - 1);
+      m = shiftdim (m, nd - 1);
+    elseif (i == 2)
+      varargout{i} = Y;
+      m = shiftdim (m, 2);
+    else
+      varargout{i} = shiftdim (Y, nd - i + 1);
+      m = shiftdim (m, 1);
+    endif
   endfor
 
-  ## Why the hell did Matlab decide to swap these two values?
-  tmp = varargout{1};
-  varargout{1} = varargout{2};
-  varargout{2} = tmp;
-
   if (transposed)
     varargout{1} = varargout{1}.';
   endif
@@ -200,7 +204,7 @@
   
   ## Calculate the gradient
   p0 = mat2cell (p0, num_points, ones (1, dim));
-  varargout = cell (1,dim);
+  varargout = cell (1, dim);
   for d = 1:dim
     s = delta (d);
     df_dx = (f (p0{1:d-1}, p0{d}+s, p0{d+1:end})
@@ -216,7 +220,40 @@
 %!test
 %! data = [1, 2, 4, 2];
 %! dx = gradient (data);
+%! dx2 = gradient (data, 0.25);
+%! dx3 = gradient (data, [0.25, 0.5, 1, 3]);
 %! assert (dx, [1, 3/2, 0, -2]);
+%! assert (dx2, [4, 6, 0, -8]);
+%! assert (dx3, [4, 4, 0, -1]);
+%! assert (size_equal(data, dx));
+
+%!test
+%! [Y,X,Z,U] = ndgrid (2:2:8,1:5,4:4:12,3:5:30);
+%! [dX,dY,dZ,dU] = gradient (X);
+%! assert (all(dX(:)==1));
+%! assert (all(dY(:)==0));
+%! assert (all(dZ(:)==0));
+%! assert (all(dU(:)==0));
+%! [dX,dY,dZ,dU] = gradient (Y);
+%! assert (all(dX(:)==0));
+%! assert (all(dY(:)==2));
+%! assert (all(dZ(:)==0));
+%! assert (all(dU(:)==0));
+%! [dX,dY,dZ,dU] = gradient (Z);
+%! assert (all(dX(:)==0));
+%! assert (all(dY(:)==0));
+%! assert (all(dZ(:)==4));
+%! assert (all(dU(:)==0));
+%! [dX,dY,dZ,dU] = gradient (U);
+%! assert (all(dX(:)==0));
+%! assert (all(dY(:)==0));
+%! assert (all(dZ(:)==0));
+%! assert (all(dU(:)==5));
+%! assert (size_equal(dX, dY, dZ, dU, X, Y, Z, U));
+%! [dX,dY,dZ,dU] = gradient (U, 5.0);
+%! assert (all(dU(:)==1));
+%! [dX,dY,dZ,dU] = gradient (U, 1.0, 2.0, 3.0, 2.5);
+%! assert (all(dU(:)==2));
 
 %!test
 %! x = 0:10;
--- a/scripts/general/isa.m
+++ b/scripts/general/isa.m
@@ -41,7 +41,14 @@
   elseif (strcmp (cname, "numeric"))
     retval = any (strcmp (class (x), fnum_classes));
   else
-    retval = strcmp (class (x), cname);
+    class_of_x = class (x);
+    retval = strcmp (class_of_x, cname);
+    if (! retval && isobject (x))
+      parent_classes_of_x = __parent_classes__ (x);
+      if (! isempty (parent_classes_of_x))
+	retval = any (strcmp (parent_classes_of_x, cname));
+      endif
+    endif
   endif
 
 endfunction
--- a/scripts/gethelp.cc
+++ b/scripts/gethelp.cc
@@ -135,20 +135,26 @@
 main (int argc, char **argv)
 {
   std::string name;
+  std::string file_name;
 
-  if (argc != 2)
+  if (argc != 3)
     {
-      std::cerr << "usage: gethelp name\n";
+      std::cerr << "usage: gethelp name file-name\n";
       return 1;
     }
   else
-    name = argv[1];
+    {
+      name = argv[1];
+      file_name = argv[2];
+    }
 
   std::string help_text = extract_help_text ();  
 
   if (! help_text.empty ())
     {
-      std::cout << "" << name << "\n" << help_text;
+      std::cout << "" << name << "\n"
+		<< "@c " << file_name << "\n"
+		<< help_text;
 
       if (help_text[help_text.length () - 1] != '\n')
 	std::cout << "\n";
--- a/scripts/mkdoc
+++ b/scripts/mkdoc
@@ -42,15 +42,17 @@
 
 EOF
   $FIND $d -name '*.m' | \
-    $PERL -ne 'm{(.*)/(.*)\.m};
-               for (qx{./gethelp $2 < $_}) {
+    $PERL -n -e 'chop;
+               $f = "$_";
+               m{(.*)/(.*)\.m};
+               for (qx{./gethelp $2 "$f" < "$f"}) {
                  s/^\s+\@/\@/ unless $i_am_in_example;
                  s/^\s+\@group/\@group/;
                  s/^\s+\@end\s+group/\@end\s+group/;
                  $i_am_in_example = 1 if /\s*\@example/;
                  $i_am_in_example = 0 if /\s*\@end\s+example/;
                  print;
-               }'  
+               }'
 else
   echo "gethelp program seems to be missing!" 1>&2
   exit 1
--- a/scripts/optimization/__dogleg__.m
+++ b/scripts/optimization/__dogleg__.m
@@ -53,6 +53,7 @@
         alpha = 0;
       endif
     else
+      alpha = delta / xn;
       snm = 0;
     endif
     ## Form the appropriate convex combination.
--- a/scripts/optimization/__fdjac__.m
+++ b/scripts/optimization/__fdjac__.m
@@ -17,21 +17,19 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{fval}, @var{fjac}]} =  __fdjac__ (@var{fcn}, @var{x}, @var{err})
+## @deftypefn{Function File} {} __fdjac__ (@var{fcn}, @var{x}, @var{fvec}, @var{err})
 ## Undocumented internal function.
 ## @end deftypefn
 
-function [fvec, fjac] = __fdjac__ (fcn, x, err = 0)
+function fjac = __fdjac__ (fcn, x, fvec, err = 0)
   err = sqrt (max (eps, err));
-  fvec = fcn (x);
-  fv = fvec(:);
   h = abs (x(:))*err;
   h(h == 0) = err;
-  fjac = zeros (length (fv), numel (x));
+  fjac = zeros (length (fvec), numel (x));
   for i = 1:numel (x)
     x1 = x;
     x1(i) += h(i);
-    fjac(:,i) = (fcn (x1)(:) - fv) / h(i);
+    fjac(:,i) = (fcn (x1)(:) - fvec) / h(i);
   endfor
 endfunction
 
--- a/scripts/optimization/fsolve.m
+++ b/scripts/optimization/fsolve.m
@@ -72,10 +72,49 @@
 ## @item -3
 ## The trust region radius became excessively small. 
 ## @end table
-##
+## 
 ## Note: If you only have a single nonlinear equation of one variable, using 
 ## @code{fzero} is usually a much better idea.
 ## @seealso{fzero, optimset}
+##
+## Note about user-supplied jacobians:
+## As an inherent property of the algorithm, jacobian is always requested for a
+## solution vector whose residual vector is already known, and it is the last
+## accepted successful step. Often this will be one of the last two calls, but
+## not always. If the savings by reusing intermediate results from residual
+## calculation in jacobian calculation are significant, the best strategy is to
+## employ OutputFcn: After a vector is evaluated for residuals, if OutputFcn is
+## called with that vector, then the intermediate results should be saved for
+## future jacobian evaluation, and should be kept until a jacobian evaluation
+## is requested or until outputfcn is called with a different vector, in which
+## case they should be dropped in favor of this most recent vector. A short
+## example how this can be achieved follows:
+##
+## @example
+## function [fvec, fjac] = my_optim_func (x, optimvalues, state)
+## persistent sav = [], sav0 = [];
+## if (nargin == 1)
+##   ## evaluation call
+##   if (nargout == 1)
+##     sav0.x = x; # mark saved vector
+##     ## calculate fvec, save results to sav0.
+##   elseif (nargout == 2)
+##     ## calculate fjac using sav.
+##   endif
+## else
+##   ## outputfcn call.
+##   if (all (x == sav0.x))
+##     sav = sav0;
+##   endif
+##   ## maybe output iteration status etc.
+## endif
+## endfunction
+##
+## ## ....
+## 
+## fsolve (@@my_optim_func, x0, optimset ("OutputFcn", @@my_optim_func, @dots{}))
+## @end example
+###
 ## @end deftypefn
 
 ## PKG_ADD: __all_opts__ ("fsolve");
@@ -125,36 +164,53 @@
   tolf = optimget (options, "TolFun", sqrt (macheps));
 
   factor = 100;
-  ## FIXME: TypicalX corresponds to user scaling (???)
-  autodg = true;
 
   niter = 1;
-  nfev = 0;
+  nfev = 1;
 
   x = x0(:);
   info = 0;
 
+  ## Initial evaluation.
+  ## Handle arbitrary shapes of x and f and remember them.
+  fvec = fcn (reshape (x, xsiz));
+  fsiz = size (fvec);
+  fvec = fvec(:);
+  fn = norm (fvec);
+  m = length (fvec);
+  n = length (x);
+
+  if (! isempty (outfcn))
+    optimvalues.iter = niter;
+    optimvalues.funccount = nfev;
+    optimvalues.fval = fn;
+    optimvalues.searchdirection = zeros (n, 1);
+    state = 'init';
+    stop = outfcn (x, optimvalues, state);
+    if (stop)
+      info = -1;
+      break;
+    endif
+  endif
+
+  nsuciter = 0;
+
   ## Outer loop.
   while (niter < maxiter && nfev < maxfev && ! info)
 
     ## Calculate function value and Jacobian (possibly via FD).
-    ## Handle arbitrary shapes of x and f and remember them.
     if (has_jac)
       [fvec, fjac] = fcn (reshape (x, xsiz));
       ## If the jacobian is sparse, disable Broyden updating.
       if (issparse (fjac))
         updating = false;
       endif
+      fvec = fvec(:);
       nfev ++;
     else
-      [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
-      nfev += 1 + length (x);
+      fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec);
+      nfev += length (x);
     endif
-    fsiz = size (fvec);
-    fvec = fvec(:);
-    fn = norm (fvec);
-    m = length (fvec);
-    n = length (x);
 
     ## For square and overdetermined systems, we update a QR
     ## factorization of the jacobian to avoid solving a full system in each
@@ -173,45 +229,34 @@
     ## Get column norms, use them as scaling factors.
     jcn = norm (fjac, 'columns').';
     if (niter == 1)
-      if (autodg)
-        dg = jcn;
-        dg(dg == 0) = 1;
-      endif
+      dg = jcn;
+      dg(dg == 0) = 1;
       xn = norm (dg .* x);
       ## FIXME: something better?
       delta = max (factor * xn, 1);
     endif
 
-    ## Rescale if necessary.
-    if (autodg)
-      ## FIXME: the original minpack used the following rescaling strategy:
-      ##   dg = max (dg, jcn);
-      ## but it seems not good if we start with a bad guess yielding jacobian
-      ## columns with large norms that later decrease, because the corresponding
-      ## variable will still be overscaled. So instead, we only give the old
-      ## scaling a small momentum, but do not honor it.
+    ## Rescale adaptively.
+    ## FIXME: the original minpack used the following rescaling strategy:
+    ##   dg = max (dg, jcn);
+    ## but it seems not good if we start with a bad guess yielding jacobian
+    ## columns with large norms that later decrease, because the corresponding
+    ## variable will still be overscaled. So instead, we only give the old
+    ## scaling a small momentum, but do not honor it.
 
-      dg = max (0.1*dg, jcn);
+    dg = max (0.1*dg, jcn);
 
-      ## It also seems that in the case of fast (and inhomogeneously) changing
-      ## jacobian, the Broyden updates are of little use, so maybe we could
-      ## skip them if a big disproportional change is expected. The question is,
-      ## of course, how to define the above terms :)
-    endif
+    ## It also seems that in the case of fast (and inhomogeneously) changing
+    ## jacobian, the Broyden updates are of little use, so maybe we could
+    ## skip them if a big disproportional change is expected. The question is,
+    ## of course, how to define the above terms :)
 
+    lastratio = 0;
     nfail = 0;
     nsuc = 0;
     ## Inner loop.
     while (niter <= maxiter && nfev < maxfev && ! info)
 
-      if (useqr)
-        tr_mat = r;
-        tr_vec = q'*fvec;
-      else
-        tr_mat = fjac;
-        tr_vec = fvec;
-      endif
-
       ## Get trust-region model (dogleg) minimizer.
       if (useqr)
         qtf = q'*fvec;
@@ -249,7 +294,7 @@
       endif
 
       ## Update delta.
-      if (ratio < 0.1)
+      if (ratio < min(max(0.1, lastratio), 0.9))
         nsuc = 0;
         nfail ++;
         delta *= 0.5;
@@ -259,6 +304,7 @@
           break;
         endif
       else
+        lastratio = ratio;
         nfail = 0;
         nsuc ++;
         if (abs (1-ratio) <= 0.1)
@@ -274,6 +320,7 @@
         xn = norm (dg .* x);
         fvec = fvec1;
         fn = fn1;
+        nsuciter ++;
       endif
 
       niter ++;
@@ -321,7 +368,7 @@
       endif
 
       ## Criterion for recalculating jacobian.
-      if (! updating || nfail == 2)
+      if (! updating || nfail == 2 || nsuciter < 2)
         break;
       endif
 
@@ -347,6 +394,7 @@
   fvec = reshape (fvec, fsiz);
 
   output.iterations = niter;
+  output.successful = nsuciter;
   output.funcCount = nfev;
 
 endfunction
@@ -358,12 +406,12 @@
     [fx, jx] = fun (x);
   else
     fx = fun (x);
-    jx = []
+    jx = [];
   endif
 
-  if (! complexeqn && ! (all (isreal (fx(:))) && all (isreal (jx(:)))))
+  if (! complexeqn && ! (isreal (fx) && isreal (jx)))
     error ("fsolve:notreal", "fsolve: non-real value encountered"); 
-  elseif (complexeqn && ! (all (isnumeric (fx(:))) && all (isnumeric(jx(:)))))
+  elseif (complexeqn && ! (isnumeric (fx) && isnumeric(jx)))
     error ("fsolve:notnum", "fsolve: non-numeric value encountered");
   elseif (any (isnan (fx(:))))
     error ("fsolve:isnan", "fsolve: NaN value encountered"); 
--- a/scripts/plot/Makefile.in
+++ b/scripts/plot/Makefile.in
@@ -47,6 +47,7 @@
   __errcomm__.m \
   __errplot__.m \
   __ezplot__.m \
+  __gnuplot_get_var__.m \
   __gnuplot_has_feature__.m \
   __go_close_all__.m \
   __go_draw_axes__.m \
new file mode 100644
--- /dev/null
+++ b/scripts/plot/__gnuplot_get_var__.m
@@ -0,0 +1,151 @@
+## Copyright (C) 2009 Ben Abbott
+##
+## 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 3 of the License, 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, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{value} =} __gnuplot_get_var__ (@var{h}, @var{name}, @var{fmt})
+## Undocumented internal function.
+## @end deftypefn
+
+## Author: Ben Abbott <bpabbott@mac.com>
+## Created: 2009-02-07
+
+function gp_var_value = __gnuplot_get_var__ (h, gp_var_name, fmt)
+
+  if (nargin == 0)
+    h = gcf ();
+  endif
+  if (nargin < 2)
+    print_usage ();
+  endif
+  if (nargin < 3)
+    fmt = '';
+  endif
+
+  if (numel (h) == 1 && isfigure (h))
+    ostream = get (h, "__plot_stream__");
+  else
+    ostream = h;
+  endif
+  if (numel (ostream) < 1)
+    error ("__gnuplot_get_var__: stream to gnuplot not open");
+  elseif (ispc ()) # || true
+    if (numel (ostream) == 1)
+      error ("__gnuplot_get_var__: Need mkfifo that is not implemented under Windows");
+    endif
+    use_mkfifo = false;
+    istream = ostream(2);
+    ostream = ostream(1);
+  else
+    use_mkfifo = true;
+    ostream = ostream(1);
+  endif
+
+  if (use_mkfifo)
+    gpin_name = tmpnam ();
+
+    ## Mode: 6*8*8 ==  0600
+    [err, msg] = mkfifo (gpin_name, 6*8*8);
+
+    if (err != 0)
+      error ("__gnuplot_get_var__: Can not open fifo (%s)", msg);
+    endif
+  endif
+
+  gp_var_name = strtrim (gp_var_name);
+  n = min (strfind (gp_var_name, " "), strfind (gp_var_name, ",")) - 1;
+  if (isempty (n))
+    n = numel (gp_var_name);
+  endif
+
+  unwind_protect
+
+    ## Notes: Variables may be undefined if user closes gnuplot by "q"
+    ## or Alt-F4. Further, this abrupt close also requires the leading
+    ## "\n" on the next line.
+    if (use_mkfifo)
+      fprintf (ostream, "\nset print \"%s\";\n", gpin_name);
+      fflush (ostream);
+      [gpin, err] = fopen (gpin_name, "r");
+      if (err != 0)
+        error ("__gnuplot_get_var__: Can not open fifo (%s)", msg);
+      endif
+      gp_cmd = sprintf ("\nif (exists(\"%s\")) print %s; else print NaN\n",
+                        gp_var_name(1:n), gp_var_name);
+      fputs (ostream, gp_cmd);
+
+      ## Close output file, to force it to be flushed
+      fputs (ostream, "set print;\n");
+      fflush (ostream);
+
+      ## Now read from fifo.
+      reading = true;
+      str = {};
+      while (reading)
+        str{end+1} = fgets (gpin);
+        if (isnumeric (str{end}) && (str{end} == -1))
+          reading = false;
+          str = str(1:(end-1));
+        endif
+      endwhile
+      str = strcat (str{:});
+      fclose (gpin);
+    else
+      ## Direct gnuplot to print to <STDOUT>
+      fprintf (ostream, "set print \"-\";\n");
+      fflush (ostream);
+      gp_cmd = sprintf ("\nif (exists(\"%s\")) print \"OCTAVE: \", %s; else print NaN\n",
+                        gp_var_name(1:n), gp_var_name);
+      fputs (ostream, gp_cmd);
+      fflush (ostream);
+      ## Direct gnuplot to print to <STDERR>
+      fputs (ostream, "set print;\n");
+      fflush (ostream);
+
+      str = {};
+      while (isempty (str))
+        str = char (fread (istream)');
+        if (! isempty (str))
+          str = regexp (str, "OCTAVE:.*", "match")
+          str = str{end}(8:end);
+        endif
+        fclear (istream);
+      endwhile
+    endif
+
+    ## Strip out EOLs and the continuation character "|"
+    str(str=="\n") = "";
+    str(str=="\r") = "";
+    n_continue = strfind (str, " \\ ");
+    if (! isempty (n_continue))
+      str(n_continue+1) = "";
+    endif
+
+    if (isempty (fmt))
+      gp_var_value = strtrim (str);
+    else
+      gp_var_value = sscanf (str, fmt);
+    endif
+
+  unwind_protect_cleanup
+    if (use_mkfifo)
+      unlink (gpin_name);
+    endif
+  end_unwind_protect
+
+endfunction
+
--- a/scripts/plot/__gnuplot_has_feature__.m
+++ b/scripts/plot/__gnuplot_has_feature__.m
@@ -29,12 +29,13 @@
               "transparent_surface",
               "epslatex_implies_eps_filesuffix",
               "epslatexstandalone_terminal",
-              "screen_coordinates_for_{lrtb}margin",};
+              "screen_coordinates_for_{lrtb}margin",
+              "variable_GPVAL_TERMINALS"};
 
   if (isempty (has_features))
     gnuplot_version = __gnuplot_version__ ();
-    versions = {"4.2.4", "4.3", "4.3", "4.2", "4.2", "4.3"};
-    operators = {">", ">=", ">=", ">=", ">=", ">="};
+    versions = {"4.2.4", "4.3", "4.3", "4.2", "4.2", "4.3", "4.3"};
+    operators = {">", ">=", ">=", ">=", ">=", ">=", ">="};
     have_features = logical (zeros (size (features)));
     for n = 1 : numel (have_features)
       has_features(n) = compare_versions (gnuplot_version, versions{n}, operators{n});
--- a/scripts/plot/gnuplot_drawnow.m
+++ b/scripts/plot/gnuplot_drawnow.m
@@ -40,17 +40,34 @@
     fid = [];
     printing = ! output_to_screen (gnuplot_trim_term (term));
     unwind_protect
-      plot_stream = open_gnuplot_stream (1, []);
-      [enhanced, implicit_margin] = gnuplot_set_term (plot_stream (1), true, h, term, file);
-      __go_draw_figure__ (h, plot_stream, enhanced, mono, printing, implicit_margin);
-      if (nargin == 5)
-        fid = fopen (debug_file, "wb");
-        enhanced = gnuplot_set_term (fid, true, h, term, file);
-        __go_draw_figure__ (h, fid, enhanced, mono, printing, implicit_margin);
+      plot_stream = open_gnuplot_stream (2, []);
+      if (__gnuplot_has_feature__ ("variable_GPVAL_TERMINALS"))
+        available_terminals = __gnuplot_get_var__ (plot_stream, "GPVAL_TERMINALS");
+        available_terminals = regexp (available_terminals, "\\b\\w+\\b", "match");
+        gnuplot_supports_term = any (strcmpi (available_terminals,
+                                              gnuplot_trim_term (term)));
+      else
+        gnuplot_supports_term = true;
+      endif
+      if (gnuplot_supports_term)
+        [enhanced, implicit_margin] = gnuplot_set_term (plot_stream (1), true,
+                                                        h, term, file);
+        __go_draw_figure__ (h, plot_stream(1), enhanced, mono, printing, implicit_margin);
+        if (nargin == 5)
+          fid = fopen (debug_file, "wb");
+          [enhanced, implicit_margin] = gnuplot_set_term (fid, true, h, term, file);
+          __go_draw_figure__ (h, fid, enhanced, mono, printing, implicit_margin);
+        endif
+      else
+        error ("gnuplot_drawnow: the gnuplot terminal, \"%s\", is not available.",
+               gnuplot_trim_term (term))
       endif
     unwind_protect_cleanup
       if (! isempty (plot_stream))
-        pclose (plot_stream);
+        pclose (plot_stream(1));
+        if (numel (plot_stream) == 2)
+          pclose (plot_stream(2));
+        endif
       endif
       if (! isempty (fid))
         fclose (fid);
--- a/scripts/plot/print.m
+++ b/scripts/plot/print.m
@@ -427,11 +427,36 @@
     endif
     set (gcf, "__pixels_per_inch__", resolution)
 
-    if (debug)
-      drawnow (new_terminal, name, mono, debug_file);
-    else
-      drawnow (new_terminal, name, mono);
-    endif
+    unwind_protect
+      if (! isempty (size))
+        size_in_pixels = sscanf (size ,"%d, %d");
+        size_in_pixels = reshape (size_in_pixels, [1, numel(size_in_pixels)]);
+        size_in_inches = size_in_pixels ./ resolution;
+        p.paperunits = get (gcf, "paperunits");
+        p.papertype = get (gcf, "papertype");
+        p.papersize = get (gcf, "papersize");
+        p.paperposition = get (gcf, "paperposition");
+        p.paperpositionmode = get (gcf, "paperpositionmode");
+        set (gcf, "paperunits", "inches");
+        set (gcf, "papertype", "<custom>");
+        set (gcf, "papersize", size_in_inches);
+        set (gcf, "paperposition", [0, 0, size_in_inches]);
+        set (gcf, "paperpositionmode", "manual");
+      endif
+      if (debug)
+        drawnow (new_terminal, name, mono, debug_file);
+      else
+        drawnow (new_terminal, name, mono);
+      endif
+    unwind_protect_cleanup
+      ## FIXME - it would be preferred to delete the added properties here.
+      if (! isempty (size))
+        props = fieldnames (p);
+        for n = 1:numel(props)
+          set (gcf, props{n}, p.(props{n}))
+        endfor
+      endif
+    end_unwind_protect
 
     if (! isempty (convertname))
       command = sprintf ("convert '%s' '%s'", name, convertname);
--- a/scripts/testfun/assert.m
+++ b/scripts/testfun/assert.m
@@ -141,7 +141,7 @@
       iserror = 1;
       coda = "Dimensions don't match";
 
-    elseif (tol == 0 && ! strcmp (typeinfo (cond), typeinfo (expected)))
+    elseif (nargin < 3 && ! strcmp (typeinfo (cond), typeinfo (expected)))
       iserror = 1;
       coda = cstrcat ("Type ", typeinfo (cond), " != ", typeinfo (expected));
 
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,133 @@
+2009-03-25  John W. Eaton  <jwe@octave.org>
+
+	* DLD-FUNCTIONS/find.cc (find_nonzero_elem_idx): Also return
+	[](0x0) if the array has 0 rows and it is not a column vector.
+
+	* oct-stream.cc (octave_stream::write (const Array<T>&,
+	octave_idx_type, oct_data_conv::data_type, octave_idx_type,
+	oct_mach_info::float_format)): Seek to skip if still inside bounds
+	of existing file.  Otherwise, write NUL to skip.
+
+	* Makefile.in (%.df : %.cc): Write source file name to output,
+	wrapped in XDEFUN_FILE_NAME macro.
+	* mkbuiltins: Provide definition for XDEFUN_FILE_NAME.
+	* mkgendoc: Likeiwse.
+	(XDEFUN_DLD_INTERNAL, XDEFUNX_DLD_INTERNAL, XDEFUN_INTERNAL,
+	XDEFCONSTFUN_INTERNAL, XDEFUNX_INTERNAL, XDEFVAR_INTERNAL,
+	XDEFCONST_INTERNAL): Pass file_name to print_doc_string.
+	(print_doc_string): New arg, FILE_NAME.  Print file name as
+	comment.
+
+2009-03-24  John W. Eaton  <jwe@octave.org>
+
+	* ov-class.cc (F__parent_classes__): New function.
+
+2009-03-24  Robert T. Short  <octave@phaselockedsystems.com>
+
+	* ov-class.h, ov-class.cc (octave_class::octave_class (const
+	Octave_map&, const std::string&, const octave_value_list&)):
+	New constructor.
+	(octave_class::find_parent_class, octave_class::parent_classes):
+	New functions.
+	(octave_class::dotref): Also look in parent class.
+	(Fclass): Handle parent class arguments.
+
+	* ov-base.h (octave_base_value::get_parent_list,
+	octave_base_value::parent_classes): New virtual functions.
+
+	* load-path.h, load-path.cc (load_path::parent_map): New data member. 
+	(load_path::add_to_parent_map): New static function.
+	(load_path::do_add_to_parent_map): New member function.
+	(load_path::do_find_method): Also look in parent classes for methods.
+	(load_path::parent_map_type, load_path::const_parent_map_iterator,
+	load_path::parent_map_iterator): New typedefs.
+
+2009-03-23  John W. Eaton  <jwe@octave.org>
+
+	* symtab.h
+	(symbol_table::symbol_record::symobl_recoord_rep::is_variable):
+	Also return true if symbol is tagged as a variable.
+	(symbol_table::symbol_record::symobl_recoord_rep::force_variable):
+	Don't set variable value.
+	(symbol_table::symbol_record::symobl_recoord_rep::clear_forced,
+	symbol_table::symbol_record::clear_forced): Delete.
+	(symbol_table::unmark_forced_variables): Rename from
+	symbol_table::clear_forced_variables.
+	(symbol_table::do_unmark_forced_variables): Rename from
+	symbol_table::do_clear_forced_variables.
+	* parse.y (make_script, finish_function): Call
+	symbol_table::unmark_forced_variables instead of
+	symbol_table::clear_forced_variables.
+	* octave.cc (unmark_forced_vars): New function.
+	(execute_eval_option_code): Add it to the unwind-protect stack.
+
+2009-03-22  Jaroslav Hajek  <highegg@gmail.com>
+
+	* pt-eval.cc (tree_evaluator::visit_simple_for_command):
+	Remove struct branch, handle structs by the generic code.
+	(tree_evaluator::visit_complex_for_command):
+	Add missing const qualifiers.
+
+2009-03-21  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-base-diag.cc: Add missing include.
+	* sparse-xdiv.cc: Ditto.
+	* DLD-FUNCTIONS/__glpk__.cc: Ditto.
+	* DLD-FUNCTIONS/__lin_interpn__.cc: Ditto.
+	* DLD-FUNCTIONS/__voronoi__.cc: Ditto.
+
+2009-03-20  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-re-mat.cc (octave_matrix::load_ascii): Simplify.
+	* ov-flt-re-mat.cc (octave_float_matrix::load_ascii): Simplify.
+	* ov-cx-mat.cc (octave_complex_matrix::load_ascii): Simplify.
+	* ov-flt-cx-mat.cc (octave_float_complex_matrix::load_ascii): Simplify.
+
+2009-03-20  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-re-mat.cc (octave_matrix::isnan, octave_matrix::isinf,
+	octave_matrix::finite): Simplify.
+	* ov-flt-re-mat.cc (octave_float_matrix::isnan,
+	octave_float_matrix::isinf, octave_float_matrix::finite): Simplify.
+	* ov-cx-mat.cc (octave_complex_matrix::isnan,
+	octave_complex_matrix::isinf, octave_complex_matrix::finite):
+	Simplify.
+	* ov-flt-cx-mat.cc (octave_float_complex_matrix::isnan,
+	octave_float_complex_matrix::isinf,
+	octave_float_complex_matrix::finite): Simplify.
+
+2009-03-19  Benjamin Lindner  <lindnerb@users.sourceforge.net>
+
+	* ls-oct-ascii.cc (extract_keyword): Replace loop with call to
+	read_until_newline to avoid leaving stray '\r' in stream when
+	reading files with CRLF line endings.
+
+2009-03-18  Jaroslav Hajek <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/cellfun.cc
+	(scalar_col_helper, scalar_col_helper_def, scalar_col_helper_nda,
+	scalar_query_helper): New classes.
+	(make_col_helper): New function.
+	(Fcellfun): Use col_helper classes for UniformOutput. Avoid copying by
+	using const variables.
+
+2009-03-17  Jaroslav Hajek  <highegg@gmail.com>
+
+	* OPERATORS/op-b-bm.cc: Add missing INSTALL_BINOPs.
+	* OPERATORS/op-bm-b.cc: Ditto.
+	* OPERATORS/op-fm-fm.cc: Ditto.
+	* OPERATORS/op-m-m.cc: Ditto.
+
+2009-03-17  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov.cc (octave_value::octave_value): Move to ov.h
+	* ov.h (octave_value::octave_value): Implement fast inline constructor
+	using a static rep.
+
+2009-03-15  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/cellfun.cc (Fcellslices): New DLD function.
+
 2009-03-13  Jaroslav Hajek  <highegg@gmail.com>
 
 	* ov.h (octave_value::compound_binary_op): Support bool compound ops.
--- a/src/DLD-FUNCTIONS/__glpk__.cc
+++ b/src/DLD-FUNCTIONS/__glpk__.cc
@@ -28,6 +28,8 @@
 #include <csetjmp>
 #include <ctime>
 
+#include "lo-ieee.h"
+
 #include "defun-dld.h"
 #include "error.h"
 #include "gripes.h"
--- a/src/DLD-FUNCTIONS/__lin_interpn__.cc
+++ b/src/DLD-FUNCTIONS/__lin_interpn__.cc
@@ -24,6 +24,7 @@
 #include <config.h>
 #endif
 
+#include "lo-ieee.h"
 #include "dNDArray.h"
 #include "oct-locbuf.h"
 
--- a/src/DLD-FUNCTIONS/__voronoi__.cc
+++ b/src/DLD-FUNCTIONS/__voronoi__.cc
@@ -34,6 +34,7 @@
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
+#include "lo-ieee.h"
 #include "oct.h"
 #include "Cell.h"
 
--- a/src/DLD-FUNCTIONS/cellfun.cc
+++ b/src/DLD-FUNCTIONS/cellfun.cc
@@ -1,6 +1,7 @@
 /*
 
 Copyright (C) 2005, 2006, 2007, 2008, 2009 Mohamed Kamoun
+Copyright (C) 2009 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -27,6 +28,7 @@
 #include <string>
 #include <vector>
 #include <list>
+#include <memory>
 
 #include "lo-mappers.h"
 #include "oct-locbuf.h"
@@ -39,6 +41,142 @@
 #include "ov-colon.h"
 #include "unwind-prot.h"
 
+// Rationale:
+// The octave_base_value::subsasgn method carries too much overhead for
+// per-element assignment strategy.
+// This class will optimize the most optimistic and most likely case
+// when the output really is scalar by defining a hierarchy of virtual
+// collectors specialized for some scalar types.
+
+class scalar_col_helper
+{
+public:
+  virtual bool collect (octave_idx_type i, const octave_value& val) = 0;
+  virtual octave_value result (void) = 0;
+  virtual ~scalar_col_helper (void) { }
+};
+
+// The default collector represents what was previously done in the main loop.
+// This reuses the existing assignment machinery via octave_value::subsasgn,
+// which can perform all sorts of conversions, but is relatively slow.
+
+class scalar_col_helper_def : public scalar_col_helper
+{
+  std::list<octave_value_list> idx_list;
+  octave_value resval;
+public:
+  scalar_col_helper_def (const octave_value& val, const dim_vector& dims)
+    : idx_list (1), resval (val)
+    {
+      idx_list.front ().resize (1);
+      if (resval.dims () != dims)
+        resval.resize (dims);
+    }
+  ~scalar_col_helper_def (void) { }
+
+  bool collect (octave_idx_type i, const octave_value& val)
+    {
+      if (val.numel () == 1)
+        {
+          idx_list.front ()(0) = static_cast<double> (i + 1);
+          resval = resval.subsasgn ("(", idx_list, val);
+        }
+      else
+        error ("cellfun: expecting all values to be scalars for UniformOutput = true");
+
+      return true;
+    }
+  octave_value result (void)
+    {
+      return resval;
+    }
+};
+
+template <class T>
+struct scalar_query_helper { };
+
+#define DEF_QUERY_HELPER(T, TEST, QUERY) \
+template <> \
+struct scalar_query_helper<T> \
+{ \
+  static bool has_value (const octave_value& val) \
+    { return TEST; } \
+  static T get_value (const octave_value& val) \
+    { return QUERY; } \
+}
+
+DEF_QUERY_HELPER (double, val.is_real_scalar (), val.scalar_value ());
+DEF_QUERY_HELPER (Complex, val.is_complex_scalar (), val.complex_value ());
+DEF_QUERY_HELPER (float, val.is_single_type () && val.is_real_scalar (), 
+                  val.float_scalar_value ());
+DEF_QUERY_HELPER (FloatComplex, val.is_single_type () && val.is_complex_scalar (), 
+                  val.float_complex_value ());
+DEF_QUERY_HELPER (bool, val.is_bool_scalar (), val.bool_value ());
+// FIXME: More?
+
+// This specializes for collecting elements of a single type, by accessing
+// an array directly. If the scalar is not valid, it returns false.
+
+template <class NDA>
+class scalar_col_helper_nda : public scalar_col_helper
+{
+  NDA arrayval;
+  typedef typename NDA::element_type T;
+public:
+  scalar_col_helper_nda (const octave_value& val, const dim_vector& dims)
+    : arrayval (dims)
+    {
+      arrayval(0) = scalar_query_helper<T>::get_value (val);
+    }
+  ~scalar_col_helper_nda (void) { }
+
+  bool collect (octave_idx_type i, const octave_value& val)
+    {
+      bool retval = scalar_query_helper<T>::has_value (val);
+      if (retval)
+        arrayval(i) = scalar_query_helper<T>::get_value (val);
+      return retval;
+    }
+  octave_value result (void)
+    {
+      return arrayval;
+    }
+};
+
+template class scalar_col_helper_nda<NDArray>;
+template class scalar_col_helper_nda<FloatNDArray>;
+template class scalar_col_helper_nda<ComplexNDArray>;
+template class scalar_col_helper_nda<FloatComplexNDArray>;
+template class scalar_col_helper_nda<boolNDArray>;
+
+// the virtual constructor.
+scalar_col_helper *
+make_col_helper (const octave_value& val, const dim_vector& dims)
+{
+  scalar_col_helper *retval;
+
+  if (val.is_bool_scalar ())
+    retval = new scalar_col_helper_nda<boolNDArray> (val, dims);
+  else if (val.is_complex_scalar ())
+    {
+      if (val.is_single_type ())
+        retval = new scalar_col_helper_nda<FloatComplexNDArray> (val, dims);
+      else
+        retval = new scalar_col_helper_nda<ComplexNDArray> (val, dims);
+    }
+  else if (val.is_real_scalar ())
+    {
+      if (val.is_single_type ())
+        retval = new scalar_col_helper_nda<FloatNDArray> (val, dims);
+      else
+        retval = new scalar_col_helper_nda<NDArray> (val, dims);
+    }
+  else
+    retval = new scalar_col_helper_def (val, dims);
+
+  return retval;
+}
+
 DEFUN_DLD (cellfun, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {} cellfun (@var{name}, @var{c})\n\
@@ -163,7 +301,7 @@
       return retval;
     }
   
-  Cell f_args = args(1).cell_value ();
+  const Cell f_args = args(1).cell_value ();
   
   octave_idx_type k = f_args.numel ();
 
@@ -202,7 +340,7 @@
         result(count) = static_cast<double> (f_args.elem(count).ndims ());
       retval(0) = result;
     }
-  else if (name == "prodofsize")
+  else if (name == "prodofsize" || name == "numel")
     {
       NDArray result (f_args.dims ());
       for (octave_idx_type count = 0; count < k ; count++)
@@ -270,7 +408,6 @@
 	error ("unknown function");
       else
 	{
-	  octave_value_list idx;
 	  octave_value_list inputlist;
 	  bool uniform_output = true;
 	  bool have_error_handler = false;
@@ -279,6 +416,8 @@
 	  int offset = 1;
 	  int i = 1;
 	  OCTAVE_LOCAL_BUFFER (Cell, inputs, nargin);
+          // This is to prevent copy-on-write.
+          const Cell *cinputs = inputs;
 
 	  while (i < nargin)
 	    {
@@ -344,19 +483,20 @@
 		}
 	    }
 
-	  inputlist.resize(nargin-offset);
+          nargin -= offset;
+	  inputlist.resize(nargin);
 
 	  if (have_error_handler)
 	    buffer_error_messages++;
 
 	  if (uniform_output)
 	    {
-	      retval.resize(nargout);
+              OCTAVE_LOCAL_BUFFER (std::auto_ptr<scalar_col_helper>, retptr, nargout);
 
 	      for (octave_idx_type count = 0; count < k ; count++)
 		{
-		  for (int j = 0; j < nargin-offset; j++)
-		    inputlist(j) = inputs[j](count);
+		  for (int j = 0; j < nargin; j++)
+		    inputlist(j) = cinputs[j](count);
 
 		  octave_value_list tmp = feval (func, inputlist, nargout);
 
@@ -390,38 +530,46 @@
 		    {
 		      for (int j = 0; j < nargout; j++)
 			{
-			  octave_value val;
-			  val = tmp(j);
+			  octave_value val = tmp(j);
 
-			  if (error_state)
-			    goto cellfun_err;
-
-			  retval(j) = val.resize(f_args.dims());
+                          if (val.numel () == 1)
+                            retptr[j].reset (make_col_helper (val, f_args.dims ()));
+                          else
+                            {
+                              error ("cellfun: expecting all values to be scalars for UniformOutput = true");
+                              break;
+                            }
 			}
 		    }
 		  else
 		    {
-		      idx(0) = octave_value (static_cast<double>(count+1));
 		      for (int j = 0; j < nargout; j++)
 			{
-			  // FIXME -- need an easier way to express
-			  // this test.
 			  octave_value val = tmp(j);
 
-			  if (val.ndims () == 2
-			      && val.rows () == 1 && val.columns () == 1)
-			    retval(j) = 
-			      retval(j).subsasgn ("(", 
-						  std::list<octave_value_list> 
-						  (1, idx(0)), val);
-			  else
-			    error ("cellfun: expecting all values to be scalars for UniformOutput = true");
-			}
+                          if (! retptr[j]->collect (count, val))
+                            {
+                              // FIXME: A more elaborate structure would allow again a virtual
+                              // constructor here.
+                              retptr[j].reset (new scalar_col_helper_def (retptr[j]->result (), 
+                                                                          f_args.dims ()));
+                              retptr[j]->collect (count, val);
+                            }
+                        }
 		    }
 
 		  if (error_state)
 		    break;
 		}
+
+              retval.resize (nargout);
+              for (int j = 0; j < nargout; j++)
+                {
+                  if (retptr[j].get ())
+                    retval(j) = retptr[j]->result ();
+                  else
+                    retval(j) = Matrix ();
+                }
 	    }
 	  else
 	    {
@@ -431,8 +579,8 @@
 
 	      for (octave_idx_type count = 0; count < k ; count++)
 		{
-		  for (int j = 0; j < nargin-offset; j++)
-		    inputlist(j) = inputs[j](count);
+		  for (int j = 0; j < nargin; j++)
+		    inputlist(j) = cinputs[j](count);
 
 		  octave_value_list tmp = feval (func, inputlist, nargout);
 
@@ -1069,6 +1217,97 @@
 %! assert(c,{empty1by0str,'abcd','ef',empty1by0str,'ghij',empty1by0str})
 
 */
+
+template <class NDA>
+Cell 
+do_cellslices_nda (const NDA& array, const idx_vector& lb, const idx_vector& ub)
+{
+  octave_idx_type n = lb.length (0);
+  Cell retval (1, n);
+  for (octave_idx_type i = 0; i < n && ! error_state; i++)
+    retval(i) = array.index (idx_vector (lb(i), ub(i) + 1));
+  return retval;
+}
+
+DEFUN_DLD (cellslices, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{sl} =} cellslices (@var{x}, @var{lb}, @var{ub})\n\
+Given a vector @var{x}, this function produces a cell array of slices from the vector\n\
+determined by the index vectors @var{lb}, @var{ub}, for lower and upper bounds, respectively.\n\
+In other words, it is equivalent to the following code:\n\
+\n\
+@example\n\
+n = length (lb);\n\
+sl = cell (1, n);\n\
+for i = 1:length (lb)\n\
+  sl@{i@} = x(lb(i):ub(i));\n\
+endfor\n\
+@end example\n\
+\n\
+If @var{X} is a matrix, linear indexing will be used.\n\
+@seealso{mat2cell}\n\
+@end deftypefn")
+{
+  octave_value retval;
+  int nargin = args.length ();
+  if (nargin == 3)
+    {
+      octave_value x = args(0);
+      idx_vector lb = args(1).index_vector (), ub = args(2).index_vector ();
+      if (! error_state)
+        {
+          if (lb.is_colon () || ub.is_colon ())
+            error ("cellslices: invalid use of colon");
+          else if (lb.length (0) != ub.length (0))
+            error ("cellslices: the lengths of lb and ub must match");
+          else
+            {
+              Cell retcell;
+              if (! x.is_sparse_type () && x.is_matrix_type ())
+                {
+                  // specialize for some dense arrays.
+                  if (x.is_bool_type ())
+                    retcell = do_cellslices_nda (x.bool_array_value (), lb, ub);
+                  else if (x.is_char_matrix ())
+                    retcell = do_cellslices_nda (x.char_array_value (), lb, ub);
+                  else if (x.is_complex_type ())
+                    {
+                      if (x.is_single_type ())
+                        retcell = do_cellslices_nda (x.float_complex_array_value (), lb, ub);
+                      else
+                        retcell = do_cellslices_nda (x.complex_array_value (), lb, ub);
+                    }
+                  else
+                    {
+                      if (x.is_single_type ())
+                        retcell = do_cellslices_nda (x.float_array_value (), lb, ub);
+                      else
+                        retcell = do_cellslices_nda (x.array_value (), lb, ub);
+                    }
+                }
+              else
+                {
+                  // generic code.
+                  octave_idx_type n = lb.length (0);
+                  retcell = Cell (1, n);
+                  octave_value_list idx (1, octave_value ());
+                  for (octave_idx_type i = 0; i < n && ! error_state; i++)
+                    {
+                      idx(0) = Range (static_cast<double> (lb(i)) + 1,
+                                      static_cast<double> (ub(i)) + 1);
+                      retcell(i) = x.do_index_op (idx);
+                    }
+                }
+              if (! error_state)
+                retval = retcell;
+            }
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
 	  
 /*
 ;;; Local Variables: ***
--- a/src/DLD-FUNCTIONS/find.cc
+++ b/src/DLD-FUNCTIONS/find.cc
@@ -89,14 +89,23 @@
   octave_idx_type result_nr = count;
   octave_idx_type result_nc = 1;
 
+  bool column_vector_arg = false;
   bool scalar_arg = false;
 
-  if (nda.ndims () == 2 && nda.rows () == 1)
+  if (nda.ndims () == 2)
     {
-      result_nr = 1;
-      result_nc = count;
+      octave_idx_type nr = nda.rows ();
+      octave_idx_type nc = nda.columns ();
 
-      scalar_arg = (nda.columns () == 1);
+      if (nr == 1)
+	{
+	  result_nr = 1;
+	  result_nc = count;
+
+	  scalar_arg = (nc == 1);
+	}
+      else if (nc == 1)
+	column_vector_arg = true;
     }
 
   Matrix idx (result_nr, result_nc);
@@ -141,7 +150,7 @@
 	  i++;
 	}
     }
-  else if (scalar_arg)
+  else if (scalar_arg || (nda.rows () == 0 && ! column_vector_arg))
     {
       idx.resize (0, 0);
 
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -40,6 +40,7 @@
 %.df : %.cc
 	@echo making $@ from $<
 	@(echo "// DO NOT EDIT!  Generated automatically by mkdefs." ; \
+	  echo " XDEFUN_FILE_NAME (\"$<\")" ; \
 	  egrep '^(///*|/\*) *PKG_ADD:' $< ; \
 	  $(CXXCPP) $(CPPFLAGS) $(CXXFLAGS_NO_PT_FLAGS) -DMAKE_BUILTINS $< \
 	    | $(srcdir)/mkdefs) > $@-t
--- a/src/OPERATORS/op-b-bm.cc
+++ b/src/OPERATORS/op-b-bm.cc
@@ -65,6 +65,8 @@
 {
   INSTALL_BINOP (op_el_and, octave_bool, octave_bool_matrix, el_and);
   INSTALL_BINOP (op_el_or, octave_bool, octave_bool_matrix, el_or);
+  INSTALL_BINOP (op_el_and_not, octave_bool, octave_bool_matrix, el_and_not);
+  INSTALL_BINOP (op_el_or_not, octave_bool, octave_bool_matrix, el_or_not);
 
   INSTALL_CATOP (octave_bool, octave_bool_matrix, b_bm);
   INSTALL_CATOP (octave_bool, octave_matrix, b_m);
--- a/src/OPERATORS/op-bm-b.cc
+++ b/src/OPERATORS/op-bm-b.cc
@@ -86,6 +86,8 @@
 {
   INSTALL_BINOP (op_el_and, octave_bool_matrix, octave_bool, el_and);
   INSTALL_BINOP (op_el_or, octave_bool_matrix, octave_bool, el_or);
+  INSTALL_BINOP (op_el_not_and, octave_bool_matrix, octave_bool, el_not_and);
+  INSTALL_BINOP (op_el_not_or, octave_bool_matrix, octave_bool, el_not_or);
 
   INSTALL_CATOP (octave_bool_matrix, octave_bool, bm_b);
   INSTALL_CATOP (octave_bool_matrix, octave_scalar, bm_s);
--- a/src/OPERATORS/op-fm-fm.cc
+++ b/src/OPERATORS/op-fm-fm.cc
@@ -202,6 +202,10 @@
   INSTALL_BINOP (op_el_ldiv, octave_float_matrix, octave_float_matrix, el_ldiv);
   INSTALL_BINOP (op_el_and, octave_float_matrix, octave_float_matrix, el_and);
   INSTALL_BINOP (op_el_or, octave_float_matrix, octave_float_matrix, el_or);
+  INSTALL_BINOP (op_el_and_not, octave_float_matrix, octave_float_matrix, el_and_not);
+  INSTALL_BINOP (op_el_or_not, octave_float_matrix, octave_float_matrix, el_or_not);
+  INSTALL_BINOP (op_el_not_and, octave_float_matrix, octave_float_matrix, el_not_and);
+  INSTALL_BINOP (op_el_not_or, octave_float_matrix, octave_float_matrix, el_not_or);
   INSTALL_BINOP (op_trans_mul, octave_float_matrix, octave_float_matrix, trans_mul);
   INSTALL_BINOP (op_mul_trans, octave_float_matrix, octave_float_matrix, mul_trans);
   INSTALL_BINOP (op_herm_mul, octave_float_matrix, octave_float_matrix, trans_mul);
--- a/src/OPERATORS/op-m-m.cc
+++ b/src/OPERATORS/op-m-m.cc
@@ -175,6 +175,10 @@
   INSTALL_BINOP (op_el_ldiv, octave_matrix, octave_matrix, el_ldiv);
   INSTALL_BINOP (op_el_and, octave_matrix, octave_matrix, el_and);
   INSTALL_BINOP (op_el_or, octave_matrix, octave_matrix, el_or);
+  INSTALL_BINOP (op_el_and_not, octave_matrix, octave_matrix, el_and_not);
+  INSTALL_BINOP (op_el_or_not, octave_matrix, octave_matrix, el_or_not);
+  INSTALL_BINOP (op_el_not_and, octave_matrix, octave_matrix, el_not_and);
+  INSTALL_BINOP (op_el_not_or, octave_matrix, octave_matrix, el_not_or);
   INSTALL_BINOP (op_trans_mul, octave_matrix, octave_matrix, trans_mul);
   INSTALL_BINOP (op_mul_trans, octave_matrix, octave_matrix, mul_trans);
   INSTALL_BINOP (op_herm_mul, octave_matrix, octave_matrix, trans_mul);
--- a/src/graphics.h.in
+++ b/src/graphics.h.in
@@ -35,6 +35,8 @@
 #include <set>
 #include <string>
 
+#include "lo-ieee.h"
+
 #include "gripes.h"
 #include "oct-map.h"
 #include "oct-mutex.h"
--- a/src/load-path.cc
+++ b/src/load-path.cc
@@ -512,6 +512,7 @@
   fcn_map.clear ();
   private_fcn_map.clear ();
   method_map.clear ();
+  parent_map.clear ();
 
   do_append (".", false);
 }
@@ -1016,7 +1017,29 @@
 
       const_fcn_map_iterator p = m.find (meth);
 
-      if (p != m.end ())
+      if (p == m.end ())
+	{
+	  // Look in parent classes.
+
+	  const_parent_map_iterator r = parent_map.find (class_name);
+
+	  if (r != parent_map.end ())
+	    {
+	      const std::list<std::string>& plist = r->second;
+	      std::list<std::string>::const_iterator it = plist.begin ();
+
+	      while (it != plist.end ())
+		{
+		  retval = do_find_method (*it, meth, dir_name, type);
+
+		  if (retval != "")
+		    break;
+
+		  it++;
+		}
+	    }
+	}
+      else
 	{
 	  const file_info_list_type& file_info_list = p->second;
 
@@ -1538,6 +1561,13 @@
 }
 
 void
+load_path::do_add_to_parent_map (const std::string& classname,
+				 const std::list<std::string>& parent_list) const
+{
+  parent_map[classname] = parent_list;
+}
+
+void
 load_path::add_to_fcn_map (const dir_info& di, bool at_end) const
 {
   std::string dir_name = di.dir_name;
--- a/src/load-path.h
+++ b/src/load-path.h
@@ -37,7 +37,8 @@
 {
 protected:
 
-  load_path (void) : dir_info_list (), fcn_map (), method_map () { }
+  load_path (void)
+    : dir_info_list (), fcn_map (), method_map (), parent_map () { }
 
 public:
 
@@ -228,6 +229,13 @@
     return instance_ok () ? instance->do_system_path () : std::string ();
   }
 
+  static void add_to_parent_map (const std::string& classname,
+				 const std::list<std::string>& parent_list)
+  {
+    if (instance_ok ())
+      instance->do_add_to_parent_map (classname, parent_list);
+  }
+
 private:
 
   static const int M_FILE = 1;
@@ -386,6 +394,12 @@
 
   typedef method_map_type::const_iterator const_method_map_iterator;
   typedef method_map_type::iterator method_map_iterator;
+ 
+  // <CLASS_NAME, PARENT_LIST>>
+  typedef std::map<std::string, std::list<std::string> > parent_map_type;
+
+  typedef parent_map_type::const_iterator const_parent_map_iterator;
+  typedef parent_map_type::iterator parent_map_iterator;
 
   mutable dir_info_list_type dir_info_list;
 
@@ -395,6 +409,8 @@
 
   mutable method_map_type method_map;
 
+  mutable parent_map_type parent_map;
+
   static load_path *instance;
 
   static hook_fcn_ptr add_hook;
@@ -493,6 +509,9 @@
 
   std::string do_get_command_line_path (void) const { return command_line_path; }
 
+  void do_add_to_parent_map (const std::string& classname,
+			     const std::list<std::string>& parent_list) const;
+
   void add_to_fcn_map (const dir_info& di, bool at_end) const;
 
   void add_to_private_fcn_map (const dir_info& di) const;
--- a/src/ls-oct-ascii.cc
+++ b/src/ls-oct-ascii.cc
@@ -110,14 +110,8 @@
 	      while (is.get (c) && (c == ' ' || c == '\t' || c == ':'))
 		; // Skip whitespace and the colon.
 
-	      if (c != '\n' && c != '\r')
-		{
-		  value << c;
-		  while (is.get (c) && c != '\n' && c != '\r')
-		    value << c;
-		}
-
-	      retval = value.str ();
+	      is.putback(c);
+	      retval = read_until_newline (is, false);
 	      break;
 	    }
 	  else if (next_only)
--- a/src/mkbuiltins
+++ b/src/mkbuiltins
@@ -65,6 +65,8 @@
 
 #endif
 
+#define XDEFUN_FILE_NAME(name)
+
 #define XDEFUN_INTERNAL(name, args_name, nargout_name, doc) \
   extern DECLARE_FUN (name, args_name, nargout_name); \
   install_builtin_function (F ## name, #name, doc); \
--- a/src/mkgendoc
+++ b/src/mkgendoc
@@ -43,42 +43,47 @@
 #include <iostream>
 #include <string>
 
+#define XDEFUN_FILE_NAME(name) \
+  std::string file_name = name;
+
 #define XDEFUN_DLD_INTERNAL(name, args_name, nargout_name, doc) \
-  print_doc_string (#name, doc);
+  print_doc_string (#name, file_name, doc);
 
 #define XDEFUNX_DLD_INTERNAL(name, fname, args_name, nargout_name, doc) \
-  print_doc_string (name, doc);
+  print_doc_string (name, file_name, doc);
 
 #define XDEFUN_INTERNAL(name, args_name, nargout_name, doc) \
-  print_doc_string (#name, doc);
+  print_doc_string (#name, file_name, doc);
 
 #define XDEFCONSTFUN_INTERNAL(name, args_name, nargout_name, doc) \
-  print_doc_string (#name, doc);
+  print_doc_string (#name, file_name, doc);
 
 #define XDEFUNX_INTERNAL(name, fname, args_name, nargout_name, doc) \
-  print_doc_string (name, doc);
+  print_doc_string (name, file_name, doc);
 
 #define XDEFALIAS_INTERNAL(alias, name)
 
 #define XDEFVAR_INTERNAL(name, sname, defn, protect, chg_fcn, doc) \
-  print_doc_string (#name, doc);
+  print_doc_string (#name, file_name, doc);
 
 #define XDEFCONST_INTERNAL(name, defn, doc) \
-  print_doc_string (#name, doc);
+  print_doc_string (#name, file_name, doc);
 
 static void
-print_doc_string (const std::string& name, const std::string& doc)
+print_doc_string (const std::string& name, const std::string& file_name,
+		  const std::string& doc)
 {
   std::cout << "";
 
   size_t len = name.length ();
 
   if (name[0] == '"' && name[len-1] == '"')
-    std::cout << name.substr (1, len-2);
+    std::cout << name.substr (1, len-2) << "\n";
   else
-    std::cout << name;
+    std::cout << name << "\n";
 
-  std::cout << "\n" << doc << "\n";
+  std::cout << "@c " << file_name << "\n"
+	    << doc << "\n";
 }
 
 EOF
--- a/src/oct-stream.cc
+++ b/src/oct-stream.cc
@@ -3595,17 +3595,36 @@
 	{
 	  std::ostream& os = *osp;
 
-	  // It seems that Matlab writes zeros instead of actually
-	  // seeking.  Hmm...
-
 	  if (skip != 0 && (i % block_size) == 0)
 	    {
-	      // FIXME -- probably should try to write larger
-	      // blocks...
-
-	      unsigned char zero = 0;
-	      for (int j = 0; j < skip; j++)
-		os.write (reinterpret_cast<const char *> (&zero), 1);
+	      // Seek to skip when inside bounds of existing file.
+	      // Otherwise, write NUL to skip.
+
+	      long orig_pos = tell ();
+
+	      seek (0, SEEK_END);
+
+	      long eof_pos = tell ();
+
+	      // Is it possible for this to fail to return us to the
+	      // original position?
+	      seek (orig_pos, SEEK_SET);
+
+	      long remaining = eof_pos - orig_pos;
+
+	      if (remaining < skip)
+		{
+		  seek (0, SEEK_END);
+
+		  // FIXME -- probably should try to write larger
+		  // blocks...
+
+		  unsigned char zero = 0;
+		  for (octave_idx_type j = 0; j < skip - remaining; j++)
+		    os.write (reinterpret_cast<const char *> (&zero), 1);
+		}
+	      else
+		seek (skip, SEEK_CUR);
 	    }
 
 	  if (os)
--- a/src/octave.cc
+++ b/src/octave.cc
@@ -367,6 +367,18 @@
   unwind_protect::run_frame ("execute_startup_files");
 }
 
+static void
+unmark_forced_vars (void *arg)
+{
+  // Unmark any symbols that may have been tagged as local variables
+  // while parsing (for example, by force_local_variable in lex.l).
+
+  symbol_table::scope_id *pscope = static_cast <symbol_table::scope_id *> (arg);
+
+  if (pscope)
+    symbol_table::unmark_forced_variables (*pscope);
+}
+
 static int
 execute_eval_option_code (const std::string& code)
 {
@@ -386,6 +398,11 @@
 
   unwind_protect_bool (interactive);
 
+  // Do this with an unwind-protect cleanup function so that the
+  // forced variables will be unmarked in the event of an interrupt.
+  symbol_table::scope_id scope = symbol_table::top_scope ();
+  unwind_protect::add (unmark_forced_vars, &scope);
+
   interactive = false;
 
   int parse_status = 0;
--- a/src/ov-base-diag.cc
+++ b/src/ov-base-diag.cc
@@ -27,6 +27,7 @@
 #include <iostream>
 
 #include "mach-info.h"
+#include "lo-ieee.h"
 
 #include "ov-base.h"
 #include "ov-base-mat.h"
--- a/src/ov-base-mat.h
+++ b/src/ov-base-mat.h
@@ -98,6 +98,8 @@
 
   dim_vector dims (void) const { return matrix.dims (); }
 
+  octave_idx_type numel (void) const { return matrix.numel (); }
+
   octave_idx_type nnz (void) const { return matrix.nnz (); }
 
   octave_value reshape (const dim_vector& new_dims) const
--- a/src/ov-base.h
+++ b/src/ov-base.h
@@ -450,6 +450,13 @@
 
   virtual string_vector map_keys (void) const;
 
+  virtual string_vector parent_class_names (void) const
+    { return string_vector (); }
+
+  // FIXME -- should this warn if called for a non-class type?
+  virtual octave_base_value *find_parent_class (const std::string&)
+    { return 0; }
+
   virtual octave_function *function_value (bool silent = false);
 
   virtual const octave_function *function_value (bool silent = false) const;
--- a/src/ov-class.cc
+++ b/src/ov-class.cc
@@ -62,6 +62,85 @@
     (octave_class::t_name, "<unknown>", octave_value (new octave_class ()));
 }
 
+octave_class::octave_class (const Octave_map& m, const std::string& id, 
+			    const octave_value_list& parents)
+  : octave_base_value (), map (m), c_name (id)
+{
+  octave_idx_type n = parents.length ();
+
+  for (octave_idx_type idx = 0; idx < n; idx++)
+    {
+      octave_value parent = parents(idx);
+
+      if (! parent.is_object ())
+	error ("parents must be objects");
+      else
+	{
+	  std::string cnm = parent.class_name ();
+
+	  parent_list.push_back (cnm);
+
+	  map.assign (cnm, parent);
+	}
+    }
+
+  load_path::add_to_parent_map (id, parent_list);
+}
+
+octave_base_value *
+octave_class::find_parent_class (const std::string& parent_class_name)
+{
+  octave_base_value* retval = 0;
+  std::string dbg_clsnm = class_name ();
+
+  if (parent_class_name == class_name ())
+    retval = this;
+  else
+    {
+      // Search in the list of immediate parents first, then in the
+      // ancestor tree.
+
+      std::list<std::string>::iterator 
+	p = find (parent_list.begin (), parent_list.end (), parent_class_name);
+
+      if (p != parent_list.end ())
+	{
+	  Octave_map::const_iterator pmap = map.seek (parent_class_name);
+
+	  if (pmap != map.end ())
+	    {
+	      const Cell& tmp = pmap->second;
+
+	      octave_value vtmp = tmp(0);
+
+	      retval = vtmp.internal_rep ();
+	    }
+	}
+      else
+	{
+	  for (std::list<std::string>::iterator pit = parent_list.begin ();
+	       pit != parent_list.end ();
+	       pit++)
+	    {
+	      Octave_map::const_iterator smap = map.seek (*pit);
+
+	      const Cell& tmp = smap->second;
+
+	      octave_value vtmp = tmp(0);
+
+	      octave_base_value *obvp = vtmp.internal_rep ();
+
+	      retval = obvp->find_parent_class (parent_class_name);
+
+	      if (retval)
+		break;
+	    }
+	}
+    }
+
+  return retval;
+}
+
 Cell
 octave_class::dotref (const octave_value_list& idx)
 {
@@ -69,12 +148,27 @@
 
   assert (idx.length () == 1);
 
+  // FIXME -- Is there a "proper" way to do this?
+  octave_function* fcn = octave_call_stack::current ();
+  std::string my_dir = fcn->dir_name ();
+  int ipos = my_dir.find_last_of ("@");
+  std::string method_class = my_dir.substr (ipos+1);
+
+  // Find the class in which this method resides before attempting to access
+  // the requested field.
+
+  octave_base_value *obvp = find_parent_class (method_class);
+
+  Octave_map my_map;
+
+  my_map = obvp ? obvp->map_value () : map;
+
   std::string nm = idx(0).string_value ();
 
-  Octave_map::const_iterator p = map.seek (nm);
+  Octave_map::const_iterator p = my_map.seek (nm);
 
-  if (p != map.end ())
-    retval = map.contents (p);
+  if (p != my_map.end ())
+    retval = my_map.contents (p);
   else
     error ("class has no member `%s'", nm.c_str ());
 
@@ -1168,18 +1262,22 @@
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} class (@var{expr})\n\
 @deftypefnx {Built-in Function} {} class (@var{s}, @var{id})\n\
-\n\
-Return the class of the expression @var{expr}, as a string or\n\
-create a class object from the structure @var{s} with name @var{id}.\n\
+@deftypefnx {Built-in Function} {} class (@var{s}, @var{id}, @var{p}, @dots{})\n\
+Return the class of the expression @var{expr} or create a class with\n\
+fields from structure @var{s} and name (string) @var{id}.  Additional\n\
+arguments name a list of parent classes from which the new class is\n\
+derived.\n\
 @end deftypefn")
 {
   octave_value retval;
 
   int nargin = args.length ();
 
-  if (nargin == 1)
+  if (nargin == 0)
+    print_usage ();
+  else if (nargin == 1)
     retval = args(0).class_name ();
-  else if (nargin == 2)
+  else
     {
       Octave_map m = args(0).map_value ();
 
@@ -1192,7 +1290,16 @@
 	      octave_function *fcn = octave_call_stack::caller ();
 
 	      if (fcn && fcn->is_class_constructor ())
-		retval = octave_value (new octave_class (m, id));
+		{
+		  if (nargin == 2)
+		    retval = octave_value (new octave_class (m, id));
+		  else
+		    {
+		      octave_value_list parents = args.slice (2, nargin-2);
+
+		      retval = octave_value (new octave_class (m, id, parents));
+		    }
+		}
 	      else
 		error ("class: invalid call from outside class constructor");
 	    }
@@ -1202,13 +1309,31 @@
       else
 	error ("class: expecting structure as first argument");
     }
+
+  return retval;
+}
+
+DEFUN (__parent_classes__, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {} __parent_classes__ (@var{x})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value retval = Cell ();
+
+  if (args.length () == 1)
+    {
+      octave_value arg = args(0);
+
+      if (arg.is_object ())
+	retval = Cell (arg.parent_class_names ());
+    }
   else
     print_usage ();
 
   return retval;
 }
 
-
 DEFUN (isobject, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} isobject (@var{x})\n\
--- a/src/ov-class.h
+++ b/src/ov-class.h
@@ -55,7 +55,11 @@
     : octave_base_value (), map (m), c_name (id) { }
 
   octave_class (const octave_class& s)
-    : octave_base_value (s), map (s.map), c_name (s.c_name) { }
+    : octave_base_value (s), map (s.map), c_name (s.c_name),
+      parent_list (s.parent_list) { }
+
+  octave_class (const Octave_map& m, const std::string& id, 
+                const octave_value_list& parents);
 
   ~octave_class (void) { }
 
@@ -118,6 +122,11 @@
 
   string_vector map_keys (void) const;
 
+  string_vector parent_class_names (void) const
+    { return string_vector (parent_list); }
+
+  octave_base_value *find_parent_class (const std::string&);
+
   void print (std::ostream& os, bool pr_as_read_syntax = false) const;
 
   void print_raw (std::ostream& os, bool pr_as_read_syntax = false) const;
@@ -165,6 +174,7 @@
 
   static const std::string t_name;
   std::string c_name;
+  std::list<std::string> parent_list;
 
   bool in_class_method (void) const;
 };
--- a/src/ov-cx-mat.cc
+++ b/src/ov-cx-mat.cc
@@ -353,20 +353,15 @@
 		{
 		  ComplexNDArray tmp(dv);
 
-		  if (tmp.is_empty ())
-		    matrix = tmp;
-		  else
-		    {
-		      is >> tmp;
+                  is >> tmp;
 
-		      if (is)
-			matrix = tmp;
-		      else
-			{
-			  error ("load: failed to load matrix constant");
-			  success = false;
-			}
-		    }
+                  if (is)
+                    matrix = tmp;
+                  else
+                    {
+                      error ("load: failed to load matrix constant");
+                      success = false;
+                    }
 		}
 	      else
 		{
@@ -856,6 +851,24 @@
   return ::imag (matrix);
 }
 
+octave_value
+octave_complex_matrix::isnan (void) const
+{
+  return matrix.isnan ();
+}
+
+octave_value
+octave_complex_matrix::isinf (void) const
+{
+  return matrix.isinf ();
+}
+
+octave_value
+octave_complex_matrix::finite (void) const
+{
+  return matrix.isfinite ();
+}
+
 DARRAY_MAPPER (erf, NDArray::dmapper, ::erf)
 DARRAY_MAPPER (erfc, NDArray::dmapper, ::erfc)
 DARRAY_MAPPER (gamma, NDArray::dmapper, xgamma)
@@ -888,10 +901,7 @@
 ARRAY_MAPPER (sqrt, ComplexNDArray::cmapper, std::sqrt)
 ARRAY_MAPPER (tan, ComplexNDArray::cmapper, std::tan)
 ARRAY_MAPPER (tanh, ComplexNDArray::cmapper, std::tanh)
-ARRAY_MAPPER (finite, ComplexNDArray::bmapper, xfinite)
-ARRAY_MAPPER (isinf, ComplexNDArray::bmapper, xisinf)
 ARRAY_MAPPER (isna, ComplexNDArray::bmapper, octave_is_NA)
-ARRAY_MAPPER (isnan, ComplexNDArray::bmapper, xisnan)
 
 /*
 ;;; Local Variables: ***
--- a/src/ov-flt-cx-mat.cc
+++ b/src/ov-flt-cx-mat.cc
@@ -342,20 +342,15 @@
 		{
 		  FloatComplexNDArray tmp(dv);
 
-		  if (tmp.is_empty ())
-		    matrix = tmp;
-		  else
-		    {
-		      is >> tmp;
+                  is >> tmp;
 
-		      if (is)
-			matrix = tmp;
-		      else
-			{
-			  error ("load: failed to load matrix constant");
-			  success = false;
-			}
-		    }
+                  if (is)
+                    matrix = tmp;
+                  else
+                    {
+                      error ("load: failed to load matrix constant");
+                      success = false;
+                    }
 		}
 	      else
 		{
@@ -823,6 +818,24 @@
   return ::imag (matrix);
 }
 
+octave_value
+octave_float_complex_matrix::isnan (void) const
+{
+  return matrix.isnan ();
+}
+
+octave_value
+octave_float_complex_matrix::isinf (void) const
+{
+  return matrix.isinf ();
+}
+
+octave_value
+octave_float_complex_matrix::finite (void) const
+{
+  return matrix.isfinite ();
+}
+
 DARRAY_MAPPER (erf, FloatNDArray::dmapper, ::erff)
 DARRAY_MAPPER (erfc, FloatNDArray::dmapper, ::erfcf)
 DARRAY_MAPPER (gamma, FloatNDArray::dmapper, xgamma)
@@ -855,10 +868,7 @@
 ARRAY_MAPPER (sqrt, FloatComplexNDArray::cmapper, std::sqrt)
 ARRAY_MAPPER (tan, FloatComplexNDArray::cmapper, std::tan)
 ARRAY_MAPPER (tanh, FloatComplexNDArray::cmapper, std::tanh)
-ARRAY_MAPPER (finite, FloatComplexNDArray::bmapper, xfinite)
-ARRAY_MAPPER (isinf, FloatComplexNDArray::bmapper, xisinf)
 ARRAY_MAPPER (isna, FloatComplexNDArray::bmapper, octave_is_NA)
-ARRAY_MAPPER (isnan, FloatComplexNDArray::bmapper, xisnan)
 
 /*
 ;;; Local Variables: ***
--- a/src/ov-flt-re-mat.cc
+++ b/src/ov-flt-re-mat.cc
@@ -372,20 +372,15 @@
 		{
 		  FloatNDArray tmp(dv);
 
-		  if (tmp.is_empty ())
-		    matrix = tmp;
-		  else
-		    {
-		      is >> tmp;
+                  is >> tmp;
 
-		      if (is)
-			matrix = tmp;
-		      else
-			{
-			  error ("load: failed to load matrix constant");
-			  success = false;
-			}
-		    }
+                  if (is)
+                    matrix = tmp;
+                  else
+                    {
+                      error ("load: failed to load matrix constant");
+                      success = false;
+                    }
 		}
 	      else
 		{
@@ -766,6 +761,24 @@
   return FloatNDArray (matrix.dims (), 0.0);
 }
 
+octave_value
+octave_float_matrix::isnan (void) const
+{
+  return matrix.isnan ();
+}
+
+octave_value
+octave_float_matrix::isinf (void) const
+{
+  return matrix.isinf ();
+}
+
+octave_value
+octave_float_matrix::finite (void) const
+{
+  return matrix.isfinite ();
+}
+
 ARRAY_MAPPER (erf, FloatNDArray::dmapper, ::erff)
 ARRAY_MAPPER (erfc, FloatNDArray::dmapper, ::erfcf)
 ARRAY_MAPPER (gamma, FloatNDArray::dmapper, xgamma)
@@ -797,10 +810,7 @@
 CD_ARRAY_MAPPER (sqrt, ::sqrtf, std::sqrt, 0.0, octave_Float_Inf)
 ARRAY_MAPPER (tan, FloatNDArray::dmapper, ::tanf)
 ARRAY_MAPPER (tanh, FloatNDArray::dmapper, ::tanhf)
-ARRAY_MAPPER (finite, FloatNDArray::bmapper, xfinite)
-ARRAY_MAPPER (isinf, FloatNDArray::bmapper, xisinf)
 ARRAY_MAPPER (isna, FloatNDArray::bmapper, octave_is_NA)
-ARRAY_MAPPER (isnan, FloatNDArray::bmapper, xisnan)
 
 DEFUN (single, args, ,
   "-*- texinfo -*-\n\
--- a/src/ov-re-mat.cc
+++ b/src/ov-re-mat.cc
@@ -380,20 +380,15 @@
 		{
 		  NDArray tmp(dv);
 
-		  if (tmp.is_empty ())
-		    matrix = tmp;
-		  else
-		    {
-		      is >> tmp;
+                  is >> tmp;
 
-		      if (is)
-			matrix = tmp;
-		      else
-			{
-			  error ("load: failed to load matrix constant");
-			  success = false;
-			}
-		    }
+                  if (is)
+                    matrix = tmp;
+                  else
+                    {
+                      error ("load: failed to load matrix constant");
+                      success = false;
+                    }
 		}
 	      else
 		{
@@ -794,6 +789,24 @@
   return NDArray (matrix.dims (), 0.0);
 }
 
+octave_value
+octave_matrix::isnan (void) const
+{
+  return matrix.isnan ();
+}
+
+octave_value
+octave_matrix::isinf (void) const
+{
+  return matrix.isinf ();
+}
+
+octave_value
+octave_matrix::finite (void) const
+{
+  return matrix.isfinite ();
+}
+
 ARRAY_MAPPER (erf, NDArray::dmapper, ::erf)
 ARRAY_MAPPER (erfc, NDArray::dmapper, ::erfc)
 ARRAY_MAPPER (gamma, NDArray::dmapper, xgamma)
@@ -825,10 +838,7 @@
 CD_ARRAY_MAPPER (sqrt, ::sqrt, std::sqrt, 0.0, octave_Inf)
 ARRAY_MAPPER (tan, NDArray::dmapper, ::tan)
 ARRAY_MAPPER (tanh, NDArray::dmapper, ::tanh)
-ARRAY_MAPPER (finite, NDArray::bmapper, xfinite)
-ARRAY_MAPPER (isinf, NDArray::bmapper, xisinf)
 ARRAY_MAPPER (isna, NDArray::bmapper, octave_is_NA)
-ARRAY_MAPPER (isnan, NDArray::bmapper, xisnan)
 
 DEFUN (double, args, ,
   "-*- texinfo -*-\n\
--- a/src/ov.cc
+++ b/src/ov.cc
@@ -474,11 +474,6 @@
   return retval;
 }
 
-octave_value::octave_value (void)
-  : rep (new octave_base_value ())
-{
-}
-
 octave_value::octave_value (short int i)
   : rep (new octave_scalar (i))
 {
--- a/src/ov.h
+++ b/src/ov.h
@@ -156,7 +156,13 @@
 
   enum magic_colon { magic_colon_t };
 
-  octave_value (void);
+  octave_value (void)
+    {
+      static octave_base_value nil_rep;
+      rep = &nil_rep;
+      rep->count++;
+    }
+
   octave_value (short int i);
   octave_value (unsigned short int i);
   octave_value (int i);
@@ -824,6 +830,13 @@
   string_vector map_keys (void) const
     { return rep->map_keys (); }
 
+  string_vector parent_class_names (void) const
+    { return rep->parent_class_names (); }
+
+  octave_base_value *
+  find_parent_class (const std::string& parent_class_name)
+    { return rep->find_parent_class (parent_class_name); }
+
   octave_function *function_value (bool silent = false);
 
   const octave_function *function_value (bool silent = false) const;
--- a/src/parse.y
+++ b/src/parse.y
@@ -2528,7 +2528,10 @@
 
   curr_fcn_ptr = script;
 
-  symbol_table::clear_forced_variables ();
+  // Unmark any symbols that may have been tagged as local variables
+  // while parsing (for example, by force_local_variable in lex.l).
+
+  symbol_table::unmark_forced_variables ();
 }
 
 // Begin defining a function.
@@ -2708,10 +2711,11 @@
 	  retval = new tree_function_def (fcn);
 	}
 
-      // Clear any local variables that may have been added while
-      // parsing (for example, by force_local_variable in lex.l). 
-
-      symbol_table::clear_forced_variables (fcn->scope ());
+      // Unmark any symbols that may have been tagged as local
+      // variables while parsing (for example, by force_local_variable
+      // in lex.l).
+
+      symbol_table::unmark_forced_variables (fcn->scope ());
     }
 
   return retval;
--- a/src/pt-eval.cc
+++ b/src/pt-eval.cc
@@ -339,8 +339,8 @@
 
 	DO_SIMPLE_FOR_LOOP_ONCE (rhs);
       }
-    else if (rhs.is_matrix_type () 
-             || rhs.is_cell () || rhs.is_string ())
+    else if (rhs.is_matrix_type () || rhs.is_cell () || rhs.is_string ()
+             || rhs.is_map ())
       {
         // A matrix or cell is reshaped to 2 dimensions and iterated by
         // columns.
@@ -384,27 +384,6 @@
               }
           }
       }
-    else if (rhs.is_map ())
-      {
-	Octave_map tmp_val (rhs.map_value ());
-
-	bool quit = false;
-
-	for (Octave_map::iterator p = tmp_val.begin ();
-	     p != tmp_val.end ();
-	     p++)
-	  {
-	    Cell val_lst = tmp_val.contents (p);
-
-	    octave_value val
-	      = (val_lst.length () == 1) ? val_lst(0) : octave_value (val_lst);
-
-	    DO_SIMPLE_FOR_LOOP_ONCE (val);
-
-	    if (quit)
-	      break;
-	  }
-      }
     else
       {
 	::error ("invalid type in for loop expression near line %d, column %d",
@@ -456,15 +435,15 @@
 
       octave_lvalue key_ref = elt->lvalue ();
 
-      Octave_map tmp_val (rhs.map_value ());
+      const Octave_map tmp_val (rhs.map_value ());
 
       tree_statement_list *loop_body = cmd.body ();
 
-      for (Octave_map::iterator q = tmp_val.begin (); q != tmp_val.end (); q++)
+      for (Octave_map::const_iterator q = tmp_val.begin (); q != tmp_val.end (); q++)
 	{
 	  octave_value key = tmp_val.key (q);
 
-	  Cell val_lst = tmp_val.contents (q);
+	  const Cell val_lst = tmp_val.contents (q);
 
 	  octave_idx_type n = tmp_val.numel ();
 
--- a/src/sparse-xdiv.cc
+++ b/src/sparse-xdiv.cc
@@ -31,6 +31,7 @@
 #include "oct-cmplx.h"
 #include "quit.h"
 #include "error.h"
+#include "lo-ieee.h"
 
 #include "dSparse.h"
 #include "dDiagMatrix.h"
--- a/src/symtab.h
+++ b/src/symtab.h
@@ -204,10 +204,7 @@
 	octave_value& val = varref (context);
 
 	if (! val.is_defined ())
-	  {
-	    val = Matrix ();
-	    mark_forced ();
-	  }
+	  mark_forced ();
       }
 
       octave_value& varref (context_id context)
@@ -293,15 +290,6 @@
 	  }
       }
 
-      void clear_forced (void)
-      {
-	if (is_forced ())
-	  {
-	    varref (xcurrent_context) = octave_value ();
-	    unmark_forced ();
-	  }
-      }
-
       bool is_defined (context_id context) const
       {
 	return varval (context).is_defined ();
@@ -309,7 +297,7 @@
 
       bool is_variable (context_id context) const
       {
-	return (storage_class != local || is_defined (context));
+	return (storage_class != local || is_defined (context) || is_forced ());
       }
 
       bool is_local (void) const { return storage_class & local; }
@@ -457,8 +445,6 @@
 
     void clear (void) { rep->clear (); }
 
-    void clear_forced (void) { rep->clear_forced (); }
-    
     bool is_defined (context_id context = xcurrent_context) const
     {
       return rep->is_defined (context);
@@ -1353,12 +1339,12 @@
       inst->do_clear_variables ();
   }
 
-  static void clear_forced_variables (scope_id scope = xcurrent_scope)
+  static void unmark_forced_variables (scope_id scope = xcurrent_scope)
   {
     symbol_table *inst = get_instance (scope);
 
     if (inst)
-      inst->do_clear_forced_variables ();
+      inst->do_unmark_forced_variables ();
   }
 
   // For unwind_protect.
@@ -2133,10 +2119,10 @@
       p->second.clear ();
   }
 
-  void do_clear_forced_variables (void)
+  void do_unmark_forced_variables (void)
   {
     for (table_iterator p = table.begin (); p != table.end (); p++)
-      p->second.clear_forced ();
+      p->second.unmark_forced ();
   }
 
   void do_clear_global (const std::string& name)