I thought I would share some code for computing sparse matrix products in Maple.  For floating point matrices this is done quickly, but for algebraic datatypes there is a performance problem with the builtin routines, as noted in http://www.mapleprimes.com/questions/205739-How-Do-I-Improve-The-Performance-Of

Download spm.txt

The code is fairly straightforward in that it uses op(1,A) to extract the dimensions and op(2,A) to extract the non-zero elements of a Matrix or Vector, and then loops over those elements.  I included a sparse map function for cases where you want to map a function (like expand) over non-zero elements only.

# sparse matrix vector product
spmv := proc(A::Matrix,V::Vector)
local m,n,Ae,Ve,Vi,R,e;
  n, m := op(1,A);
  if op(1,V) <> m then error "incompatible dimensions"; end if;
  Ae := op(2,A);
  Ve := op(2,V);
  Vi := map2(op,1,Ve);
  R := Vector(n, storage=sparse);
  for e in Ae do
    n, m := op(1,e);
    if member(m, Vi) then R[n] := R[n] + A[n,m]*V[m]; end if;
  end do;
  return R;
end proc:

# sparse matrix product
spmm := proc(A::Matrix, B::Matrix)
local m,n,Ae,Be,Bi,R,l,e,i;
  n, m := op(1,A);
  i, l := op(1,B);
  if i <> m then error "incompatible dimensions"; end if;
  Ae := op(2,A);
  Be := op(2,B);
  R := Matrix(n,l,storage=sparse);
  for i from 1 to l do
    Bi, Be := selectremove(type, Be, (anything,i)=anything);
    Bi := map2(op,[1,1],Bi);
    for e in Ae do
      n, m := op(1,e);
      if member(m, Bi) then R[n,i] := R[n,i] + A[n,m]*B[m,i]; end if;
    end do;
  end do;
  return R;
end proc:

# sparse map
smap := proc(f, A::{Matrix,Vector})
local B, Ae, e;
  if A::Vector then
    B := Vector(op(1,A),storage=sparse):
  else
    B := Matrix(op(1,A),storage=sparse):
  end if;
  Ae := op(2,A);
  for e in Ae do
    B[op(1,e)] := f(op(2,e),args[3..nargs]);
  end do;
  return B;
end proc:

As for how it performs, here is a demo inspired by the original post.

n := 674;
k := 6;

A := Matrix(n,n,storage=sparse):
for i to n do
  for j to k do
    A[i,irem(rand(),n)+1] := randpoly(x):
  end do:
end do:

V := Vector(n):
for i to k do
  V[irem(rand(),n)+1] := randpoly(x):
end do:

C := CodeTools:-Usage( spmv(A,V) ):  # 7ms, 25x faster
CodeTools:-Usage( A.V ):  # 174 ms

B := Matrix(n,n,storage=sparse):
for i to n do
  for j to k do
    B[i,irem(rand(),n)+1] := randpoly(x):
  end do:
end do:

C := CodeTools:-Usage( spmm(A,B) ):  # 2.74 sec, 50x faster
CodeTools:-Usage( A.B ):  # 2.44 min

# expand and collect like terms
C := CodeTools:-Usage( smap(expand, C) ):