[PATCH resend] Eliminate the workspace in sparse transpose.

[Jason Riedy]

I expect it was missed in the middle of the diag/perm bits.

# HG changeset patch
# User Jason Riedy <address@hidden>
# Date 1237237387 14400
# Node ID ce94a9dfdd97c0d46157fc476efe3e7eccee15aa
# Parent  193804a4f82f14a492daf2fe97c6bcffdac12475
Eliminate the workspace in sparse transpose.

The output's cidx (column start offset array) can serve as the
workspace, so the routines operate in the space of their output.

diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc
--- 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++)
-    {
-      retval.xcidx(i) = nz;
-      nz += w[i];
-      w[i] = retval.xcidx(i);
-    }
-  retval.xcidx(nr) = nz;
-  w[nr] = nz;
+  for (octave_idx_type i = 1; i <= nr; i++)
+    {
+      const octave_idx_type tmp = retval.xcidx (i);
+      retval.xcidx (i) = nz;
+      nz += tmp;
+    }
+  // 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;
 }
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,10 @@
+2009-03-16  Jason Riedy  <address@hidden>
+
+       * 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  <address@hidden>
 
        * mx-op-decl.h (NDS_BOOL_OP_DECLS, SND_BOOL_OP_DECLS, 
NDND_BOOL_OP_DECLS): Support compound binary ops.
diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc
--- 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;
 }