# HG changeset patch
# User Max Brister <max@2bass.com>
# Date 1340653221 18000
# Node ID 0f156affb036037bda38c89de140c9132d6253cd
# Parent  bbeef7b8ea2ef49e1d5e16cff33874e0cfeef677# Parent  3c5e6971064cc4bd5e12898f2a679bfe8dc108f1
maint: periodic merge of default to jit

diff --git a/examples/embedded.cc b/examples/embedded.cc
--- a/examples/embedded.cc
+++ b/examples/embedded.cc
@@ -13,30 +13,21 @@
   octave_main (2, argv.c_str_vec(), 1);
 
   octave_idx_type n = 2;
-  Matrix a_matrix = Matrix (1, 2);
+  octave_value_list in;
 
-  std::cout << "GCD of [";
-  for (octave_idx_type i = 0; i < n; i++)
-    {
-      a_matrix (i) = 5 * (i + 1);
-      if (i != 0)
-        std::cout << ", " << 5 * (i + 2);
-      else
-        std::cout << 5 * (i + 2);
-    }
-  std::cout << "] is ";
-
-  octave_value_list in = octave_value (a_matrix);
+  for (octave_idx_type i = 0; i < n; i++)  
+    in(i) = octave_value (5 * (i + 1));
+  
   octave_value_list out = feval ("gcd", in, 1);
 
+  
   if (!error_state && out.length () > 0)
-    {
-      a_matrix = out(0).matrix_value ();
-      if (a_matrix.numel () == 1)
-        std::cout << a_matrix(0) << "\n";
-      else
-        std::cout << "invalid\n";
-    }
+    std::cout << "GCD of [" 
+              << in(0).int_value () 
+              << ", " 
+              << in(1).int_value ()
+              << "] is " << out(0).int_value () 
+              << std::endl;
   else
     std::cout << "invalid\n";
 
diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc
--- a/liboctave/Sparse.cc
+++ b/liboctave/Sparse.cc
@@ -1240,15 +1240,31 @@
         gripe_del_index_out_of_range (false, idx_j.extent (nc), nc);
       else if (idx_j.is_cont_range (nc, lb, ub))
         {
-          const Sparse<T> tmp = *this;
-          octave_idx_type lbi = tmp.cidx(lb), ubi = tmp.cidx(ub), new_nz = nz - (ubi - lbi);
-          *this = Sparse<T> (nr, nc - (ub - lb), new_nz);
-          copy_or_memcpy (lbi, tmp.data (), data ());
-          copy_or_memcpy (lbi, tmp.ridx (), ridx ());
-          copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi);
-          copy_or_memcpy (nz - ubi, tmp.ridx () + ubi, xridx () + lbi);
-          copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1);
-          mx_inline_sub (nc - ub, xcidx () + 1, tmp.cidx () + ub + 1, ubi - lbi);
+          if (lb == 0 && ub == nc)
+            {
+              // Delete all rows and columns.
+              *this = Sparse<T> (nr, 0);
+            }
+          else if (nz == 0)
+            {
+              // No elements to preserve; adjust dimensions.
+              *this = Sparse<T> (nr, nc - (ub - lb));
+            }
+          else
+            {
+              const Sparse<T> tmp = *this;
+              octave_idx_type lbi = tmp.cidx(lb), ubi = tmp.cidx(ub),
+                new_nz = nz - (ubi - lbi);
+
+              *this = Sparse<T> (nr, nc - (ub - lb), new_nz);
+              copy_or_memcpy (lbi, tmp.data (), data ());
+              copy_or_memcpy (lbi, tmp.ridx (), ridx ());
+              copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi);
+              copy_or_memcpy (nz - ubi, tmp.ridx () + ubi, xridx () + lbi);
+              copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1);
+              mx_inline_sub (nc - ub, xcidx () + lb + 1,
+                             tmp.cidx () + ub + 1, ubi - lbi);
+            }
         }
       else
         *this = index (idx_i, idx_j.complement (nc));
@@ -1261,24 +1277,40 @@
         gripe_del_index_out_of_range (false, idx_i.extent (nr), nr);
       else if (idx_i.is_cont_range (nr, lb, ub))
         {
-          // This is more memory-efficient than the approach below.
-          const Sparse<T> tmpl = index (idx_vector (0, lb), idx_j);
-          const Sparse<T> tmpu = index (idx_vector (ub, nr), idx_j);
-          *this = Sparse<T> (nr - (ub - lb), nc, tmpl.nnz () + tmpu.nnz ());
-          for (octave_idx_type j = 0, k = 0; j < nc; j++)
+          if (lb == 0 && ub == nr)
+            {
+              // Delete all rows and columns.
+              *this = Sparse<T> (0, nc);
+            }
+          else if (nz == 0)
             {
-              for (octave_idx_type i = tmpl.cidx(j); i < tmpl.cidx(j+1); i++)
+              // No elements to preserve; adjust dimensions.
+              *this = Sparse<T> (nr - (ub - lb), nc);
+            }
+          else
+            {
+              // This is more memory-efficient than the approach below.
+              const Sparse<T> tmpl = index (idx_vector (0, lb), idx_j);
+              const Sparse<T> tmpu = index (idx_vector (ub, nr), idx_j);
+              *this = Sparse<T> (nr - (ub - lb), nc,
+                                 tmpl.nnz () + tmpu.nnz ());
+              for (octave_idx_type j = 0, k = 0; j < nc; j++)
                 {
-                  xdata(k) = tmpl.data(i);
-                  xridx(k++) = tmpl.ridx(i);
+                  for (octave_idx_type i = tmpl.cidx(j); i < tmpl.cidx(j+1);
+                       i++)
+                    {
+                      xdata(k) = tmpl.data(i);
+                      xridx(k++) = tmpl.ridx(i);
+                    }
+                  for (octave_idx_type i = tmpu.cidx(j); i < tmpu.cidx(j+1);
+                       i++)
+                    {
+                      xdata(k) = tmpu.data(i);
+                      xridx(k++) = tmpu.ridx(i) + lb;
+                    }
+
+                  xcidx(j+1) = k;
                 }
-              for (octave_idx_type i = tmpu.cidx(j); i < tmpu.cidx(j+1); i++)
-                {
-                  xdata(k) = tmpu.data(i);
-                  xridx(k++) = tmpu.ridx(i) + lb;
-                }
-
-              xcidx(j+1) = k;
             }
         }
       else
@@ -1583,6 +1615,10 @@
     {
       // It's actually vector indexing. The 1D index is specialized for that.
       retval = index (idx_i);
+
+      // If nr == 1 then the vector indexing will return a column vector!!
+      if (nr == 1)
+        retval.transpose();
     }
   else if (idx_i.is_scalar ())
     {
@@ -2736,6 +2772,18 @@
 
 %!assert (speye (3,1)(3:-1:1), sparse ([0; 0; 1]))
 
+## Test removing columns (bug #36656)
+
+%!test
+%! s = sparse (magic (5));
+%! s(:,2:4) = [];
+%! assert (s, sparse (magic (5)(:, [1,5])));
+
+%!test
+%! s = sparse([], [], [], 1, 1);
+%! s(1,:) = [];
+%! assert (s, sparse ([], [], [], 0, 1));
+
 */
 
 template <class T>
diff --git a/liboctave/Sparse.h b/liboctave/Sparse.h
--- a/liboctave/Sparse.h
+++ b/liboctave/Sparse.h
@@ -71,33 +71,34 @@
     octave_idx_type ncols;
     octave_refcount<int> count;
 
-    SparseRep (void) : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0),
-                       ncols (0), count (1) { c[0] = 0; }
+    SparseRep (void)
+      : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0),
+      ncols (0), count (1)
+      {
+        c[0] = 0;
+      }
 
-    SparseRep (octave_idx_type n) : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n),
+    SparseRep (octave_idx_type n)
+      : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n),
       ncols (n), count (1)
       {
         for (octave_idx_type i = 0; i < n + 1; i++)
           c[i] = 0;
       }
 
-    SparseRep (octave_idx_type nr, octave_idx_type nc) : d (0), r (0), c (new octave_idx_type [nc+1]), nzmx (0),
-      nrows (nr), ncols (nc), count (1)
-      {
-        for (octave_idx_type i = 0; i < nc + 1; i++)
-          c[i] = 0;
-      }
-
-    SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : d (new T [nz]),
-      r (new octave_idx_type [nz]), c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr),
+    SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz = 0)
+      : d (new T [nz]), r (new octave_idx_type [nz]),
+      c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr),
       ncols (nc), count (1)
       {
-        for (octave_idx_type i = 0; i < nc + 1; i++)
+        c[nc] = nz;
+        for (octave_idx_type i = 0; i < nc; i++)
           c[i] = 0;
       }
 
     SparseRep (const SparseRep& a)
-      : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]), c (new octave_idx_type [a.ncols + 1]),
+      : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]),
+      c (new octave_idx_type [a.ncols + 1]),
       nzmx (a.nzmx), nrows (a.nrows), ncols (a.ncols), count (1)
       {
         octave_idx_type nz = a.nnz ();
@@ -144,17 +145,17 @@
   //--------------------------------------------------------------------
 
   void make_unique (void)
-    {
-      if (rep->count > 1)
-        {
-          SparseRep *r = new SparseRep (*rep);
+  {
+    if (rep->count > 1)
+      {
+        SparseRep *r = new SparseRep (*rep);
 
-          if (--rep->count == 0)
-            delete rep;
+        if (--rep->count == 0)
+          delete rep;
 
-          rep = r;
-        }
-    }
+        rep = r;
+      }
+  }
 
 public:
 
@@ -168,18 +169,18 @@
 private:
 
   typename Sparse<T>::SparseRep *nil_rep (void) const
-    {
-      static typename Sparse<T>::SparseRep nr;
-      return &nr;
-    }
+  {
+    static typename Sparse<T>::SparseRep nr;
+    return &nr;
+  }
 
 public:
 
   Sparse (void)
     : rep (nil_rep ()), dimensions (dim_vector(0,0))
-    {
-      rep->count++;
-    }
+  {
+    rep->count++;
+  }
 
   explicit Sparse (octave_idx_type n)
     : rep (new typename Sparse<T>::SparseRep (n)),
@@ -193,7 +194,7 @@
 
   Sparse (const dim_vector& dv, octave_idx_type nz)
     : rep (new typename Sparse<T>::SparseRep (dv(0), dv(1), nz)),
-    dimensions (dv) { }
+      dimensions (dv) { }
 
   Sparse (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz)
     : rep (new typename Sparse<T>::SparseRep (nr, nc, nz)),
@@ -207,20 +208,20 @@
   template <class U>
   Sparse (const Sparse<U>& a)
     : rep (new typename Sparse<T>::SparseRep (a.rep->nrows, a.rep->ncols, a.rep->nzmx)),
-    dimensions (a.dimensions)
-    {
-      octave_idx_type nz = a.nnz ();
-      std::copy (a.rep->d, a.rep->d + nz, rep->d);
-      copy_or_memcpy (nz, a.rep->r, rep->r);
-      copy_or_memcpy (rep->ncols + 1, a.rep->c, rep->c);
-    }
+      dimensions (a.dimensions)
+  {
+    octave_idx_type nz = a.nnz ();
+    std::copy (a.rep->d, a.rep->d + nz, rep->d);
+    copy_or_memcpy (nz, a.rep->r, rep->r);
+    copy_or_memcpy (rep->ncols + 1, a.rep->c, rep->c);
+  }
 
   // No type conversion case.
   Sparse (const Sparse<T>& a)
     : rep (a.rep), dimensions (a.dimensions)
-    {
-      rep->count++;
-    }
+  {
+    rep->count++;
+  }
 
 public:
 
@@ -249,9 +250,9 @@
   // Querying the number of elements (incl. zeros) may overflow the index type,
   // so don't do it unless you really need it.
   octave_idx_type numel (void) const
-    {
-      return dimensions.safe_numel ();
-    }
+  {
+    return dimensions.safe_numel ();
+  }
 
   octave_idx_type nelem (void) const { return capacity (); }
   octave_idx_type length (void) const { return numel (); }
@@ -265,18 +266,19 @@
 
   octave_idx_type get_row_index (octave_idx_type k) { return ridx (k); }
   octave_idx_type get_col_index (octave_idx_type k)
-    {
-      octave_idx_type ret = 0;
-      while (cidx(ret+1) < k)
-        ret++;
-      return ret;
-    }
+  {
+    octave_idx_type ret = 0;
+    while (cidx(ret+1) < k)
+      ret++;
+    return ret;
+  }
 
   size_t byte_size (void) const
-    {
-      return (static_cast<size_t>(cols () + 1) * sizeof (octave_idx_type)
-              + static_cast<size_t> (capacity ()) * (sizeof (T) + sizeof (octave_idx_type)));
-    }
+  {
+    return (static_cast<size_t>(cols () + 1) * sizeof (octave_idx_type)
+            + static_cast<size_t> (capacity ())
+            * (sizeof (T) + sizeof (octave_idx_type)));
+  }
 
   dim_vector dims (void) const { return dimensions; }
 
@@ -296,145 +298,189 @@
   // No checking, even for multiple references, ever.
 
   T& xelem (octave_idx_type n)
-    {
-      octave_idx_type i = n % rows (), j = n / rows();
-      return xelem (i, j);
-    }
+  {
+    octave_idx_type i = n % rows (), j = n / rows();
+    return xelem (i, j);
+  }
 
   T xelem (octave_idx_type n) const
-    {
-      octave_idx_type i = n % rows (), j = n / rows();
-      return xelem (i, j);
-    }
+  {
+    octave_idx_type i = n % rows (), j = n / rows();
+    return xelem (i, j);
+  }
 
   T& xelem (octave_idx_type i, octave_idx_type j) { return rep->elem (i, j); }
-  T xelem (octave_idx_type i, octave_idx_type j) const { return rep->celem (i, j); }
+  T xelem (octave_idx_type i, octave_idx_type j) const
+  {
+    return rep->celem (i, j);
+  }
 
   T& xelem (const Array<octave_idx_type>& ra_idx)
-    { return xelem (compute_index (ra_idx)); }
+  { return xelem (compute_index (ra_idx)); }
 
   T xelem (const Array<octave_idx_type>& ra_idx) const
-    { return xelem (compute_index (ra_idx)); }
+  { return xelem (compute_index (ra_idx)); }
 
   // FIXME -- would be nice to fix this so that we don't
   // unnecessarily force a copy, but that is not so easy, and I see no
   // clean way to do it.
 
   T& checkelem (octave_idx_type n)
-    {
-      if (n < 0 || n >= numel ())
-        return range_error ("T& Sparse<T>::checkelem", n);
-      else
-        {
-          make_unique ();
-          return xelem (n);
-        }
-    }
+  {
+    if (n < 0 || n >= numel ())
+      return range_error ("T& Sparse<T>::checkelem", n);
+    else
+      {
+        make_unique ();
+        return xelem (n);
+      }
+  }
 
   T& checkelem (octave_idx_type i, octave_idx_type j)
-    {
-      if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
-        return range_error ("T& Sparse<T>::checkelem", i, j);
-      else
-        {
-          make_unique ();
-          return xelem (i, j);
-        }
-    }
+  {
+    if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
+      return range_error ("T& Sparse<T>::checkelem", i, j);
+    else
+      {
+        make_unique ();
+        return xelem (i, j);
+      }
+  }
 
   T& checkelem (const Array<octave_idx_type>& ra_idx)
-    {
-      octave_idx_type i = compute_index (ra_idx);
+  {
+    octave_idx_type i = compute_index (ra_idx);
 
-      if (i < 0)
-        return range_error ("T& Sparse<T>::checkelem", ra_idx);
-      else
-        return elem (i);
-    }
+    if (i < 0)
+      return range_error ("T& Sparse<T>::checkelem", ra_idx);
+    else
+      return elem (i);
+  }
 
   T& elem (octave_idx_type n)
-    {
-      make_unique ();
-      return xelem (n);
-    }
+  {
+    make_unique ();
+    return xelem (n);
+  }
 
   T& elem (octave_idx_type i, octave_idx_type j)
-    {
-      make_unique ();
-      return xelem (i, j);
-    }
+  {
+    make_unique ();
+    return xelem (i, j);
+  }
 
   T& elem (const Array<octave_idx_type>& ra_idx)
-    { return Sparse<T>::elem (compute_index (ra_idx)); }
+  { return Sparse<T>::elem (compute_index (ra_idx)); }
 
 #if defined (BOUNDS_CHECKING)
-  T& operator () (octave_idx_type n) { return checkelem (n); }
-  T& operator () (octave_idx_type i, octave_idx_type j) { return checkelem (i, j); }
-  T& operator () (const Array<octave_idx_type>& ra_idx) { return checkelem (ra_idx); }
+  T& operator () (octave_idx_type n)
+  {
+    return checkelem (n);
+  }
+
+  T& operator () (octave_idx_type i, octave_idx_type j)
+  {
+    return checkelem (i, j);
+  }
+
+  T& operator () (const Array<octave_idx_type>& ra_idx)
+  {
+    return checkelem (ra_idx);
+  }
+
 #else
-  T& operator () (octave_idx_type n) { return elem (n); }
-  T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); }
-  T& operator () (const Array<octave_idx_type>& ra_idx) { return elem (ra_idx); }
+  T& operator () (octave_idx_type n)
+  {
+    return elem (n);
+  }
+
+  T& operator () (octave_idx_type i, octave_idx_type j)
+  {
+    return elem (i, j);
+  }
+
+  T& operator () (const Array<octave_idx_type>& ra_idx)
+  {
+    return elem (ra_idx);
+  }
+
 #endif
 
   T checkelem (octave_idx_type n) const
-    {
-      if (n < 0 || n >= numel ())
-        return range_error ("T Sparse<T>::checkelem", n);
-      else
-        return xelem (n);
-    }
+  {
+    if (n < 0 || n >= numel ())
+      return range_error ("T Sparse<T>::checkelem", n);
+    else
+      return xelem (n);
+  }
 
   T checkelem (octave_idx_type i, octave_idx_type j) const
-    {
-      if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
-        return range_error ("T Sparse<T>::checkelem", i, j);
-      else
-        return xelem (i, j);
-    }
+  {
+    if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
+      return range_error ("T Sparse<T>::checkelem", i, j);
+    else
+      return xelem (i, j);
+  }
 
   T checkelem (const Array<octave_idx_type>& ra_idx) const
-    {
-      octave_idx_type i = compute_index (ra_idx);
+  {
+    octave_idx_type i = compute_index (ra_idx);
 
-      if (i < 0)
-        return range_error ("T Sparse<T>::checkelem", ra_idx);
-      else
-        return Sparse<T>::elem (i);
-    }
+    if (i < 0)
+      return range_error ("T Sparse<T>::checkelem", ra_idx);
+    else
+      return Sparse<T>::elem (i);
+  }
 
   T elem (octave_idx_type n) const { return xelem (n); }
 
   T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); }
 
   T elem (const Array<octave_idx_type>& ra_idx) const
-    { return Sparse<T>::elem (compute_index (ra_idx)); }
+  { return Sparse<T>::elem (compute_index (ra_idx)); }
 
 #if defined (BOUNDS_CHECKING)
   T operator () (octave_idx_type n) const { return checkelem (n); }
-  T operator () (octave_idx_type i, octave_idx_type j) const { return checkelem (i, j); }
-  T operator () (const Array<octave_idx_type>& ra_idx) const { return checkelem (ra_idx); }
+  T operator () (octave_idx_type i, octave_idx_type j) const
+  {
+    return checkelem (i, j);
+  }
+
+  T operator () (const Array<octave_idx_type>& ra_idx) const
+  {
+    return checkelem (ra_idx);
+  }
+
 #else
   T operator () (octave_idx_type n) const { return elem (n); }
-  T operator () (octave_idx_type i, octave_idx_type j) const { return elem (i, j); }
-  T operator () (const Array<octave_idx_type>& ra_idx) const { return elem (ra_idx); }
+  T operator () (octave_idx_type i, octave_idx_type j) const
+  {
+    return elem (i, j);
+  }
+
+  T operator () (const Array<octave_idx_type>& ra_idx) const
+  {
+    return elem (ra_idx);
+  }
 #endif
 
   Sparse<T> maybe_compress (bool remove_zeros = false)
-    {
-      if (remove_zeros)
-        make_unique (); // Needs to unshare because elements are removed.
+  {
+    if (remove_zeros)
+      make_unique (); // Needs to unshare because elements are removed.
 
-      rep->maybe_compress (remove_zeros);
-      return (*this);
-    }
+    rep->maybe_compress (remove_zeros);
+    return (*this);
+  }
 
   Sparse<T> reshape (const dim_vector& new_dims) const;
 
   Sparse<T> permute (const Array<octave_idx_type>& vec, bool inv = false) const;
 
   Sparse<T> ipermute (const Array<octave_idx_type>& vec) const
-    { return permute (vec, true); }
+  {
+    return permute (vec, true);
+  }
 
   void resize1 (octave_idx_type n);
 
@@ -443,11 +489,11 @@
   void resize (const dim_vector& dv);
 
   void change_capacity (octave_idx_type nz)
-    {
-      if (nz < nnz ())
-        make_unique (); // Unshare now because elements will be truncated.
-      rep->change_length (nz);
-    }
+  {
+    if (nz < nnz ())
+      make_unique (); // Unshare now because elements will be truncated.
+    rep->change_length (nz);
+  }
 
   Sparse<T>& insert (const Sparse<T>& a, octave_idx_type r, octave_idx_type c);
   Sparse<T>& insert (const Sparse<T>& a, const Array<octave_idx_type>& idx);
@@ -468,7 +514,11 @@
   T* data (void) const { return rep->d; }
 
   octave_idx_type* ridx (void) { make_unique (); return rep->r; }
-  octave_idx_type& ridx (octave_idx_type i) { make_unique (); return rep->ridx (i); }
+  octave_idx_type& ridx (octave_idx_type i)
+  {
+    make_unique (); return rep->ridx (i);
+  }
+
   octave_idx_type* xridx (void) { return rep->r; }
   octave_idx_type& xridx (octave_idx_type i) { return rep->ridx (i); }
 
@@ -477,7 +527,11 @@
   octave_idx_type* ridx (void) const { return rep->r; }
 
   octave_idx_type* cidx (void) { make_unique (); return rep->c; }
-  octave_idx_type& cidx (octave_idx_type i) { make_unique (); return rep->cidx (i); }
+  octave_idx_type& cidx (octave_idx_type i)
+  {
+    make_unique (); return rep->cidx (i);
+  }
+
   octave_idx_type* xcidx (void) { return rep->c; }
   octave_idx_type& xcidx (octave_idx_type i) { return rep->cidx (i); }
 
@@ -495,7 +549,8 @@
 
   Sparse<T> index (const idx_vector& i, bool resize_ok = false) const;
 
-  Sparse<T> index (const idx_vector& i, const idx_vector& j, bool resize_ok = false) const;
+  Sparse<T> index (const idx_vector& i, const idx_vector& j,
+                   bool resize_ok = false) const;
 
   void assign (const idx_vector& i, const Sparse<T>& rhs);
 
@@ -507,13 +562,19 @@
   // You should not use them anywhere else.
   void *mex_get_data (void) const { return const_cast<T *> (data ()); }
 
-  octave_idx_type *mex_get_ir (void) const { return const_cast<octave_idx_type *> (ridx ()); }
+  octave_idx_type *mex_get_ir (void) const
+  {
+    return const_cast<octave_idx_type *> (ridx ());
+  }
 
-  octave_idx_type *mex_get_jc (void) const { return const_cast<octave_idx_type *> (cidx ()); }
+  octave_idx_type *mex_get_jc (void) const
+  {
+    return const_cast<octave_idx_type *> (cidx ());
+  }
 
   Sparse<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const;
   Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
-                 sortmode mode = ASCENDING) const;
+                  sortmode mode = ASCENDING) const;
 
   Sparse<T> diag (octave_idx_type k = 0) const;
 
diff --git a/liboctave/lo-specfun.cc b/liboctave/lo-specfun.cc
--- a/liboctave/lo-specfun.cc
+++ b/liboctave/lo-specfun.cc
@@ -3152,7 +3152,7 @@
 // This algorithm is due to P. J. Acklam.
 // See http://home.online.no/~pjacklam/notes/invnorm/
 // The rational approximation has relative accuracy 1.15e-9 in the whole region.
-// For doubles, it is refined by a single step of Higham's 3rd order method.
+// For doubles, it is refined by a single step of Halley's 3rd order method.
 // For single precision, the accuracy is already OK, so we skip it to get
 // faster evaluation.
 
@@ -3175,7 +3175,7 @@
     {  7.784695709041462e-03,  3.224671290700398e-01,
        2.445134137142996e+00,  3.754408661907416e+00 };
 
-  static const double spi2 =  8.862269254527579e-01; // sqrt(pi)/2.
+  static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
   static const double pbreak = 0.95150;
   double ax = fabs (x), y;
 
@@ -3204,7 +3204,7 @@
   if (refine)
     {
       // One iteration of Halley's method gives full precision.
-      double u = (erf(y) - x) * spi2 * exp (y*y);
+      double u = (erf (y) - x) * spi2 * exp (y*y);
       y -= u / (1 + y*u);
     }
 
@@ -3221,6 +3221,10 @@
   return do_erfinv (x, false);
 }
 
+// The algorthim for erfcinv is an adaptation of the erfinv algorithm above
+// from P. J. Acklam.  It has been modified to run over the different input
+// domain of erfcinv.  See the notes for erfinv for an explanation.
+
 static double do_erfcinv (double x, bool refine)
 {
   // Coefficients of rational approximation.
@@ -3241,12 +3245,12 @@
        2.445134137142996e+00,  3.754408661907416e+00 };
 
   static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
-  static const double pi = 3.14159265358979323846;
-  static const double pbreak = 0.95150;
+  static const double pbreak_lo = 0.04850;  // 1-pbreak
+  static const double pbreak_hi = 1.95150;  // 1+pbreak
   double y;
 
   // Select case.
-  if (x <= 1+pbreak && x >= 1-pbreak)
+  if (x >= pbreak_lo && x <= pbreak_hi)
     {
       // Middle region.
       const double q = 0.5*(1-x), r = q*q;
@@ -3254,15 +3258,15 @@
       const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
       y = yn / yd;
     }
-  else if (x < 2.0 && x > 0.0)
+  else if (x > 0.0 && x < 2.0)
     {
       // Tail region.
       const double q = x < 1 ? sqrt (-2*log (0.5*x)) : sqrt (-2*log (0.5*(2-x)));
       const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
       const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
       y = yn / yd;
-      if (x < 1-pbreak)
-        y *= -1;
+      if (x < pbreak_lo)
+        y = -y;
     }
   else if (x == 0.0)
     return octave_Inf;
diff --git a/scripts/io/strread.m b/scripts/io/strread.m
--- a/scripts/io/strread.m
+++ b/scripts/io/strread.m
@@ -154,7 +154,21 @@
 ##
 ## @end table
 ##
-## @seealso{textscan, textread, load, dlmread, fscanf}
+## When the number of words in @var{str} doesn't match an exact multiple
+## of the number of format conversion specifiers, strread's behavior
+## depends on the last character of @var{str}:
+##
+## @table @asis
+## @item last character = "\n"
+## Data columns are padded with empty fields or Nan so that all columns
+## have equal length 
+##
+## @item last character is not "\n"
+## Data columns are not padded; strread returns columns of unequal length
+##
+## @end table
+##
+# @seealso{textscan, textread, load, dlmread, fscanf}
 ## @end deftypefn
 
 function varargout = strread (str, format = "%f", varargin)
@@ -445,7 +459,8 @@
       ## Alternative below goes by simply parsing a first grab of words
       ## and counting words until the fmt_words array is exhausted:
       iwrd = 1; iwrdp = 0; iwrdl = length (words{iwrd});
-      for ii = 1:numel (fmt_words)
+      ii = 1;
+      while ii <= numel (fmt_words)
 
         nxt_wrd = 0;
 
@@ -507,12 +522,16 @@
 
         if (nxt_wrd)
           ++iwrd; iwrdp = 0;
-          if (ii < numel (fmt_words))
+          if (iwrd > numel (words))
+            ## Apparently EOF; assume incomplete row already at L.1 of data
+            ii = numel (fmt_words);
+          elseif (ii < numel (fmt_words))
             iwrdl = length (words{iwrd});
           endif
         endif
+        ++ii;
 
-      endfor
+      endwhile
       ## Done
       words_period = max (iwrd - 1, 1);
       num_lines = ceil (num_words / words_period);
@@ -896,6 +915,16 @@
 %! assert (c, [0.94; 0.87], 0.01)
 
 %!test
+%! [a, b] = strread (['Empty 1' char(10)], 'Empty%s %f');
+%! assert (a{1}, '1');
+%! assert (b, NaN);
+
+%!test
+%! [a, b] = strread (['Empty' char(10)], 'Empty%f %f');
+%! assert (a, NaN);
+%! assert (b, NaN);
+
+%!test
 %! # Bug #35999
 %! [a, b, c] = strread ("", "%f");
 %! assert (isempty (a));
diff --git a/scripts/io/textscan.m b/scripts/io/textscan.m
--- a/scripts/io/textscan.m
+++ b/scripts/io/textscan.m
@@ -25,13 +25,13 @@
 ## @deftypefnx {Function File} {[@var{C}, @var{position}] =} textscan (@var{fid}, @dots{})
 ## Read data from a text file or string.
 ##
-## The file associated with @var{fid} is read and parsed according to
-## @var{format}.  The function behaves like @code{strread} except it works by
-## parsing a file instead of a string.  See the documentation of
-## @code{strread} for details.
+## The string @var{str} or file associated with @var{fid} is read from and
+## parsed according to @var{format}. The function behaves like @code{strread}
+## except it can also read from file instead of a string. See the documentation
+## of @code{strread} for details.
 ##
-## In addition to the options supported by
-## @code{strread}, this function supports a few more:
+## In addition to the options supported by @code{strread}, this function
+## supports a few more:
 ##
 ## @itemize
 ## @item "collectoutput":
@@ -50,16 +50,19 @@
 ## @item "returnonerror":
 ## If set to numerical 1 or true (default), return normally when read errors
 ## have been encountered.  If set to 0 or false, return an error and no data.
+## As the string or file is read by columns rather than by rows, and because
+## textscan is fairly forgiving as regards read errors, setting this option
+## may have little or no actual effect.
 ## @end itemize
 ##
 ## When reading from a character string, optional input argument @var{n}
 ## specifies the number of times @var{format} should be used (i.e., to limit
 ## the amount of data read).
-## When reading fro file, @var{n} specifies the number of data lines to read;
+## When reading from file, @var{n} specifies the number of data lines to read;
 ## in this sense it differs slightly from the format repeat count in strread.
 ##
-## The output @var{C} is a cell array whose length is given by the number
-## of format specifiers.
+## The output @var{C} is a cell array whose second dimension is determined
+## by the number of format specifiers.
 ##
 ## The second output, @var{position}, provides the position, in characters,
 ## from the beginning of the file.
@@ -80,14 +83,18 @@
     format = "%f";
   endif
 
+  if (! ischar (format))
+    error ("textscan: FORMAT must be a string");
+  endif
+
+  ## Determine the number of data fields & initialize output array
+  num_fields = numel (strfind (format, "%")) - numel (strfind (format, "%*"));
+  C = cell (1, num_fields);
+
   if (! (isa (fid, "double") && fid > 0) && ! ischar (fid))
     error ("textscan: first argument must be a file id or character string");
   endif
 
-  if (! ischar (format))
-    error ("textscan: FORMAT must be a string");
-  endif
-
   args = varargin;
   if (nargin > 2 && isnumeric (args{1}))
     nlines = args{1};
@@ -96,7 +103,6 @@
   endif
   if (nlines < 1)
     printf ("textscan: N = 0, no data read\n");
-    C = [];
     return
   endif
 
@@ -174,7 +180,6 @@
   ## Check for empty result
   if (isempty (str))
     warning ("textscan: no data read");
-    C = [];
     return;
   endif
 
@@ -249,10 +254,18 @@
   ## Determine the number of data fields
   num_fields = numel (strfind (format, "%")) - numel (strfind (format, "%*"));
 
-  ## Strip trailing EOL to avoid returning stray missing values (f. strread)
-  if (strcmp (str(end-length (eol_char) + 1 : end), eol_char));
-    str(end-length (eol_char) + 1 : end) = "";
-  endif
+  ## Strip trailing EOL to avoid returning stray missing values (f. strread).
+  ## However, in case of CollectOutput request, presence of EOL is required
+  eol_at_end = strcmp (str(end-length (eol_char) + 1 : end), eol_char);
+  if (collop)
+    if (! eol_at_end)
+      str(end+1 : end+length (eol_char)) = eol_char;
+    endif
+  elseif (eol_at_end)
+     str(end-length (eol_char) + 1 : end) = "";
+    ## A corner case: str may now be empty....
+    if (isempty (str)); return; endif
+   endif
 
   ## Call strread to make it do the real work
   C = cell (1, num_fields);
@@ -316,14 +329,14 @@
 %! assert (b(1,:)', c{1}, 1e-5);
 %! assert (b(2,:)', c{2}, 1e-5);
 
-#%!test
-#%! str = "13, 72, NA, str1, 25\r\n// Middle line\r\n36, na, 05, str3, 6";
-#%! a = textscan (str, "%d %n %f %s %n", "delimiter", ",","treatAsEmpty", {"NA", "na"},"commentStyle", "//");
-#%! assert (a{1}, int32([13; 36]));
-#%! assert (a{2}, [72; NaN]);
-#%! assert (a{3}, [NaN; 5]);
-#%! assert (a{4}, {"str1"; "str3"});
-#%! assert (a{5}, [25; 6]);
+%!test
+%! str = "13, 72, NA, str1, 25\r\n// Middle line\r\n36, na, 05, str3, 6";
+%! a = textscan (str, "%d %n %f %s %n", "delimiter", ",","treatAsEmpty", {"NA", "na"},"commentStyle", "//");
+%! assert (a{1}, int32([13; 36]));
+%! assert (a{2}, [72; NaN]);
+%! assert (a{3}, [NaN; 5]);
+%! assert (a{4}, {"str1"; "str3"});
+%! assert (a{5}, [25; 6]);
 
 %!test
 %! str = "Km:10 = hhhBjjj miles16hour\r\n";
@@ -362,6 +375,21 @@
 %! assert (size(c{3}), [10, 2]);
 %! assert (size(c{2}), [10, 2]);
 
+%!test
+%% CollectOutput test with uneven column length files
+%! b = [10:10:100];
+%! b = [b; 8*b/5; 8*b*1000/5];
+%! str = sprintf ("%g miles/hr = %g (%g) kilometers (meters)/hr\n", b);
+%! str = [str "110 miles/hr"];
+%! fmt = "%f miles%s %s %f (%f) kilometers %*s";
+%! c = textscan (str, fmt, "collectoutput", 1);
+%! assert (size(c{1}), [11, 1]);
+%! assert (size(c{3}), [11, 2]);
+%! assert (size(c{2}), [11, 2]);
+%! assert (c{3}(end), NaN);
+%! assert (c{2}{11, 1}, "/hr");
+%! assert (isempty (c{2}{11, 2}), true);
+
 %% Test input validation
 %!error textscan ()
 %!error textscan (single (4))
@@ -370,3 +398,7 @@
 %!error <cannot provide position information> [C, pos] = textscan ("Hello World")
 %!error <character value required> textscan ("Hello World", '%s', 'EndOfLine', 3)
 
+%! Test incomplete first data line
+%! R = textscan (['Empty1' char(10)], 'Empty%d %f');
+%! assert (R{1}, int32(1));
+%! assert (isempty(R{2}), true);
diff --git a/scripts/optimization/lsqnonneg.m b/scripts/optimization/lsqnonneg.m
--- a/scripts/optimization/lsqnonneg.m
+++ b/scripts/optimization/lsqnonneg.m
@@ -21,6 +21,7 @@
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d})
 ## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0})
+## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0}, @var{options})
 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{})
 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{})
 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{})
@@ -29,6 +30,9 @@
 ## Minimize @code{norm (@var{c}*@var{x} - d)} subject to
 ## @code{@var{x} >= 0}.  @var{c} and @var{d} must be real.  @var{x0} is an
 ## optional initial guess for @var{x}.
+## Currently, @code{lsqnonneg}
+## recognizes these options: @code{"MaxIter"}, @code{"TolX"}.
+## For a description of these options, see @ref{doc-optimset,,optimset}.
 ##
 ## Outputs:
 ##
@@ -147,7 +151,8 @@
     ## compute the gradient.
     w = c'*(d - c*x);
     w(p) = [];
-    if (! any (w > 0))
+    tolx = optimget (options, "TolX", 10*eps*norm (c, 1)*length (c));
+    if (! any (w > tolx))
       if (useqr)
         ## verify the solution achieved using qr updating.
         ## in the best case, this should only take a single step.
diff --git a/scripts/plot/plotyy.m b/scripts/plot/plotyy.m
--- a/scripts/plot/plotyy.m
+++ b/scripts/plot/plotyy.m
@@ -91,6 +91,7 @@
       ax = ax(1:2);
     elseif (length (ax) == 1)
       ax(2) = axes ();
+      set (ax(2), "nextplot", get (ax(1), "nextplot"))
     elseif (isempty (ax))
       ax(1) = axes ();
       ax(2) = axes ();
@@ -156,6 +157,7 @@
     axes (ax(2));
   else
     ax(2) = axes ();
+    set (ax(2), "nextplot", get (ax(1), "nextplot"))
   endif
   newplot ();
 
@@ -172,7 +174,6 @@
   set (ax(2), "yaxislocation", "right");
   set (ax(2), "ycolor", getcolor (h2(1)));
 
-
   if (strcmp (get(ax(1), "activepositionproperty"), "position"))
     set (ax(2), "position", get (ax(1), "position"));
   else
@@ -208,6 +209,8 @@
   addlistener (ax(2), "plotboxaspectratio", {@update_position, ax(1)});
   addlistener (ax(1), "plotboxaspectratiomode", {@update_position, ax(2)});
   addlistener (ax(2), "plotboxaspectratiomode", {@update_position, ax(1)});
+  addlistener (ax(1), "nextplot", {@update_nextplot, ax(2)});
+  addlistener (ax(2), "nextplot", {@update_nextplot, ax(1)});
 
   ## Store the axes handles for the sister axes.
   if (ishandle (ax(1)) && ! isprop (ax(1), "__plotyy_axes__"))
@@ -261,10 +264,24 @@
 %! clf;
 %! x = linspace (-1, 1, 201);
 %! hax = plotyy (x, sin (pi*x), x, cos (pi*x));
-%! ylabel ('Blue on the Left');
+%! ylabel (hax(1), 'Blue on the Left');
 %! ylabel (hax(2), 'Green on the Right');
 %! xlabel ('xlabel');
 
+%!demo
+%! clf
+%! hold on
+%! t = (0:0.1:9);
+%! x = sin (t);
+%! y = 5 * cos (t);
+%! [hax, h1, h2] = plotyy (t, x, t, y);
+%! [~, h3, h4] = plotyy (t+1, x, t+1, y);
+%! set ([h3, h4], "linestyle", "--")
+%! xlabel (hax(1), 'xlabel')
+%! title (hax(2), 'title')
+%! ylabel (hax(1), 'Left axis is Blue')
+%! ylabel (hax(2), 'Right axis is Green')
+
 function deleteplotyy (h, d, ax2, t2)
   if (ishandle (ax2) && strcmp (get (ax2, "type"), "axes")
       && (isempty (gcbf()) || strcmp (get (gcbf(), "beingdeleted"),"off"))
@@ -274,6 +291,19 @@
   endif
 endfunction
 
+function update_nextplot (h, d, ax2)
+  persistent recursion = false;
+  prop = "nextplot";
+  if (! recursion)
+    unwind_protect
+      recursion = true;
+      set (ax2, prop, get (h, prop));
+    unwind_protect_cleanup
+      recursion = false;
+    end_unwind_protect
+  endif
+endfunction
+
 function update_position (h, d, ax2)
   persistent recursion = false;
 
diff --git a/scripts/plot/private/__scatter__.m b/scripts/plot/private/__scatter__.m
--- a/scripts/plot/private/__scatter__.m
+++ b/scripts/plot/private/__scatter__.m
@@ -268,7 +268,8 @@
     ## Does gnuplot only support triangles with different vertex colors ?
     ## TODO - Verify gnuplot can only support one color. If RGB triplets
     ##        can be assigned to each vertex, then fix __go_draw_axe__.m
-    gnuplot_hack = numel (x) > 1 && strcmp (toolkit, "gnuplot");
+    gnuplot_hack = (numel (x) > 1 && size(c, 2) == 3
+                    && strcmp (toolkit, "gnuplot"));
     if (ischar (c) || ! isflat || gnuplot_hack)
       if (filled)
         h = __go_patch__ (hg, "xdata", x, "ydata", y, "zdata", z,
diff --git a/scripts/plot/scatter.m b/scripts/plot/scatter.m
--- a/scripts/plot/scatter.m
+++ b/scripts/plot/scatter.m
@@ -92,6 +92,14 @@
 %! clf;
 %! x = randn (100, 1);
 %! y = randn (100, 1);
+%! c = x .* y;
+%! scatter (x, y, 20, c, "filled");
+%! title ('Scatter with colors');
+
+%!demo
+%! clf;
+%! x = randn (100, 1);
+%! y = randn (100, 1);
 %! scatter (x, y, [], sqrt (x.^2 + y.^2));
 %! title ('Scatter plot with bubble color determined by distance from origin');
 
diff --git a/src/graphics.cc b/src/graphics.cc
--- a/src/graphics.cc
+++ b/src/graphics.cc
@@ -4284,7 +4284,7 @@
   xticklabelmode = "auto";
   yticklabelmode = "auto";
   zticklabelmode = "auto";
-  color = "none";
+  color = color_values ("white");
   xcolor = color_values ("black");
   ycolor = color_values ("black");
   zcolor = color_values ("black");