### Usage:
###
### X = sp_multiply(A,B,X_template)
###
### This function multiplies two sparce matrices so
### that only those non-zero elements in X_template
### are returned in the sparse matrix X.
###
function X = sp_multiply(A,B,X_template)

    if (nargin ~= 3)
        usage('X = sp_multiply(A,B,X_template)\n');
    endif

    if (rows(X_template) ~= rows(A) ||
        columns(X_template) ~= columns(B))
        error('Template matrix has wrong dimensions');
    endif

    [i,j]=find(X_template);
    tmp=zeros(length(i),1);

    for k=1:length(i)
        ## The "full" command is invoked because some
        ## sort of overflow error occurs in Octave
        ## 2.1.57 if A and B are both sparse
        tmp(k) = full(A(i(k),:)) * full(B(:,j(k)));       
    endfor

    [r,c]=size(X_template);
    X=sparse(i,j,tmp,r,c);

endfunction