benoitseron / Permanents.jl

Functions to compute matrix permanents
MIT License
8 stars 3 forks source link

code is not multi threading #6

Open benkj opened 2 years ago

benkj commented 2 years ago

I'm starting to play with multi-threading, with the long term goal of using the GPU for computing the permanent. So far I've only constructed a "naive" Ryser algorithm without the gray ordering

using Base.Threads

function naive_ryser(A::AbstractMatrix{T}) where T
    """ Ryser algorithm without gray ordering """
    res = Atomic{T}(0)
    n,m = size(A)
    n==m || throw(ArgumentError("perm: argument must be a square matrix"))
    m=2^(n-1)
    @inbounds @threads for d=0:m-1
        v = one(T)
        for j=1:n
            w = A[1,j]
            for i=2:n
                if d & 1<<(i-2) != 0
                    w += A[i,j]
                else
                    w -= A[i,j]
                end
            end
            v *= w
        end
        if iseven(count_ones(d))
            atomic_sub!(res,v)
        else
            atomic_add!(res,v)
        end
    end
    res[]/m
end

this code has complexity O(2^n n^2) rather than O(2^n n), but in my machine with a 30x30 matrix it is "only" 8 times slower (rather than 30 times). This shows the advantage of multi-threading.

To make a proper code with proper scaling for computing the permanent, we should use gray ordering. The difficulty is to parallelize the sum over gray sequences. One possibility would be to divide the sum over gray patterns into a fixed number of batches, where the sums in the different batches are done by different threads. In each thread, we can use the same code of ryser, but with a proper initial value of v. However, computing the initial value of v might be complex and kill all the benefits of multi-threading. Other ideas?

benoitseron commented 2 years ago

In fact we discussed this just a few hours ago with @AntoineRestivo. We can certainly take inspiration from https://arxiv.org/abs/1904.06229

Otherwise the case of the tensor permanent may give even more useful results.

benkj commented 2 years ago

Great, that's exactly what I had in mind.

Il gio 17 feb 2022, 20:48 Benoît Seron @.***> ha scritto:

In fact we discussed this just a few hours ago with @AntoineRestivo https://github.com/AntoineRestivo. We can certainly take inspiration from https://arxiv.org/abs/1904.06229

Otherwise the case of the tensor permanent may give even more useful results.

— Reply to this email directly, view it on GitHub https://github.com/benoitseron/Permanents.jl/issues/6#issuecomment-1043355868, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD7GTCTH3EUB3WUHBXNZUDTU3VGKJANCNFSM5OVOS63A . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>