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) ):