# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1271163541 -7200
# Node ID f0266ee4aabe34de4e309b1499d5aaa91cc72daa
# Parent  189274f6c7c42dc037d77278ab2ef53579ccd2da
optimize some special indexing & assignment cases

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,12 @@
+2010-04-13  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Sparse.cc (Sparse<T>::index): If S is a sparse column vector,
+	forward S(I,1) and S(I,:) to 1D indexing. Handle permutation indexing
+	in the 1D case.
+	(Sparse<T>::assign): If S is a sparse column vector,
+	forward S(I,1) = X and S(I,:) =X to 1D indexed assignment.
+	* idx-vector.cc (idx_vector::inverse_permutation): Add missing break.
+
 2010-04-13  Jaroslav Hajek  <highegg@gmail.com>
 
 	* Array-util.cc (gripe_invalid_assignment_size,
diff --git a/liboctave/MSparse.h b/liboctave/MSparse.h
--- a/liboctave/MSparse.h
+++ b/liboctave/MSparse.h
@@ -55,6 +55,9 @@
 
   MSparse (const Sparse<T>& a) : Sparse<T> (a) { }
 
+  template <class U>
+  MSparse (const Sparse<U>& a) : Sparse<T> (a) { }
+
   MSparse (const Array<T>& a, const idx_vector& r, const idx_vector& c,
            octave_idx_type nr = -1, octave_idx_type nc = -1, bool sum_terms = true)
     : Sparse<T> (a, r, c, nr, nc, sum_terms) { }
diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc
--- a/liboctave/Sparse.cc
+++ b/liboctave/Sparse.cc
@@ -171,29 +171,6 @@
 }
 
 template <class T>
-template <class U>
-Sparse<T>::Sparse (const Sparse<U>& a)
-  : dimensions (a.dimensions)
-{
-  if (a.nnz () == 0)
-    rep = new typename Sparse<T>::SparseRep (rows (), cols());
-  else
-    {
-      rep = new typename Sparse<T>::SparseRep (rows (), cols (), a.nnz ());
-      
-      octave_idx_type nz = a.nnz ();
-      octave_idx_type nc = cols ();
-      for (octave_idx_type i = 0; i < nz; i++)
-        {
-          xdata (i) = T (a.data (i));
-          xridx (i) = a.ridx (i);
-        }
-      for (octave_idx_type i = 0; i < nc + 1; i++)
-        xcidx (i) = a.cidx (i);
-    }
-}
-
-template <class T>
 Sparse<T>::Sparse (octave_idx_type nr, octave_idx_type nc, T val)
   : dimensions (dim_vector (nr, nc))
 { 
@@ -1395,6 +1372,21 @@
           mx_inline_sub (nz_new, retval.xridx (), ridx () + li, lb);
           retval.xcidx(1) = nz_new;
         }
+      else if (idx.is_permutation (nel) && idx.is_vector ())
+        {
+          if (idx.is_range () && idx.increment () == -1)
+            {
+              retval = Sparse<T> (nr, 1, nz);
+              std::reverse_copy (ridx (), ridx () + nz, retval.ridx ());
+              std::reverse_copy (data (), data () + nz, retval.data ());
+            }
+          else
+            {
+              Array<T> tmp = array_value ();
+              tmp = tmp.index (idx);
+              retval = Sparse<T> (tmp);
+            }
+        }
       else
         {
           // If indexing a sparse column vector by a vector, the result is a
@@ -1557,6 +1549,11 @@
             }
         }
     }
+  else if (nc == 1 && idx_j.is_colon_equiv (nc) && idx_i.is_vector ())
+    {
+      // It's actually vector indexing. The 1D index is specialized for that.
+      retval = index (idx_i);
+    }
   else if (idx_i.is_scalar ())
     {
       octave_idx_type ii = idx_i(0);
@@ -1652,10 +1649,8 @@
       else
         {
           // Get inverse permutation.
-          OCTAVE_LOCAL_BUFFER (octave_idx_type, iinv, nr);
-          const Array<octave_idx_type> ia = idx_i.as_array ();
-          for (octave_idx_type i = 0; i < nr; i++)
-            iinv[ia(i)] = i;
+          idx_vector idx_iinv = idx_i.inverse_permutation (nr);
+          const octave_idx_type *iinv = idx_iinv.raw ();
 
           // Scatter buffer.
           OCTAVE_LOCAL_BUFFER (T, scb, nr);
@@ -1682,14 +1677,20 @@
         }
 
     }
+  else if (idx_j.is_colon ())
+    {
+      // This requires  uncompressing columns, which is generally costly,
+      // so we rely on the efficient transpose to handle this.
+      // It may still make sense to optimize some cases here.
+      retval = transpose ();
+      retval = retval.index (idx_vector::colon, idx_i);
+      retval = retval.transpose ();
+    }
   else
     {
-      // This is the most general case, where all optimizations failed.
-      // I suppose this is a relatively rare case, so it will be done
-      // as s(i,j) = ((s(:,j).')(:,i)).'
-      // Note that if j is :, the first indexing expr. is a shallow copy.
-      retval = index (idx_vector::colon, idx_j).transpose ();
-      retval = retval.index (idx_vector::colon, idx_i).transpose ();
+      // A(I, J) is decomposed into A(:, J)(I, :).
+      retval = index (idx_vector::colon, idx_j);
+      retval = retval.index (idx_i, idx_vector::colon);
     }
 
   return retval;
@@ -1988,6 +1989,11 @@
 
             }
         }
+      else if (nc == 1 && idx_j.is_colon_equiv (nc) && idx_i.is_vector ())
+        {
+          // It's actually vector indexing. The 1D assign is specialized for that.
+          assign (idx_i, rhs);
+        }
       else if (idx_j.is_colon ())
         {
           if (idx_i.is_permutation (nr))
diff --git a/liboctave/idx-vector.cc b/liboctave/idx-vector.cc
--- a/liboctave/idx-vector.cc
+++ b/liboctave/idx-vector.cc
@@ -1178,6 +1178,7 @@
         for (octave_idx_type i = 0; i < n; i++)
           idx.xelem(ri[i]) = i;
         retval = new idx_vector_rep (idx, r->extent (0), DIRECT);
+        break;
       }
     default:
       retval = *this;
@@ -1259,6 +1260,12 @@
 {
   return rep->as_array ();
 }
+
+bool
+idx_vector::is_vector (void) const
+{
+  return idx_class () != class_vector || orig_dimensions ().is_vector ();
+}
     
 octave_idx_type 
 idx_vector::freeze (octave_idx_type z_len, const char *, bool resize_ok)
diff --git a/liboctave/idx-vector.h b/liboctave/idx-vector.h
--- a/liboctave/idx-vector.h
+++ b/liboctave/idx-vector.h
@@ -1005,6 +1005,8 @@
   // Raw pointer to index array.  This is non-const because it may be
   // necessary to mutate the index.
   const octave_idx_type *raw (void);
+
+  bool is_vector (void) const;
     
   // FIXME -- these are here for compatibility.  They should be removed
   // when no longer in use.