# HG changeset patch
# User dbateman
# Date 1135801010 0
# Node ID d37b96139376513997b3c214b76b3d08138c69c9
# Parent  6dbb3f6d005488ef8f960b1124335f13a2e0dc4f
[project @ 2005-12-28 20:16:50 by dbateman]

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,8 @@
+2005-12-28  David Bateman  <dbateman@free.fr>
+
+	* Sparse-op-defs.h (SPARSE_SPARSE_MUL): Improved algorithm that is faster in
+	all cases, and significantly so for low density or small oder problems.
+
 2005-11-30  John W. Eaton  <jwe@octave.org>
 
 	* LSODE.cc (LSODE::do_integrate (double)): Resize iwork and rwork
diff --git a/liboctave/Sparse-op-defs.h b/liboctave/Sparse-op-defs.h
--- a/liboctave/Sparse-op-defs.h
+++ b/liboctave/Sparse-op-defs.h
@@ -1547,24 +1547,26 @@
     } \
   else \
     { \
-      OCTAVE_LOCAL_BUFFER (EL_TYPE, Xcol, nr); \
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \
+      for (octave_idx_type i = 0; i < nr; i++) \
+	w[i] = 0; \
       \
       octave_idx_type nel = 0; \
       \
       for (octave_idx_type i = 0; i < a_nc; i++) \
         { \
-          OCTAVE_QUIT; \
-          for (octave_idx_type k = 0; k < nr; k++) \
-	    Xcol[k]= 0.; \
           for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
             { \
               octave_idx_type  col = a.ridx(j); \
               for (octave_idx_type k = m.cidx(col) ; k < m.cidx(col+1); k++) \
-		if (Xcol[m.ridx(k)] == 0.) \
-                  { \
-		    Xcol[m.ridx(k)] = 1.; \
-		    nel++; \
-		  } \
+		{ \
+		  OCTAVE_QUIT; \
+		  if (w[m.ridx(k)] < i + 1) \
+                    { \
+		      w[m.ridx(k)] = i + 1; \
+		      nel++; \
+		    } \
+		} \
 	    } \
 	} \
       \
@@ -1572,34 +1574,41 @@
 	return RET_TYPE (nr, a_nc); \
       else \
 	{  \
+          for (octave_idx_type i = 0; i < nr; i++) \
+	    w[i] = 0; \
+	  \
+          OCTAVE_LOCAL_BUFFER (EL_TYPE, Xcol, nr); \
+          \
 	  RET_TYPE retval (nr, a_nc, nel); \
 	  \
 	  octave_idx_type ii = 0; \
 	  \
-	  retval.cidx(0) = 0; \
+	  retval.xcidx(0) = 0; \
 	  for (octave_idx_type i = 0; i < a_nc ; i++) \
 	    { \
-              OCTAVE_QUIT; \
-	      for (octave_idx_type k = 0; k < nr; k++) \
-		Xcol[k]= 0.; \
 	      for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
 		{ \
                   octave_idx_type col = a.ridx(j); \
                   EL_TYPE tmpval = a.data(j); \
 		  for (octave_idx_type k = m.cidx(col) ; k < m.cidx(col+1); k++) \
-		    Xcol[m.ridx(k)] += tmpval * m.data(k); \
-		} \
-	      for (octave_idx_type k = 0; k < nr; k++) \
-		{ \
-		  if (Xcol[k] !=0. ) \
 		    { \
-		      retval.ridx (ii) = k; \
-		      retval.data (ii++) = Xcol[k]; \
+                      OCTAVE_QUIT; \
+		      octave_idx_type row = m.ridx(k); \
+		      if (w[row] < i + 1) \
+		        { \
+		          w[row] = i + 1; \
+		          retval.xridx(ii++) = row; \
+		          Xcol[row] = tmpval * m.data(k); \
+		        } \
+		      else \
+		        Xcol[row] += tmpval * m.data(k); \
 		    } \
 		} \
-	      retval.cidx(i+1) = ii; \
+	      retval.xcidx(i+1) = ii; \
+	      for (octave_idx_type k = retval.xcidx(i); k < retval.xcidx(i+1); k++) \
+		retval.xdata(k) = Xcol[retval.xridx(k)]; \
 	    } \
-	  retval.maybe_compress (); \
+	  retval.maybe_compress (true); \
 	  return retval; \
 	} \
     }