# HG changeset patch
# User jwe
# Date 1139425014 0
# Node ID 4b33d802ef3cb3ee1fbc72d481f09b3ec38bbdb7
# Parent  70ef31ebe1568c2648531ba3223092c4c134c8ab
[project @ 2006-02-08 18:56:54 by jwe]

diff --git a/liboctave/Array-util.cc b/liboctave/Array-util.cc
--- a/liboctave/Array-util.cc
+++ b/liboctave/Array-util.cc
@@ -451,25 +451,6 @@
   return retval;
 }
 
-Array<octave_idx_type>
-calc_permutated_idx (const Array<octave_idx_type>& old_idx, 
-		     const Array<octave_idx_type>& perm_vec, bool inv)
-{
-  octave_idx_type n_el = old_idx.length ();
-
-  Array<octave_idx_type> retval (n_el);
-
-  for (octave_idx_type i = 0; i < n_el; i++)
-    {
-      if (inv)
-	retval(perm_vec(i)) = old_idx(i);
-      else
-	retval(i) = old_idx(perm_vec(i));
-    }
-
-  return retval;
-}
-
 void
 gripe_nonconformant (const char *op, int op1_len, int op2_len)
 {
diff --git a/liboctave/Array-util.h b/liboctave/Array-util.h
--- a/liboctave/Array-util.h
+++ b/liboctave/Array-util.h
@@ -79,9 +79,6 @@
 				const dim_vector& dimensions,
 				int resize_ok);
 
-extern Array<octave_idx_type> calc_permutated_idx (const Array<octave_idx_type>& old_idx, 
-				       const Array<octave_idx_type>& perm_vec, bool inv);
-
 extern void gripe_nonconformant (const char *op, int op1_len, int op2_len);
 
 extern void gripe_nonconformant (const char *op, int op1_nr, int op1_nc,
diff --git a/liboctave/Array.cc b/liboctave/Array.cc
--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -30,6 +30,7 @@
 #include <climits>
 
 #include <iostream>
+#include <vector>
 
 #include "Array.h"
 #include "Array-flags.h"
@@ -423,12 +424,30 @@
   return retval;
 }
 
+struct
+permute_vector
+{
+  octave_idx_type pidx;
+  octave_idx_type iidx;
+};
+
+static int
+permute_vector_compare (const void *a, const void *b)
+{
+  const permute_vector *pva = static_cast<const permute_vector *> (a);
+  const permute_vector *pvb = static_cast<const permute_vector *> (b);
+
+  return pva->pidx > pvb->pidx;
+}
+
 template <class T>
 Array<T>
-Array<T>::permute (const Array<octave_idx_type>& perm_vec, bool inv) const
+Array<T>::permute (const Array<octave_idx_type>& perm_vec_arg, bool inv) const
 {
   Array<T> retval;
 
+  Array<octave_idx_type> perm_vec = perm_vec_arg;
+
   dim_vector dv = dims ();
   dim_vector dv_new;
 
@@ -443,8 +462,6 @@
   // Append singleton dimensions as needed.
   dv.resize (perm_vec_len, 1);
 
-  const Array<T> tmp = reshape (dv);
-
   // Need this array to check for identical elements in permutation array.
   Array<bool> checked (perm_vec_len, false);
 
@@ -476,24 +493,76 @@
       dv_new(i) = dv(perm_elt);
     }
 
+  int nd = dv.length ();
+
+  // XXX FIXME XXX -- it would be nice to have a sort method in the
+  // Array class that also returns the sort indices.
+
+  if (inv)
+    {
+      OCTAVE_LOCAL_BUFFER (permute_vector, pvec, nd);
+
+      for (int i = 0; i < nd; i++)
+	{
+	  pvec[i].pidx = perm_vec(i);
+	  pvec[i].iidx = i;
+	}
+
+      octave_qsort (pvec, static_cast<size_t> (nd),
+		    sizeof (permute_vector), permute_vector_compare);
+
+      for (int i = 0; i < nd; i++)
+	{
+	  perm_vec(i) = pvec[i].iidx;
+	  dv_new(i) = dv(perm_vec(i));
+	}
+    }
+
   retval.resize (dv_new);
 
-  // Index array to the original array.
-  Array<octave_idx_type> old_idx (perm_vec_len, 0);
-
-  // Number of elements in Array (should be the same for
-  // both the permuted array and original array).
-  octave_idx_type n = retval.length ();
-
-  // Permute array.
+  Array<octave_idx_type> cp (nd+1, 1);
+  for (octave_idx_type i = 1; i < nd+1; i++)
+    cp(i) = cp(i-1) * dv(i-1);
+
+  octave_idx_type incr = cp(perm_vec(0));
+
+  Array<octave_idx_type> base_delta (nd-1, 0);
+  Array<octave_idx_type> base_delta_max (nd-1);
+  Array<octave_idx_type> base_incr (nd-1);
+  for (octave_idx_type i = 0; i < nd-1; i++)
+    {
+      base_delta_max(i) = dv_new(i+1);
+      base_incr(i) = cp(perm_vec(i+1));
+    }
+
+  octave_idx_type nr_new = dv_new(0);
+  octave_idx_type nel_new = dv_new.numel ();
+  octave_idx_type n = nel_new / nr_new;
+
+  octave_idx_type k = 0;
+
   for (octave_idx_type i = 0; i < n; i++)
     {
-      // Get the idx of permuted array.
-      Array<octave_idx_type> new_idx = calc_permutated_idx (old_idx, perm_vec, inv);
-
-      retval.elem (new_idx) = tmp.elem (old_idx);
-
-      increment_index (old_idx, dv);
+      octave_idx_type iidx = 0;
+      for (octave_idx_type kk = 0; kk < nd-1; kk++)
+	iidx += base_delta(kk) * base_incr(kk);
+
+      for (octave_idx_type j = 0; j < nr_new; j++)
+	{
+	  retval(k++) = elem(iidx);
+	  iidx += incr;
+	}
+
+      base_delta(0)++;
+
+      for (octave_idx_type kk = 0; kk < nd-2; kk++)
+	{
+	  if (base_delta(kk) == base_delta_max(kk))
+	    {
+	      base_delta(kk) = 0;
+	      base_delta(kk+1)++;
+	    }
+	}
     }
 
   retval.chop_trailing_singletons ();
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,11 @@
+2006-02-08  John W. Eaton  <jwe@octave.org>
+
+	* Array-util.h (calc_permutated_idx): Delete.
+	* Array.cc (permute_vector): New data structure.
+	(permute_vector_compare): New function.
+	(Array<T>::permute): Rewrite to avoid calc_permutated_index for
+	improved performance.
+
 2006-02-04  David Bateman  <dbateman@free.fr>
 
 	* COLAMD: Remove all files, as now unused.
diff --git a/src/data.cc b/src/data.cc
--- a/src/data.cc
+++ b/src/data.cc
@@ -872,7 +872,7 @@
     {
       Array<int> vec = args(1).int_vector_value ();
 
-      // XXX FIXME XXX -- maybe we shoudl create an idx_vector object
+      // XXX FIXME XXX -- maybe we should create an idx_vector object
       // here and pass that to permute?
 
       int n = vec.length ();