Sparse Transpose

[David Bateman]

Find attached a small patch that vastly accelerates the sparse transpose
function. Running

n=1e4;for
i=-5:-2,a=sprandn(n,n,10^i);t=cputime();b=a';t0=cputime()-t;printf("%6g
%6g %f\n", n, 10^i, t0); fflush(stdout); endfor

I get

   n          d        time
 10000  1e-05 1.092834
 10000 0.0001 2.678592
 10000  0.001 10.210447
  (Ctrl-c since my patience ran out)

before and after the patch I get

 10000  1e-05 0.000999
 10000 0.0001 0.001000
 10000  0.001 0.006000
 10000   0.01 0.091985

for a factor of between 100 to 1000 speed-up.

Regards
David

2006-03-07  David Bateman  <address@hidden>

    * Sparse.cc (Sparse<T>::transpose (void) const): Accelerate algorithm.
    * CSparse.cc (SparseComplexMatrix::transpose (void) const): ditto.


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

Index: liboctave/Sparse.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/Sparse.cc,v
retrieving revision 1.8
diff -c -r1.8 Sparse.cc
*** liboctave/Sparse.cc 20 Feb 2006 22:05:31 -0000      1.8
--- liboctave/Sparse.cc 7 Mar 2006 22:54:23 -0000
***************
*** 1034,1054 ****
  
    octave_idx_type nr = rows ();
    octave_idx_type nc = cols ();
!   octave_idx_type nz = nzmax ();
    Sparse<T> retval (nc, nr, nz);
  
!   retval.cidx(0) = 0;
!   for (octave_idx_type i = 0, iidx = 0; i < nr; i++)
!     {
!       for (octave_idx_type j = 0; j < nc; j++)
!       for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
!         if (ridx(k) == i)
!           {
!             retval.data(iidx) = data(k);
!             retval.ridx(iidx++) = j;
!           }
!       retval.cidx(i+1) = iidx;
!     }
  
    return retval;
  }
--- 1034,1064 ----
  
    octave_idx_type nr = rows ();
    octave_idx_type nc = cols ();
!   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)]++;
!   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 j = 0; j < nc; j++)
!     for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
!       {
!       octave_idx_type q = w [ridx(k)]++;
!       retval.xridx (q) = j;
!       retval.xdata (q) = data (k);
!       }
  
    return retval;
  }
Index: liboctave/CSparse.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/CSparse.cc,v
retrieving revision 1.18
diff -c -r1.18 CSparse.cc
*** liboctave/CSparse.cc        20 Feb 2006 22:05:30 -0000      1.18
--- liboctave/CSparse.cc        7 Mar 2006 22:55:20 -0000
***************
*** 549,566 ****
    octave_idx_type nz = nzmax ();
    SparseComplexMatrix retval (nc, nr, nz);
  
!   retval.cidx(0) = 0;
!   for (octave_idx_type i = 0, iidx = 0; i < nr; i++)
!     {
!       for (octave_idx_type j = 0; j < nc; j++)
!       for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
!         if (ridx(k) == i)
!           {
!             retval.data(iidx) = conj (data(k));
!             retval.ridx(iidx++) = j;
!           }
!       retval.cidx(i+1) = iidx;
!     }
  
    return retval;
  }
--- 549,576 ----
    octave_idx_type nz = nzmax ();
    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)]++;
!   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 j = 0; j < nc; j++)
!     for (octave_idx_type k = cidx(j); k < cidx(j+1); k++)
!       {
!       octave_idx_type q = w [ridx(k)]++;
!       retval.xridx (q) = j;
!       retval.xdata (q) = conj (data (k));
!       }
  
    return retval;
  }