function P = prodmat(A,B)

if ~isequal(size(A,2),size(B,1))
    error('Aaa!!!');
end

MAX_ROWS = 1000;

P = sparse(size(A,1),size(B,2));

blocksA = ceil(size(A,1)/MAX_ROWS);
blocksB = ceil(size(B,2)/MAX_ROWS);
for i=1:blocksA
    for j=1:blocksB
        if i==blocksA
            linesA = (i-1)*MAX_ROWS+1:size(A,1);
        else
            linesA = (i-1)*MAX_ROWS+1:i*MAX_ROWS;
        end
        
        if j==blocksB
            linesB = (j-1)*MAX_ROWS+1:size(B,2);
        else
            linesB = (j-1)*MAX_ROWS+1:j*MAX_ROWS;
        end
        
        P(linesA,linesB) = A(linesA,:)*B(:,linesB);
    end
end

end

