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