Boolean matrix multiplication in Julia

Following Oscar's comment of adding two for loops around your code, but without the LoopVectorization improvement, although with not allocating the full array inside the any call (so that the any stops on the first occurrence), this is decently fast (edit: replaced standard AND & with short-circuit &&):

function bool_mul2(A, B)
    mA, nA = size(A)
    mB, nB = size(B)
    nA ≠ mB && error()
    AB = BitArray(undef, mA, nB)
    for i in 1:mA, j in 1:nB
        AB[i,j] = any(A[i,k] && B[k,j] for k in 1:nA)
    end
    AB
end

(Note I removed the [ and ] inside the any to not allocate there.

E.g., with A and B of size 1000×1000, I get

julia> @btime bool_mul2($A, $B) ;
  16.128 ms (3 allocations: 122.25 KiB)

compared to

julia> @btime bool_mul($A, $B) ;
  346.374 ms (12 allocations: 7.75 MiB)

---

EDIT: For squaring the matrix, maybe try

function bool_square(A)
    m, n = size(A)
    m ≠ n && error()
    A² = BitArray(undef, n, n)
    for i in 1:n, j in 1:n
        A²[i,j] = any(A[i,k] && A[k,j] for k in 1:n)
    end
    A²
end

for which I get

julia> A = rand(Bool, 500, 500) ;

julia> @btime $A * $A .!= 0 ;
  42.483 ms (12 allocations: 1.94 MiB)

julia> @btime bool_square($A) ;
  4.653 ms (3 allocations: 30.69 KiB)

One very simple solution is

function bool_mul(A,B)
    return A*B .!= 0
end

This won't be the most efficient since it will allocate a matrix for A*B, but might end up being one of the fastest solutions available.