# sparse matrix routines

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