diff -r 530803d4f65f liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc	Wed Jun 10 15:21:24 2015 -0700
+++ b/liboctave/array/dSparse.cc	Tue Jun 23 20:31:39 2015 +0200
@@ -7560,7 +7560,209 @@
 SparseMatrix
 operator * (const SparseMatrix& m, const SparseMatrix& a)
 {
-  SPARSE_SPARSE_MUL (SparseMatrix, double, double);
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.cols ();
+
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (nr == 1 && nc == 1)
+   {
+     double s = m.elem (0,0);
+     octave_idx_type nz = a.nnz ();
+     SparseMatrix r (a_nr, a_nc, nz);
+
+     for (octave_idx_type i = 0; i < nz; i++)
+       {
+         octave_quit ();
+         r.data (i) = s * a.data (i);
+         r.ridx (i) = a.ridx (i);
+       }
+     for (octave_idx_type i = 0; i < a_nc + 1; i++)
+       {
+         octave_quit ();
+         r.cidx (i) = a.cidx (i);
+       }
+
+     r.maybe_compress (true);
+     return r;
+   }
+  else if (a_nr == 1 && a_nc == 1)
+   {
+     double s = a.elem (0,0);
+     octave_idx_type nz = m.nnz ();
+     SparseMatrix r (nr, nc, nz);
+
+     for (octave_idx_type i = 0; i < nz; i++)
+       {
+         octave_quit ();
+         r.data (i) = m.data (i) * s;
+         r.ridx (i) = m.ridx (i);
+       }
+     for (octave_idx_type i = 0; i < nc + 1; i++)
+       {
+         octave_quit ();
+         r.cidx (i) = m.cidx (i);
+       }
+
+     r.maybe_compress (true);
+     return r;
+   }
+  else if (nc != a_nr)
+    {
+      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
+      return SparseMatrix ();
+    }
+  else
+  {
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr);
+      SparseMatrix retval (nr, a_nc, static_cast<octave_idx_type> (0));
+      for (octave_idx_type i = 0; i < nr; i++)
+          w[i] = 0;
+      retval.xcidx (0) = 0;
+
+      octave_idx_type nel = 0;
+
+      for (octave_idx_type i = 0; i < a_nc; i++)
+      {
+          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 (w[m.ridx (k)] < i + 1)
+                  {
+                      w[m.ridx (k)] = i + 1;
+                      nel++;
+                  }
+                  octave_quit ();
+              }
+          }
+          retval.xcidx (i+1) = nel;
+      }
+
+      if (nel == 0)
+      {
+         return SparseMatrix (nr, a_nc);
+      }
+      else
+      {
+          CXSPARSE_DNAME () csm, csa;
+
+          csm.m = m.rows();
+          csa.m = a.rows();
+
+          csm.n = m.cols();
+          csa.n = a.cols();
+
+          // use -1 to indicate compressed column space
+          csm.nz = -1;
+          csa.nz = -1;
+
+          csm.nzmax = m.nnz ();
+          csa.nzmax = a.nnz ();
+
+          // Cast away const on A, with full knowledge that CSparse won't touch it.
+          // Prevents the methods below making a copy of the data.
+          csm.p = const_cast<octave_idx_type *>(m.cidx ());
+          csa.p = const_cast<octave_idx_type *>(a.cidx ());
+          csm.i = const_cast<octave_idx_type *>(m.ridx ());
+          csa.i = const_cast<octave_idx_type *>(a.ridx ());
+          csm.x = const_cast<double *>(m.data ());
+          csa.x = const_cast<double *>(a.data ());
+
+          CXSPARSE_DNAME () * csret_p = CXSPARSE_DNAME (_multiply)(&csm, &csa);
+          if (csret_p == NULL)
+          {
+              (*current_liboctave_error_handler)
+                  ("SuiteSparse_interface: multiply returned null");
+          }
+
+          // sometime suitesparse returns zero entries
+          // for some of the weights. we drop those
+          CXSPARSE_DNAME (_dropzeros) (csret_p);
+
+          // this sorts the content of the row and weights matrix
+          // sparse<T> has a sort method, but simply calling that does not
+          // pass tests
+          // we need auxiliary vectors to store the original weights and
+          // a range, to sort the weights in the same way the rows have been
+          // sorted
+
+          std::vector<octave_idx_type> perm_vec;
+          std::vector<double> orig_weight_vec;
+
+          //the maximum size of the workspace is m (numbers of rows)
+          perm_vec.reserve(csret_p->m);
+          orig_weight_vec.reserve(csret_p->m);
+
+          // we operate in place, so const this for performance and clearness
+          double * const ret_data = csret_p->x;
+
+          for (octave_idx_type i = 0; i < csret_p->n; i++)
+          {
+              perm_vec.clear();
+              orig_weight_vec.clear();
+
+              const octave_idx_type start_idx = (csret_p->p)[i];
+              const octave_idx_type end_idx = (csret_p->p)[i+1];
+
+              // stupid way to generate a range and save the original weights,
+              // maybe something better exists in the codebase (could not
+              // really see how to use octave_range)
+              octave_idx_type k = 0;
+              for (octave_idx_type j = start_idx; j < end_idx; j++ , k++)
+              {
+                  perm_vec.push_back(k);
+                  orig_weight_vec.push_back(ret_data[j]);
+              }
+
+              octave_sort<octave_idx_type> lsort;
+              lsort.set_compare (ASCENDING);
+
+              // shuffle the permutation indexes as well as the row indexes
+              // Two notes:
+              // 1) you cannot pass directly the weights as second argument
+              // or they will be casted as integers, while they are float.
+              // this might call for a reimplementation or a templated
+              // version
+              // template <typename T1, typename T2> sort(T1 *, T2 *, point)
+              // and then the current implementation would be
+              // < T1, octave_idx_type > sort()
+              // 2) not so nice hack to access the beginning address of the
+              // vector storage. It works due to the definition of vector,
+              // and the fact that the vector is not grown during the
+              // sorting and so it does not need to reallocate.
+              lsort.sort(csret_p->i + start_idx, &(*perm_vec.begin()), end_idx - start_idx);
+
+              for (octave_idx_type j = start_idx; j < end_idx; j++ )
+              {
+                  ret_data[j] = orig_weight_vec[perm_vec[j-start_idx]];
+              }
+          }
+
+          // this is pretty much the only way I found to avoid copying over the
+          // data. ret_rep.rep is declared public, but a comment in the
+          // source code says it should be protected. In both cases this would
+          // still compile, but it is very coupling.
+          Sparse<double> ret_rep = Sparse<double>(csret_p->m, csret_p->n);
+          free(ret_rep.rep->c);
+          ret_rep.rep->c = csret_p->p;
+          free(ret_rep.rep->r);
+          ret_rep.rep->r = csret_p->i;
+          free(ret_rep.rep->d);
+          ret_rep.rep->d = csret_p->x;
+          ret_rep.rep->nzmx = csret_p->nzmax;
+
+          CXSPARSE_DNAME (_free) (csret_p);
+
+          SparseMatrix ret = SparseMatrix(static_cast< const Sparse<double> >(ret_rep));
+
+
+          return ret;
+      }
+  }
 }
 
 Matrix
