Speed Problem of Many Body Operator

[Umut-Can-Physics]

Many body operator in the module manybodyoperator(basis::ManyBodyBasis, op) works slowly. My method is faster than module's method such as
for i in 1:N for j in 1:N mb_op += op.data[i,j] * transition(ManyBodyBasis, i, j) end end
As you can see even having nested loop (dimension is N^2), my method still faster than original one. What is the fastest way to construct many body operator? Note that my problem dimension generally is from N=50 up to N=100.

[Krastanov]

Could you give a bit more details, e.g. a minimum working example for each approach (your hand-rolled one and the one where you use manybodyoperator). If we can just copy paste your code and have it run, it would be much easier for us to answer questions. You can use @time to report how long the runs take as well.

[Umut-Can-Physics]

Here is an example

using BenchmarkTools
using QuantumOptics
# FIRST METHOD
N = 64
b = NLevelBasis(N)
op = randoperator(b)
M = 2
s = bosonstates(b, [M])
mb = ManyBodyBasis(b, s)
mb_op = manybodyoperator(mb, op) # MB OPERATOR 1
# SECOND METHOD
sparse_mb_op = SparseOperator(mb) # MB OPERATOR 2
for i in 1:N 
    for j in 1:N 
        sparse_mb_op += op.data[i,j] * transition(mb, i, j)  
    end 
end 

Note that two operator MB OPERATOR 1 and MB OPERATOR 2 are the same but first method is so slow, WHY? Also, I need to speed up the second method.

[Umut-Can-Physics]

Update: Now, I use the following methode that is faster than the other methods.

using BenchmarkTools
using QuantumOptics
using LinearAlgebra
using Base.Threads

# Define the system and operator
N = 64
b = NLevelBasis(N)
op = randoperator(b)
M = 2
s = bosonstates(b, [M])
mb = ManyBodyBasis(b, s)

function MB_Op(mb, op)
    # Prepare a SparseOperator for each thread
    n_threads = Threads.nthreads()
    thread_results = [SparseOperator(mb) for _ in 1:n_threads]

    # Parallel loop with independent storage for each thread
    @threads for i in 1:N
        local_op = thread_results[threadid()]
        for j in 1:N
            if op.data[i, j] != 0  # Skip zero entries for efficiency
                local_op += op.data[i, j] * transition(mb, i, j)
            end
        end
        thread_results[threadid()] = local_op  # Ensure the result is stored in thread_results
    end

    # Combine the results from all threads
    sparse_mb_op = sum(thread_results)

    return sparse_mb_op
end