Taking advantage of rowwise-constant sparse matrix

[Sukera]

Apologies, I was cooking my dinner. That’s why I didn’t answer right away in the other thread.

@lrnv

function my_prod(r,b,M,x)
    @inbounds for i in size(M,1)
        r[i] = b[i]*sum(x[M[i,:]])
    end
    return r
end

This function doesn’t do what you think it does.

• size(M,1) is a single integer - the loop is only running for a single iteration because numbers are iterable in julia and thus the result must be wrong. You probably meant to write axes(M, 1) or 1:size(M,1). The former gives an iterable over the first axis, the latter constructs the same manually. Your check still passes because only the last entry is overwritten and r still holds the correct result from P*x.

julia> r_ = zeros(size(r)); 
                            
julia> my_prod(r_,b,M,x) ≈ r
false                       
                            
julia> findall(!=(0), r_)   
1-element Vector{Int64}:    
 1000                       

• Since M[i,:] is a vector, x[M[i,:]] will slice (and thus allocate). This doesn’t hurt too bad for one iteration, but it will kill performance for all appropriate iterations:

julia> function my_prod(r,b,M,x)         
           for i in 1:size(M,1) # fixed iterations         
               r[i] = b[i]*sum(x[M[i,:]])
           end                           
           return r                      
       end                               

my_prod (generic function with 1 method) 
julia> @btime P*x;                        
  29.900 μs (1 allocation: 7.94 KiB)      
                                          
julia> @btime my_prod(r,b,M,$x);          
  39.710 ms (7038 allocations: 952.67 KiB)

julia> r_ = zeros(size(r)); 
                            
julia> my_prod(r_,b,M,x) ≈ r
true