[PATCH] Eliminate the workspace in sparse transpose.

[Jason Riedy]

# HG changeset patch
# User Jason Riedy <address@hidden>
# Date 1236708494 14400
# Node ID fa3bc420d7ddb9fc1df3e095751567cbbc5c15bd
# Parent  42e24f4ebc8c70db07d0fcfd08d70b7a40e5519c
[PATCH] Eliminate the workspace in sparse transpose.

>From 226faf8f1680a9eb5ee8a736dfc9244248ba2984 Mon Sep 17 00:00:00 2001
Date: Mon, 9 Mar 2009 23:35:15 -0400
The output's cidx (column start offset array) can serve as the
workspace, so the routines operate in the space of their output.

Signed-off-by: Jason Riedy <address@hidden>
---
 liboctave/CSparse.cc |   22 +++++++++++-----------
 liboctave/ChangeLog  |    7 +++++++
 liboctave/Sparse.cc  |   22 +++++++++++-----------
 3 files changed, 29 insertions(+), 22 deletions(-)

diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -640,28 +640,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-10  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-08  Jaroslav Hajek  <address@hidden>
 
        * idx-vector.h (idx_vector::bloop): loop --> bloop.
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;
 }