Sparse Transpose

[John W. Eaton]

On  7-Mar-2006, David Bateman wrote:

| 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.

Excellent.  Please check in this patch.

jwe


| 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;
|   }