From cb21b8da3d0f434a23904d8f77baee983e9d2a91 Mon Sep 17 00:00:00 2001 From: Sungho Shin Date: Mon, 12 Oct 2026 23:53:25 -0400 Subject: [PATCH] ordering: GPU-hybrid METIS NodeND bench with packed-sort contraction and task-parallel recursion Adds bench/ordering (the exact METIS 5.1.0 NodeND port, the device matching and contraction kernels, and the hybrid ordering developed alongside #126) with two performance changes: the contraction groups edges with one radix sort! of (crow<<44 | ccol<<24 | edge index) instead of sortperm + gathers (element-identical output, 1.8x on the GV100 coarsening chain, guarded with the sortperm path as fallback for larger graphs), and the sub-threshold recursion spawns tasks above SPAWN_MIN instead of running one serial walk per top subtree. Opt-in PARALLEL_TOP overlaps the top sibling bisections' host work; on multi-GPU hosts TASK_SETUP must pin each task's device (documented). Case78484 (n = 674,562, GV100, 16 Julia threads): MadNLP init_time 8.62s -> 7.42s (packed contraction) -> 5.64s (+ parallel top), identical convergence (101 iterations, same objective) and fill (1.48e7 nnz_L vs METIS 1.476e7). Gates: lap2d fill 1.017x vs the exact port; 78k fill 1.003x/1.012x, valid permutations. The recursion gains require a multi-threaded Julia session; single-threaded sessions should keep algo4. Co-Authored-By: Claude Fable 5 --- bench/ordering/GpuBatch.jl | 177 +++ bench/ordering/GpuND.jl | 399 ++++++ bench/ordering/GpuRefine.jl | 267 ++++ bench/ordering/Hybrid.jl | 339 +++++ bench/ordering/MetisExact.jl | 2411 ++++++++++++++++++++++++++++++++ bench/ordering/gpu_nodend.jl | 31 + bench/ordering/nodend_smoke.jl | 51 + 7 files changed, 3675 insertions(+) create mode 100644 bench/ordering/GpuBatch.jl create mode 100644 bench/ordering/GpuND.jl create mode 100644 bench/ordering/GpuRefine.jl create mode 100644 bench/ordering/Hybrid.jl create mode 100644 bench/ordering/MetisExact.jl create mode 100644 bench/ordering/gpu_nodend.jl create mode 100644 bench/ordering/nodend_smoke.jl diff --git a/bench/ordering/GpuBatch.jl b/bench/ordering/GpuBatch.jl new file mode 100644 index 0000000..0861c1d --- /dev/null +++ b/bench/ordering/GpuBatch.jl @@ -0,0 +1,177 @@ +module GpuBatch + +# Level-batched sub-bisections: the ~80% of NodeND time that is thousands of +# small host bisections, re-expressed so a whole recursion level of blocks is +# ONE block-diagonal device graph, coarsened in single kernel passes. Blocks +# are disconnected components of the batch, so GpuND's matching/contraction +# apply unchanged: matching proposals follow edges and cannot cross blocks, +# and contraction preserves the block structure by construction. Per-block +# control (CoarsenTo, freezing a finished block) rides on the per-vertex +# maxvwgt vector gpu_match already takes (cap 0 = no matches = identity). +# +# The 2026-10-09 Batched.jl attempt is the documented dead-end this replaces: +# it re-sliced every block to a host GraphT at every level (12x regression). +# Here the batch stays device-resident across levels; the host sees only +# per-block counters (nb Int32s per level). + +using SparseDirectSolver.KernelAbstractions +const KA = KernelAbstractions +using SparseDirectSolver.Atomix +using ..MetisExact, ..GpuND +const ME = MetisExact + +export build_batch, coarsen_batch!, batch_invariants + +struct Batch + n::Int # total vertices + nb::Int # blocks + dx; da; dw; dv # device CSR, GpuND conventions (0-based ids) + bid # device Int32, 1-based block id per vertex + nv_b::Vector{Int} # host per-block vertex counts + tvw_b::Vector{Int} # host per-block total vertex weight +end + +""" + build_batch(backend, graphs::Vector{ME.GraphT}) -> Batch + +Concatenate independent subgraphs into one block-diagonal device CSR +(vertex ids shifted; no edges are added, so blocks stay disconnected). +""" +function build_batch(backend, graphs::Vector{ME.GraphT}) + nb = length(graphs) + n = sum(g.nvtxs for g in graphs) + ne = sum(g.nedges for g in graphs) + xadj = Vector{Int32}(undef, n + 1) + adj = Vector{Int32}(undef, ne) + adjw = Vector{Int32}(undef, ne) + vw = Vector{Int32}(undef, n) + bid = Vector{Int32}(undef, n) + nv_b = Vector{Int}(undef, nb) + tvw_b = zeros(Int, nb) + xadj[1] = 0 + off = 0 + voff = 0 + for (b, g) in enumerate(graphs) + nv_b[b] = g.nvtxs + for i in 1:g.nvtxs + vw[voff + i] = Int32(g.vwgt[i]) + tvw_b[b] += g.vwgt[i] + bid[voff + i] = Int32(b) + xadj[voff + i + 1] = Int32(off + g.xadj[i + 1]) + end + for p in 1:g.nedges + adj[off + p] = Int32(g.adjncy[p] + voff) + adjw[off + p] = Int32(g.adjwgt[p]) + end + off += g.nedges + voff += g.nvtxs + end + dev(x) = (y = KA.allocate(backend, Int32, length(x)); copyto!(y, x); y) + Batch(n, nb, dev(xadj), dev(adj), dev(adjw), dev(vw), dev(bid), nv_b, tvw_b) +end + +@kernel function _mv_kernel!(mv, @Const(bid), @Const(caps)) + v = @index(Global) + @inbounds mv[v] = caps[bid[v]] +end + +# every fine vertex of a coarse vertex carries the same block id +@kernel function _bidmap_kernel!(bidc, @Const(bid), @Const(cmap)) + v = @index(Global) + @inbounds bidc[cmap[v] + 1] = bid[v] +end + +@kernel function _bcount_kernel!(cnt, @Const(bid)) + v = @index(Global) + @inbounds Atomix.@atomic cnt[bid[v]] += Int32(1) +end + +@kernel function _bwsum_kernel!(ws, @Const(bid), @Const(vwgt)) + v = @index(Global) + @inbounds Atomix.@atomic ws[bid[v]] += Int32(vwgt[v]) +end + +@kernel function _crossedge_kernel!(bad, @Const(xadj), @Const(adjncy), @Const(bid)) + v = @index(Global) + @inbounds for p in (xadj[v] + 1):xadj[v + 1] + bid[adjncy[p] + 1] != bid[v] && (bad[1] = Int32(1)) + end +end + +function _percount(backend, bid, nb, n) + cnt = KA.allocate(backend, Int32, nb) + fill!(cnt, Int32(0)) + _bcount_kernel!(backend)(cnt, bid; ndrange = n) + KA.synchronize(backend) + return Int.(Array(cnt)) +end + +""" + batch_invariants(backend, B) -> (cross_block_edges::Bool, nv, wsum) + +Structural checks: no edge crosses blocks; per-block vertex counts and weight +sums (the weight sums must be invariant under every contraction level). +""" +function batch_invariants(backend, B::Batch) + bad = KA.allocate(backend, Int32, 1) + fill!(bad, Int32(0)) + _crossedge_kernel!(backend)(bad, B.dx, B.da, B.bid; ndrange = B.n) + ws = KA.allocate(backend, Int32, B.nb) + fill!(ws, Int32(0)) + _bwsum_kernel!(backend)(ws, B.bid, B.dv; ndrange = B.n) + KA.synchronize(backend) + return (Array(bad)[1] != 0, _percount(backend, B.bid, B.nb, B.n), Int.(Array(ws))) +end + +""" + coarsen_batch!(backend, B; seed = 1, levels = nothing, maxlevels = 40) + -> (B_coarse, frozen_done) + +Coarsen every block of the batch simultaneously until each reaches its METIS +CoarsenTo (clamp(nvtxs/8, 40, 100)) or stalls (shrink > COARSEN_FRACTION), +freezing finished blocks via a zero matching cap. `levels`, when a vector, is +filled with (dx, da, dv, n, cmap, bid) per level, finest first, for device +uncoarsening. +""" +function coarsen_batch!(backend, B::Batch; seed::Integer = 1, + levels::Union{Nothing, Vector{Any}} = nothing, + maxlevels::Int = 40) + nb = B.nb + ct_b = [clamp(B.nv_b[b] ÷ 8, 40, 100) for b in 1:nb] + frozen = [B.nv_b[b] <= ct_b[b] for b in 1:nb] + nv = copy(B.nv_b) + dx, da, dw, dv, bid, n = B.dx, B.da, B.dw, B.dv, B.bid, B.n + caps = KA.allocate(backend, Int32, nb) + caps_h = Vector{Int32}(undef, nb) + level = 0 + while !all(frozen) && level < maxlevels + for b in 1:nb + caps_h[b] = frozen[b] ? Int32(0) : + Int32(max(trunc(Int, 1.5 * B.tvw_b[b] / ct_b[b]), 1)) + end + copyto!(caps, caps_h) + mv = KA.allocate(backend, Int32, n) + _mv_kernel!(backend)(mv, bid, caps; ndrange = n) + KA.synchronize(backend) + match = GpuND.gpu_match(backend, dx, da, dw, dv, mv; seed = seed + level) + cn, cx, ca, cw, cv, cmap = GpuND.gpu_contract(backend, dx, da, dw, dv, match) + levels === nothing || + push!(levels, (dx = dx, da = da, dv = dv, n = n, cmap = cmap, bid = bid)) + bidc = KA.allocate(backend, Int32, cn) + _bidmap_kernel!(backend)(bidc, bid, cmap; ndrange = n) + KA.synchronize(backend) + nvc = _percount(backend, bidc, nb, cn) + for b in 1:nb + frozen[b] && continue + if nvc[b] <= ct_b[b] || nvc[b] > ME.COARSEN_FRACTION * nv[b] + frozen[b] = true # done, or stalled at this size + end + end + nv = nvc + dx, da, dw, dv, bid, n = cx, ca, cw, cv, bidc, cn + level += 1 + end + return (Batch(n, nb, dx, da, dw, dv, bid, nv, copy(B.tvw_b)), frozen) +end + +end # module diff --git a/bench/ordering/GpuND.jl b/bench/ordering/GpuND.jl new file mode 100644 index 0000000..8d120ad --- /dev/null +++ b/bench/ordering/GpuND.jl @@ -0,0 +1,399 @@ +module GpuND + +# M2: device kernels for the NodeND pipeline, written on KernelAbstractions +# against the exact host port (MetisExact.jl) as the stage oracle. +# +# Parallelization contract, stage by stage: +# matching — hashed-handshake HEM: heaviest-edge preference, maxvwgt +# respected, per-round hash tie-breaking. Deterministic given +# (graph, seed, rounds) but NOT the sequential SHEM order, so +# it is validated by invariants + end quality, not bitwise. +# contraction — given a matching it is order-free: validated ROW-SET +# identical (sorted (id, wgt) rows) against the host +# CreateCoarseGraph fed the same match array. +# edge order — GPU coarse rows are sorted by coarse id (canonical), +# where METIS's are hash-insertion-ordered; downstream GPU +# stages define their own deterministic order. +# +# All indices 0-based in device arrays, matching the host port. + +using SparseDirectSolver.KernelAbstractions +const KA = KernelAbstractions + +export gpu_match, gpu_contract + +# --------------------------------------------------------------------------- +# matching: propose / resolve rounds + +@kernel function propose_kernel!(prop, @Const(xadj), @Const(adjncy), @Const(adjwgt), + @Const(vwgt), @Const(match), @Const(maxvwgt), salt::UInt64) + gi = @index(Global) + v0 = gi - 1 # 0-based vertex + v = Int32(v0) + if match[v + 1] == Int32(-2) # -2 = unmatched (device sentinel) + best = Int32(-1) + bw = Int32(-1) + bh = UInt64(0) + mv = maxvwgt[v + 1] # per-vertex cap (batched: per subgraph) + if mv > Int32(0) && vwgt[v + 1] < mv + for p in (xadj[v + 1] + 1):xadj[v + 2] + u = adjncy[p] + w = adjwgt[p] + if u != v && match[u + 1] == Int32(-2) && vwgt[v + 1] + vwgt[u + 1] <= mv + hu = (UInt64(u) * 0x9e3779b97f4a7c15 + salt) >> 1 + if w > bw || (w == bw && (best == Int32(-1) || hu > bh)) + best = u + bw = w + bh = hu + end + end + end + end + prop[v + 1] = best + else + prop[v + 1] = Int32(-1) + end +end + +@kernel function resolve_kernel!(match, @Const(prop)) + gi = @index(Global) + v0 = gi - 1 + v = Int32(v0) + u = prop[v + 1] + if u != Int32(-1) && prop[u + 1] == v && match[v + 1] == Int32(-2) + # mutual proposal; both sides see it symmetrically, no atomics needed + match[v + 1] = u + end +end + +""" + gpu_match(backend, xadj, adjncy, adjwgt, vwgt, maxvwgt; rounds = 8, seed = 1) -> match + +Heavy-edge matching by hashed mutual proposals. `match[v+1]` is the partner +(0-based) or `v` for a singleton. Input arrays are device vectors (Int32, +0-based ids); the result is a device vector. +""" +function gpu_match(backend, xadj, adjncy, adjwgt, vwgt, maxvwgt; + rounds::Int = 8, seed::Integer = 1) + n = length(xadj) - 1 + if maxvwgt isa Integer # scalar cap -> per-vertex vector + mv = KA.allocate(backend, Int32, n) + fill!(mv, Int32(maxvwgt)) + maxvwgt = mv + end + match = KA.allocate(backend, Int32, n) + fill!(match, Int32(-2)) + prop = KA.allocate(backend, Int32, n) + pk = propose_kernel!(backend) + rk = resolve_kernel!(backend) + for r in 1:rounds + salt = UInt64(seed) * 0xbf58476d1ce4e5b9 + UInt64(r) * 0x94d049bb133111eb + pk(prop, xadj, adjncy, adjwgt, vwgt, match, maxvwgt, salt; ndrange = n) + rk(match, prop; ndrange = n) + end + # 2-hop matching for hub-dense graphs (METIS Match_2HopAny semantics): + # unmatched low-degree vertices pair through a common neighbour. Parallel + # form: per-neighbour atomic-min "two smallest unmatched" registers, then + # the same mutual-proposal resolve as the handshake. Gated like METIS at + # >10% unmatched. + nun = count_unmatched(backend, match) + if nun > n ÷ 10 + minu = KA.allocate(backend, Int32, n) + min2u = KA.allocate(backend, Int32, n) + for (rnds, maxdeg) in ((2, 2), (2, 3), (3, typemax(Int32))) + for _ in 1:rnds + fill!(minu, typemax(Int32)) + fill!(min2u, typemax(Int32)) + minu_kernel!(backend)(minu, xadj, adjncy, match, maxvwgt, Int32(maxdeg); ndrange = n) + min2_kernel!(backend)(min2u, minu, xadj, adjncy, match, maxvwgt, Int32(maxdeg); ndrange = n) + twohop_propose_kernel!(backend)(prop, minu, min2u, xadj, adjncy, match, maxvwgt, + Int32(maxdeg); ndrange = n) + resolve_kernel!(backend)(match, prop; ndrange = n) + end + nun = count_unmatched(backend, match) + nun > n ÷ 10 || break + end + end + # leftovers become singletons [KERNEL] + singles_kernel!(backend)(match; ndrange = n) + KA.synchronize(backend) + return match +end + +function count_unmatched(backend, match) + f = KA.allocate(backend, Int32, length(match)) + unmatched_flag_kernel!(backend)(f, match; ndrange = length(match)) + return Int(sum(f)) # the reduction synchronizes +end + +@kernel function unmatched_flag_kernel!(f, @Const(match)) + gi = @index(Global) + f[gi] = match[gi] == Int32(-2) ? Int32(1) : Int32(0) +end + +@kernel function minu_kernel!(minu, @Const(xadj), @Const(adjncy), @Const(match), @Const(maxvwgt), maxdeg::Int32) + gi = @index(Global) + v0 = gi - 1 + v = Int32(v0) + if match[v + 1] == Int32(-2) && maxvwgt[v + 1] > Int32(0) && xadj[v + 2] - xadj[v + 1] < maxdeg + for p in (xadj[v + 1] + 1):xadj[v + 2] + c = adjncy[p] + KA.@atomic minu[c + 1] min v + end + end +end + +@kernel function min2_kernel!(min2u, @Const(minu), @Const(xadj), @Const(adjncy), + @Const(match), @Const(maxvwgt), maxdeg::Int32) + gi = @index(Global) + v0 = gi - 1 + v = Int32(v0) + if match[v + 1] == Int32(-2) && maxvwgt[v + 1] > Int32(0) && xadj[v + 2] - xadj[v + 1] < maxdeg + for p in (xadj[v + 1] + 1):xadj[v + 2] + c = adjncy[p] + if minu[c + 1] != v + KA.@atomic min2u[c + 1] min v + end + end + end +end + +@kernel function twohop_propose_kernel!(prop, @Const(minu), @Const(min2u), @Const(xadj), + @Const(adjncy), @Const(match), @Const(maxvwgt), maxdeg::Int32) + gi = @index(Global) + v0 = gi - 1 + v = Int32(v0) + best = Int32(-1) + if match[v + 1] == Int32(-2) && maxvwgt[v + 1] > Int32(0) && xadj[v + 2] - xadj[v + 1] < maxdeg + for p in (xadj[v + 1] + 1):xadj[v + 2] + c = adjncy[p] + u = minu[c + 1] == v ? min2u[c + 1] : minu[c + 1] + if u != typemax(Int32) && u != v && match[u + 1] == Int32(-2) + (best == Int32(-1) || u < best) && (best = u) + end + end + end + prop[v + 1] = best +end + +@kernel function singles_kernel!(match) + gi = @index(Global) + v0 = gi - 1 + if match[v0 + 1] == Int32(-2) + match[v0 + 1] = Int32(v0) + end +end + +# --------------------------------------------------------------------------- +# contraction: coarse ids by scan over pair representatives, then +# pack (coarse row, coarse col) into UInt64 keys, sort, reduce duplicates. + +@kernel function repflag_kernel!(flag, @Const(match)) + gi = @index(Global) + v0 = gi - 1 + flag[v0 + 1] = Int32(v0) <= match[v0 + 1] ? Int32(1) : Int32(0) +end + +@kernel function cmap_kernel!(cmap, @Const(match), @Const(repid)) + gi = @index(Global) + v0 = gi - 1 + u = match[v0 + 1] + if Int32(v0) <= u + cmap[v0 + 1] = repid[v0 + 1] + else + cmap[v0 + 1] = repid[u + 1] + end +end + +@kernel function packkeys_kernel!(keys, vals, @Const(xadj), @Const(adjncy), + @Const(adjwgt), @Const(cmap)) + gi = @index(Global) + v0 = gi - 1 + cv = cmap[v0 + 1] + for p in (xadj[v0 + 1] + 1):xadj[v0 + 2] + cu = cmap[adjncy[p] + 1] + if cu == cv + keys[p] = typemax(UInt64) # self edge: sorts to the end + vals[p] = Int32(0) + else + keys[p] = (UInt64(cv) << 32) | UInt64(cu) + vals[p] = adjwgt[p] + end + end +end + +@kernel function boundary_kernel!(isnew, @Const(keys)) + i = @index(Global) + k = keys[i] + if k == typemax(UInt64) + isnew[i] = Int32(0) + elseif i == 1 + isnew[i] = Int32(1) + else + isnew[i] = keys[i - 1] == k ? Int32(0) : Int32(1) + end +end + +@kernel function cvwgt_kernel!(cvwgt, @Const(vwgt), @Const(match), @Const(cmap)) + gi = @index(Global) + v0 = gi - 1 + u = match[v0 + 1] + if Int32(v0) <= u + w = vwgt[v0 + 1] + u != Int32(v0) && (w += vwgt[u + 1]) + cvwgt[cmap[v0 + 1] + 1] = w + end +end + +@kernel function scatter_edges_kernel!(cadjncy, cadjwgt, rowcount, @Const(keys), + @Const(segsum), @Const(isnew), @Const(pos)) + i = @index(Global) + if isnew[i] == Int32(1) + j = pos[i] # output slot (1-based) + k = keys[i] + cadjncy[j] = Int32(k & 0xffffffff) + cadjwgt[j] = segsum[i] + KA.@atomic rowcount[Int(k >> 32) + 2] += Int32(1) + end +end + +""" + gpu_contract(backend, xadj, adjncy, adjwgt, vwgt, match) + -> (cn, cxadj, cadjncy, cadjwgt, cvwgt, cmap) + +Contract a matching. Coarse rows come out sorted by coarse id; fully +device-resident, the representative scan and segment scans are device cumsums, +segmented reduction by Int32 atomics (deterministic: integer adds commute). + +Two grouping paths, element-identical output: + +* packed (the common case, `cn < 2^20` and `m < 2^24`): one radix `sort!` of + `(crow << 44) | (ccol << 24) | edge index` — the free low bits carry the + edge index, so no `sortperm` and no gather; the weight is read back through + the embedded index. `sort!` of a plain `UInt64` vector is ~3x faster than + `sortperm` on the same keys (GV100, 7.3M edges: 50ms vs 167ms), which makes + the whole contraction ~1.8x faster. +* fallback (larger graphs): the original `sortperm` over `(crow << 32) | ccol` + keys with a value gather. +""" +function gpu_contract(backend, xadj, adjncy, adjwgt, vwgt, match) + n = length(xadj) - 1 + m = length(adjncy) + flag = KA.allocate(backend, Int32, n) + repflag_kernel!(backend)(flag, match; ndrange = n) + KA.synchronize(backend) + incl = cumsum(flag) # inclusive scan (GPUArrays / Base) + cn = Int(Array(view(incl, n:n))[1]) + repid = KA.allocate(backend, Int32, n) + exclusive_kernel!(backend)(repid, incl, flag; ndrange = n) + cmap = KA.allocate(backend, Int32, n) + cmap_kernel!(backend)(cmap, match, repid; ndrange = n) + local isnew, segid, nseg, cadjncy, cadjwgt, rowcount + if cn < (1 << 20) && m < (1 << 24) + packed = KA.allocate(backend, UInt64, m) + packed_kernel!(backend)(packed, xadj, adjncy, cmap; ndrange = n) + KA.synchronize(backend) + sort!(packed) # plain radix sort, no permutation + isnew = KA.allocate(backend, Int32, m) + boundary_packed_kernel!(backend)(isnew, packed; ndrange = m) + KA.synchronize(backend) + segid = cumsum(isnew) # segment index per element (1-based) + nseg = Int(Array(view(segid, m:m))[1]) + cadjncy = KA.allocate(backend, Int32, nseg) + cadjwgt = KA.allocate(backend, Int32, nseg) + fill!(cadjwgt, Int32(0)) + rowcount = KA.allocate(backend, Int32, Int(cn) + 1) + fill!(rowcount, Int32(0)) + segfill_packed_kernel!(backend)(cadjncy, cadjwgt, rowcount, packed, adjwgt, + isnew, segid; ndrange = m) + else + keys = KA.allocate(backend, UInt64, m) + vals = KA.allocate(backend, Int32, m) + packkeys_kernel!(backend)(keys, vals, xadj, adjncy, adjwgt, cmap; ndrange = n) + KA.synchronize(backend) + p = sortperm(keys) + keys = keys[p] + vals = vals[p] + isnew = KA.allocate(backend, Int32, m) + boundary_kernel!(backend)(isnew, keys; ndrange = m) + KA.synchronize(backend) + segid = cumsum(isnew) + nseg = Int(Array(view(segid, m:m))[1]) + cadjncy = KA.allocate(backend, Int32, nseg) + cadjwgt = KA.allocate(backend, Int32, nseg) + fill!(cadjwgt, Int32(0)) + rowcount = KA.allocate(backend, Int32, Int(cn) + 1) + fill!(rowcount, Int32(0)) + segfill_kernel!(backend)(cadjncy, cadjwgt, rowcount, keys, vals, isnew, segid; ndrange = m) + end + KA.synchronize(backend) + cxadj = cumsum(rowcount) # rowcount[1] = 0 stays 0 + cvwgt = KA.allocate(backend, Int32, Int(cn)) + cvwgt_kernel!(backend)(cvwgt, vwgt, match, cmap; ndrange = n) + KA.synchronize(backend) + return Int(cn), cxadj, cadjncy, cadjwgt, cvwgt, cmap +end + +# packed-path kernels: the key is (crow << 44) | (ccol << 24) | edge index, +# self edges sort to the end as typemax +@kernel function packed_kernel!(packed, @Const(xadj), @Const(adjncy), @Const(cmap)) + gi = @index(Global) + v0 = gi - 1 + cv = cmap[v0 + 1] + for p in (xadj[v0 + 1] + 1):xadj[v0 + 2] + cu = cmap[adjncy[p] + 1] + if cu == cv + packed[p] = typemax(UInt64) + else + packed[p] = (UInt64(cv) << 44) | (UInt64(cu) << 24) | UInt64(p - 1) + end + end +end + +@kernel function boundary_packed_kernel!(isnew, @Const(packed)) + i = @index(Global) + k = packed[i] + if k == typemax(UInt64) + isnew[i] = Int32(0) + elseif i == 1 + isnew[i] = Int32(1) + else + isnew[i] = (packed[i - 1] >> 24) == (k >> 24) ? Int32(0) : Int32(1) + end +end + +@kernel function segfill_packed_kernel!(cadjncy, cadjwgt, rowcount, @Const(packed), + @Const(adjwgt), @Const(isnew), @Const(segid)) + i = @index(Global) + k = packed[i] + if k != typemax(UInt64) + sid = segid[i] + idx = Int(k & 0xffffff) # original edge position (0-based) + KA.@atomic cadjwgt[sid] += adjwgt[idx + 1] + if isnew[i] == Int32(1) + cadjncy[sid] = Int32((k >> 24) & 0xfffff) + KA.@atomic rowcount[Int(k >> 44) + 2] += Int32(1) + end + end +end + +@kernel function exclusive_kernel!(excl, @Const(incl), @Const(flag)) + gi = @index(Global) + excl[gi] = incl[gi] - flag[gi] +end + +@kernel function segfill_kernel!(cadjncy, cadjwgt, rowcount, @Const(keys), @Const(vals), + @Const(isnew), @Const(segid)) + gi = @index(Global) + k = keys[gi] + if k != typemax(UInt64) + sid = segid[gi] + KA.@atomic cadjwgt[sid] += vals[gi] + if isnew[gi] == Int32(1) + cadjncy[sid] = Int32(k & 0xffffffff) + KA.@atomic rowcount[Int(k >> 32) + 2] += Int32(1) + end + end +end + +end # module diff --git a/bench/ordering/GpuRefine.jl b/bench/ordering/GpuRefine.jl new file mode 100644 index 0000000..16ca1b2 --- /dev/null +++ b/bench/ordering/GpuRefine.jl @@ -0,0 +1,267 @@ +module GpuRefine + +# Device node-separator refinement: a round-based, data-parallel approximation +# of METIS's FM_2WayNodeRefine1Sided (MetisExact.jl:1911), so the uncoarsening +# walk of a bisection can stay on the device instead of round-tripping to the +# host FM at every level. This is the prerequisite for level-batched +# sub-bisections (the Batched.jl dead-end failed on exactly these host trips). +# +# Semantics relative to the serial FM: +# * a PASS moves separator vertices to ONE side (`to`), alternating sides +# across passes, exactly like the 1-sided FM; +# * within a pass, moves are selected in parallel: every separator vertex +# with gain = vwgt[v] - Σ_{u∈N(v), where[u]==other} vwgt[u] ≥ 0 moves, +# subject to a balance budget (badmaxpwgt, METIS's cap). The serial code's +# priority order, incremental rinfo and negative-gain hill-climb with +# best-prefix rollback are replaced by whole-pass accept/revert: the pass +# is kept only if it reduces the separator weight (ties: improves balance). +# * separator validity is preserved for ANY subset of movers: all movers go +# to the same side, and every other-side neighbor of a mover is pulled +# into the separator in the same pass, so no side-to-side edge can appear. +# +# Known approximation costs (accepted; the quality gate is end nnz_L): +# * concurrent movers overcount shared pulled neighbors in their gains; +# * no negative-gain hill-climbing; +# * under a binding balance budget the atomic claim makes the selected +# subset nondeterministic (the budget rarely binds in the ND chain). + +using SparseDirectSolver.KernelAbstractions +const KA = KernelAbstractions +using SparseDirectSolver.Atomix + +export refine_rounds!, project_where! + +@kernel function _project_kernel!(fwhere, @Const(cwhere), @Const(cmap)) + v = @index(Global) + @inbounds fwhere[v] = cwhere[cmap[v] + 1] +end + +# gain of moving separator vertex v to `to` (other = 1 - to); non-separator +# vertices get typemin so the selection predicate is a single compare +@kernel function _gains_kernel!(gain, @Const(xadj), @Const(adjncy), @Const(vwgt), + @Const(where_), other::Int32) + v = @index(Global) + @inbounds begin + if where_[v] == Int32(2) + g = Int32(vwgt[v]) + for p in (xadj[v] + 1):xadj[v + 1] + u = adjncy[p] + 1 + where_[u] == other && (g -= vwgt[u]) + end + gain[v] = g + else + gain[v] = typemin(Int32) + end + end +end + +# candidate filter: gain ≥ 0 AND locally maximal among separator neighbors +# (ties by hashed id). Adjacent separator vertices therefore never move in +# the same pass, which removes the dominant gain-overcount (shared pulled +# neighbors of adjacent movers) and is the standard conflict rule of +# round-based refinement (Jet / mt-KaHyPar style). +@kernel function _filter_kernel!(cand, @Const(gain), @Const(xadj), @Const(adjncy), + @Const(where_), seed::UInt32, thr::Int32) + v = @index(Global) + @inbounds begin + c = Int32(0) + g = gain[v] + if g >= thr + hv = (UInt32(v) * 0x9e3779b9) ⊻ seed + c = Int32(1) + for p in (xadj[v] + 1):xadj[v + 1] + u = adjncy[p] + 1 + where_[u] == Int32(2) || continue + gu = gain[u] + gu < thr && continue + hu = (UInt32(u) * 0x9e3779b9) ⊻ seed + if gu > g || (gu == g && hu > hv) + c = Int32(0) + break + end + end + end + cand[v] = c + end +end + +# select movers among candidates under the balance budget. budget[1] starts +# at badmaxpwgt - pwgts[to]; an atomic claim keeps pwgts[to] under cap. +@kernel function _select_kernel!(mask, budget, nmov, @Const(cand), @Const(vwgt)) + v = @index(Global) + @inbounds begin + m = Int32(0) + if cand[v] == Int32(1) + left = (Atomix.@atomic budget[1] += Int32(-vwgt[v])) + if left >= Int32(0) + m = Int32(1) + Atomix.@atomic nmov[1] += Int32(1) + end + end + mask[v] = m + end +end + +@kernel function _apply_kernel!(where_, @Const(mask), to::Int32) + v = @index(Global) + @inbounds mask[v] == Int32(1) && (where_[v] = to) +end + +# pull: any `other`-side vertex adjacent to a mover enters the separator +@kernel function _pull_kernel!(where_, @Const(xadj), @Const(adjncy), @Const(mask), + other::Int32) + v = @index(Global) + @inbounds if where_[v] == other + for p in (xadj[v] + 1):xadj[v + 1] + if mask[adjncy[p] + 1] == Int32(1) + where_[v] = Int32(2) + break + end + end + end +end + +@kernel function _pwgts_kernel!(pwgts, @Const(where_), @Const(vwgt)) + v = @index(Global) + @inbounds Atomix.@atomic pwgts[where_[v] + 1] += Int32(vwgt[v]) +end + +function _pwgts!(backend, dpw, dwhere, dvw, n) + fill!(dpw, Int32(0)) + _pwgts_kernel!(backend)(dpw, dwhere, dvw; ndrange = n) + KA.synchronize(backend) + return Int.(Array(dpw)) +end + +""" + project_where!(backend, fwhere, cwhere, cmap, n) + +`fwhere[v] = cwhere[cmap[v]+1]` on the device (METIS Project2WayNodePartition +without the host parameter rebuild). +""" +function project_where!(backend, fwhere, cwhere, cmap, n::Integer) + n > 0 && _project_kernel!(backend)(fwhere, cwhere, cmap; ndrange = Int(n)) + KA.synchronize(backend) + return fwhere +end + +""" + refine_rounds!(backend, dx, da, dv, dwhere, n, ubfactors; niter = 10) + -> (pwgts, npasses) + +Round-based 1-sided node refinement of `dwhere` (0/1/2 Int32, device) against +the device CSR (`dx`, `da` 0-based ids, `dv` weights). Runs up to `2niter` +alternating passes; each pass is accepted only if the separator weight drops +(ties: balance improves), reverted otherwise; stops after an even/odd pair +without improvement, like the serial FM's termination. +""" +function refine_rounds!(backend, dx, da, dv, dwhere, n::Integer, ubfactors::Real; + niter::Int = 10, patience::Int = 4) + n == 0 && return (Int[0, 0, 0], 0) + gain = KA.allocate(backend, Int32, n) + cand = KA.allocate(backend, Int32, n) + mask = KA.allocate(backend, Int32, n) + best = KA.allocate(backend, Int32, n) + budget = KA.allocate(backend, Int32, 1) + dpw = KA.allocate(backend, Int32, 3) + pw = _pwgts!(backend, dpw, dwhere, dv, n) + badmaxpwgt = trunc(Int, 0.5 * ubfactors * (pw[1] + pw[2] + pw[3])) + to = pw[1] < pw[2] ? 1 : 0 + copyto!(best, dwhere) # global best, restored at exit + best_pw = pw + since_best = 0 + npasses = 0 + nmov = KA.allocate(backend, Int32, 1) + nmov_h = Int32[0] + for pass in 0:(2 * niter - 1) + other, to = to, (to + 1) % 2 + # a PASS is a wavefront sweep at fixed direction, like the serial + # 1-sided FM where pulled vertices join the queue mid-pass: repeat + # micro-rounds until no admissible mover remains + for micro in 0:199 + _gains_kernel!(backend)(gain, dx, da, dv, dwhere, Int32(other); ndrange = Int(n)) + _filter_kernel!(backend)(cand, gain, dx, da, dwhere, + UInt32(2654435761) + UInt32(64 * pass + micro), + Int32(0); ndrange = Int(n)) + pwt = _pwgts!(backend, dpw, dwhere, dv, n) + copyto!(budget, Int32[max(badmaxpwgt - pwt[to + 1], 0)]) + fill!(nmov, Int32(0)) + _select_kernel!(backend)(mask, budget, nmov, cand, dv; ndrange = Int(n)) + KA.synchronize(backend) + copyto!(nmov_h, nmov) + nmov_h[1] == Int32(0) && break + _apply_kernel!(backend)(dwhere, mask, Int32(to); ndrange = Int(n)) + _pull_kernel!(backend)(dwhere, dx, da, mask, Int32(other); ndrange = Int(n)) + KA.synchronize(backend) + end + pw = _pwgts!(backend, dpw, dwhere, dv, n) + npasses += 1 + # hill-climb with a global-best record, like the serial FM's rollback: + # keep walking through non-improving passes (bounded by `patience` and + # the serial code's 1.10x escape hatch), remember the best state seen + if pw[3] < best_pw[3] || + (pw[3] == best_pw[3] && abs(pw[1] - pw[2]) < abs(best_pw[1] - best_pw[2])) + copyto!(best, dwhere) + best_pw = pw + since_best = 0 + else + since_best += 1 + (since_best > patience || pw[3] > 1.10 * best_pw[3]) && break + end + end + copyto!(dwhere, best) + return (best_pw, npasses) +end + +""" + balance_rounds!(backend, dx, da, dv, dwhere, n, ubfactors) -> pwgts + +Parallel FM_2WayNodeBalance: rounds of moves toward the LIGHT side. Positive +gains move freely; when none remain and the heavy side still exceeds METIS's +cap (0.5·ubfactors·(pw1+pw2)), the maximum-gain class moves even at negative +gain. Stops when the sides cross, the cap is satisfied with no free moves, or +nothing is admissible. +""" +function balance_rounds!(backend, dx, da, dv, dwhere, n::Integer, ubfactors::Real; + maxrounds::Int = 64) + gain = KA.allocate(backend, Int32, n) + cand = KA.allocate(backend, Int32, n) + mask = KA.allocate(backend, Int32, n) + budget = KA.allocate(backend, Int32, 1) + nmov = KA.allocate(backend, Int32, 1) + nmov_h = Int32[0] + dpw = KA.allocate(backend, Int32, 3) + pw = _pwgts!(backend, dpw, dwhere, dv, n) + badmax = trunc(Int, 0.5 * ubfactors * (pw[1] + pw[2])) + max(pw[1], pw[2]) < badmax && return pw + to = pw[1] < pw[2] ? 0 : 1 + other = 1 - to + for r in 1:maxrounds + pw[to + 1] >= pw[other + 1] && break + _gains_kernel!(backend)(gain, dx, da, dv, dwhere, Int32(other); ndrange = Int(n)) + KA.synchronize(backend) + gmax = maximum(gain) + gmax == typemin(Int32) && break # no separator vertices + thr = Int32(0) + if gmax < Int32(0) + pw[other + 1] <= badmax && break # balanced enough, only bad moves left + thr = gmax # admit the best (negative) class only + end + _filter_kernel!(backend)(cand, gain, dx, da, dwhere, + UInt32(0x85ebca6b) + UInt32(r), thr; ndrange = Int(n)) + b = max(badmax - pw[to + 1], 0) + b == 0 && break + copyto!(budget, Int32[b]) + fill!(nmov, Int32(0)) + _select_kernel!(backend)(mask, budget, nmov, cand, dv; ndrange = Int(n)) + KA.synchronize(backend) + copyto!(nmov_h, nmov) + nmov_h[1] == Int32(0) && break + _apply_kernel!(backend)(dwhere, mask, Int32(to); ndrange = Int(n)) + _pull_kernel!(backend)(dwhere, dx, da, mask, Int32(other); ndrange = Int(n)) + pw = _pwgts!(backend, dpw, dwhere, dv, n) + end + return pw +end + +end # module diff --git a/bench/ordering/Hybrid.jl b/bench/ordering/Hybrid.jl new file mode 100644 index 0000000..7f01992 --- /dev/null +++ b/bench/ordering/Hybrid.jl @@ -0,0 +1,339 @@ +# Hybrid NodeND v1: the coarsening chain runs through the GpuND kernels +# (device-formulated matching + contraction), everything else through the +# bitwise-exact host port (initial separator, FM node refinement, split, MMD, +# recursion). Quality gate: nnz_L / schedule depth vs MetisExact on the same +# graphs — the cost of parallel matching is the ONLY delta by construction. +module Hybrid + +using SparseDirectSolver.KernelAbstractions +const KA = KernelAbstractions +using ..MetisExact +using ..GpuND +using Random +const ME = MetisExact + +export nodend_hybrid + +# device coarsening chain following CoarsenGraph's control flow. The device +# arrays stay resident across levels (each level's contraction output feeds +# the next level's matching directly); only the per-level host GraphT copies +# needed by the host refinement are downloaded. +function coarsen_gpu!(backend, ctrl::ME.CtrlT, graph::ME.GraphT; seed::Integer = 1, + keep_levels::Union{Nothing, Vector{Any}} = nothing) + ctrl.maxvwgt = trunc(Int, 1.5 * graph.tvwgt / ctrl.CoarsenTo) + dev(x, T) = (y = KA.allocate(backend, T, length(x)); copyto!(y, T.(x)); y) + n = graph.nvtxs + dx = dev(view(graph.xadj, 1:(n + 1)), Int32) + da = dev(view(graph.adjncy, 1:graph.nedges), Int32) + dw = dev(view(graph.adjwgt, 1:graph.nedges), Int32) + dv = dev(view(graph.vwgt, 1:n), Int32) + level = 0 + while true + match = GpuND.gpu_match(backend, dx, da, dw, dv, ctrl.maxvwgt; + seed = seed + level) + cn, cx, ca, cw, cv, cmap = GpuND.gpu_contract(backend, dx, da, dw, dv, match) + # the finer level's CSR + its fine->coarse map, for device uncoarsening + keep_levels === nothing || + push!(keep_levels, (dx = dx, da = da, dv = dv, n = graph.nvtxs, cmap = cmap)) + cgraph = ME.CreateGraph() + cgraph.nvtxs = cn + cgraph.nedges = length(ca) + cgraph.xadj = Vector{Int}(Array(cx)) + cgraph.adjncy = Vector{Int}(Array(ca)) + cgraph.adjwgt = Vector{Int}(Array(cw)) + cgraph.vwgt = Vector{Int}(Array(cv)) + ME.SetupGraph_tvwgt!(cgraph) + graph.cmap = Vector{Int}(Array(cmap)) + cgraph.finer = graph + graph.coarser = cgraph + graph = cgraph + dx, da, dw, dv = cx, ca, cw, cv # stay on the device + level += 1 + (graph.nvtxs > ctrl.CoarsenTo && + graph.nvtxs < ME.COARSEN_FRACTION * (graph.finer::ME.GraphT).nvtxs && + graph.nedges > graph.nvtxs ÷ 2) || break + end + return graph +end + +# below this size a subgraph's bisection runs entirely on the host exact +# port: small graphs pay GPU launch/transfer latency without filling the +# device, and the deep recursion has thousands of them (batched device +# coarsening of these blocks measured 4.2x SLOWER than 16 host threads on the +# case78484 population, so they stay on the CPU deliberately) +const GPU_MIN_NVTXS = Ref(50_000) + +# the sub-threshold subtree recursion is itself task-parallel: a child larger +# than SPAWN_MIN spawns a task (own CtrlT copy, task-local RNG, disjoint order +# ranges), smaller children recurse serially. Without this, each top subtree +# runs single-threaded and the imbalanced largest subtree sets the wall +# (case78484: ordering 6.3s -> 4.3s at identical fill) +const SPAWN_MIN = Ref(500) + +# also spawn the >= GPU_MIN_NVTXS siblings of the top recursion as tasks, so +# their host work (initial separator, FM, split, downloads) overlaps while the +# device coarsening serializes per device (case78484: another ~1.7x on the +# ordering at identical fill). On a MULTI-GPU host this must be combined with +# a TASK_SETUP that pins every task to one device — unpinned tasks default to +# device 0 while the caller may be elsewhere, and that cross-device mix was +# observed to corrupt the ordering. Off by default; enable per session: +# Hybrid.PARALLEL_TOP[] = true +# Hybrid.TASK_SETUP[] = let dev = CUDA.device(); () -> CUDA.device!(dev) end +const PARALLEL_TOP = Ref(false) +const TASK_SETUP = Ref{Any}(nothing) # per-task setup hook (or nothing) + +# device uncoarsening: project + round-based parallel refinement per level +# (GpuRefine.jl) instead of the host FM walk. Off by default until the fill +# quality gate passes on the benchmark set. +const GPU_REFINE = Ref(false) +# :rounds — device wavefront refinement at every level; +# :project — device projection only, one exact host balance+FM at the finest +const GPU_REFINE_MODE = Ref(:rounds) +const GPU_REFINE_TRACE = Ref(false) # print per-level separator weights + +function refine_chain_gpu!(backend, ctrl::ME.CtrlT, graph::ME.GraphT, cgraph::ME.GraphT, + levels::Vector{Any}) + # chain[k] is the graph whose CSR/cmap is levels[k]; chain[end] = cgraph + chain = ME.GraphT[graph] + while chain[end] !== cgraph + push!(chain, chain[end].coarser::ME.GraphT) + end + # coarse half: EXACT host project/balance/refine below the size threshold + # (serial FM wins where graphs are tiny and every vertex is an aggregate; + # job_35: device rounds left the n=108 level at sides 90k/55k) + k = length(chain) - 1 + while k >= 1 && chain[k].nvtxs < REFINE_HOST_BELOW[] + g = chain[k] + ME.Project2WayNodePartition!(ctrl, g) + ME.FM_2WayNodeBalance!(ctrl, g) + ME.FM_2WayNodeRefine1Sided!(ctrl, g, ctrl.niter) + GPU_REFINE_TRACE[] && + println(" host n=", g.nvtxs, " sep=", g.pwgts[3], " sides=", + g.pwgts[1], "/", g.pwgts[2]) + k -= 1 + end + k == 0 && return nothing # whole chain refined on host + # fine half: device projection + balance + wavefront rounds + src = chain[k + 1] + cw = KA.allocate(backend, Int32, src.nvtxs) + copyto!(cw, Int32.(view(src.where_, 1:src.nvtxs))) + for j in k:-1:1 + L = levels[j] + fw = KA.allocate(backend, Int32, L.n) + Main.GpuRefine.project_where!(backend, fw, cw, L.cmap, L.n) + if GPU_REFINE_MODE[] == :rounds + Main.GpuRefine.balance_rounds!(backend, L.dx, L.da, L.dv, fw, L.n, + Float64(ctrl.ubfactors)) + pw, np = Main.GpuRefine.refine_rounds!(backend, L.dx, L.da, L.dv, fw, L.n, + Float64(ctrl.ubfactors); niter = ctrl.niter) + GPU_REFINE_TRACE[] && + println(" dev n=", L.n, " sep=", pw[3], " sides=", pw[1], "/", pw[2], + " passes=", np) + end + cw = fw + end + ME.Allocate2WayNodePartitionMemory!(graph) + copyto!(view(graph.where_, 1:graph.nvtxs), Int.(Array(cw))) + graph.coarser = nothing + ME.Compute2WayNodePartitionParams!(ctrl, graph) + if GPU_REFINE_MODE[] == :project + ME.FM_2WayNodeBalance!(ctrl, graph) + ME.FM_2WayNodeRefine1Sided!(ctrl, graph, ctrl.niter) + end + return nothing +end + +# device refinement runs only at levels at least this large; smaller levels +# use the exact serial walk above (they cost microseconds on the host) +const REFINE_HOST_BELOW = Ref(50_000) + +# L1-shaped bisection with the device coarsening +function node_bisection!(backend, ctrl::ME.CtrlT, graph::ME.GraphT, niparts::Int) + if graph.nvtxs < GPU_MIN_NVTXS[] + ME.MlevelNodeBisectionL2!(ctrl, graph, niparts) + return nothing + end + ctrl.CoarsenTo = graph.nvtxs ÷ 8 + ctrl.CoarsenTo > 100 && (ctrl.CoarsenTo = 100) + ctrl.CoarsenTo < 40 && (ctrl.CoarsenTo = 40) + levels = GPU_REFINE[] ? Any[] : nothing + cgraph = coarsen_gpu!(backend, ctrl, graph; keep_levels = levels) + niparts = max(1, cgraph.nvtxs <= ctrl.CoarsenTo ? niparts ÷ 2 : niparts) + ME.InitSeparator!(ctrl, cgraph, niparts) + if levels === nothing + ME.Refine2WayNode!(ctrl, graph, cgraph) + else + refine_chain_gpu!(backend, ctrl, graph, cgraph, levels) + end + return nothing +end + +# MlevelNodeBisectionMultiple semantics: nseps tries, keep the lightest separator +function node_bisection_multiple!(backend, ctrl::ME.CtrlT, graph::ME.GraphT) + if ctrl.nseps == 1 || graph.nvtxs < (ctrl.compress ? 1000 : 2000) + node_bisection!(backend, ctrl, graph, ME.LARGENIPARTS) + return nothing + end + bestwhere = Vector{Int}(undef, graph.nvtxs) + mincut = graph.tvwgt + for i in 0:(ctrl.nseps - 1) + node_bisection!(backend, ctrl, graph, ME.LARGENIPARTS) + if i == 0 || graph.mincut < mincut + mincut = graph.mincut + i < ctrl.nseps - 1 && copyto!(bestwhere, graph.where_) + end + mincut == 0 && break + end + if mincut != graph.mincut + copyto!(graph.where_, bestwhere) + ME.Compute2WayNodePartitionParams!(ctrl, graph) + end + return nothing +end + +_copy_ctrl(c::ME.CtrlT) = ME.CtrlT(c.ctype, c.rtype, c.iptype, c.nseps, c.niter, c.ufactor, + c.compress, c.ccorder, c.pfactor, c.no2hop, c.seed, + c.CoarsenTo, c.maxvwgt, c.ubfactors, copy(c.pijbm)) + +# host nested dissection of a sub-threshold subtree, task-parallel above +# SPAWN_MIN: order ranges are disjoint by construction, every spawned task gets +# its own CtrlT copy and a seeded task-local RNG (the invariants the existing +# top-level spawns already rely on) +function mlnd!(ctrl::ME.CtrlT, graph::ME.GraphT, order::Vector{Int}, lastvtx::Int) + ME.MlevelNodeBisectionMultiple!(ctrl, graph) + nbnd = graph.nbnd + bndind = graph.bndind + label = graph.label + for i in 0:(nbnd - 1) + lastvtx -= 1 + order[label[bndind[i + 1] + 1] + 1] = lastvtx + end + lgraph, rgraph = ME.SplitGraphOrder(ctrl, graph) + tasks = Task[] + for (g, lv) in ((lgraph, lastvtx - rgraph.nvtxs), (rgraph, lastvtx)) + if g.nvtxs > ME.MMDSWITCH && g.nedges > 0 + if g.nvtxs > SPAWN_MIN[] && Threads.nthreads() > 1 + c = _copy_ctrl(ctrl) + let g = g, lv = lv, c = c + push!(tasks, Threads.@spawn begin + Random.seed!(Random.TaskLocalRNG(), UInt32(hash((g.nvtxs, lv)) % UInt32)) + mlnd!(c, g, order, lv) + end) + end + else + mlnd!(ctrl, g, order, lv) # serial: ctrl mutation is sequential + end + else + ME.MMDOrder!(ctrl, g, order, lv) + end + end + foreach(wait, tasks) + return nothing +end + +function dissect!(backend, ctrl::ME.CtrlT, graph::ME.GraphT, order::Vector{Int}, lastvtx::Int) + node_bisection_multiple!(backend, ctrl, graph) + nbnd = graph.nbnd + bndind = graph.bndind + label = graph.label + for i in 0:(nbnd - 1) + lastvtx -= 1 + order[label[bndind[i + 1] + 1] + 1] = lastvtx + end + lgraph, rgraph = ME.SplitGraphOrder(ctrl, graph) + tasks = Task[] + for (g, lv) in ((lgraph, lastvtx - rgraph.nvtxs), (rgraph, lastvtx)) + if g.nvtxs > ME.MMDSWITCH && g.nedges > 0 + if g.nvtxs < GPU_MIN_NVTXS[] && ME.RNG_TASK[] && Threads.nthreads() > 1 + c = _copy_ctrl(ctrl) + let g = g, lv = lv, c = c + push!(tasks, Threads.@spawn begin + Random.seed!(Random.TaskLocalRNG(), UInt32(hash((g.nvtxs, lv)) % UInt32)) + mlnd!(c, g, order, lv) + end) + end + elseif PARALLEL_TOP[] && Threads.nthreads() > 1 + c = _copy_ctrl(ctrl) + setup = TASK_SETUP[] + let g = g, lv = lv, c = c, setup = setup + push!(tasks, Threads.@spawn begin + setup === nothing || setup() + Random.seed!(Random.TaskLocalRNG(), UInt32(hash((g.nvtxs, lv)) % UInt32)) + dissect!(backend, c, g, order, lv) + end) + end + else + dissect!(backend, ctrl, g, order, lv) + end + else + ME.MMDOrder!(ctrl, g, order, lv) + end + end + foreach(wait, tasks) + return nothing +end + +""" + nodend_hybrid(backend, A; seed = -1) -> (perm, iperm) 1-based + +NodeND with the coarsening chain on the GpuND kernels and the rest on the +exact host port. Compression as in METIS_NodeND. +""" +function nodend_hybrid(backend, A; seed::Int = -1) + n = size(A, 1) + xadj = Vector{Int}(undef, n + 1) + adjncy = Int[] + rws = MetisExact.SparseArrays.rowvals(A) + xadj[1] = 0 + for j in 1:n + for p in MetisExact.SparseArrays.nzrange(A, j) + i = rws[p] + i != j && push!(adjncy, i - 1) + end + xadj[j + 1] = length(adjncy) + end + ctrl = ME.CtrlT(ME.METIS_CTYPE_SHEM, ME.METIS_RTYPE_SEP1SIDED, ME.METIS_IPTYPE_EDGE, + 1, 10, 200, true, false, 0.0, false, seed, 100, 0, + Float32(Float64(Float32(1.0 + 0.001 * 200)) + 0.0000499), Float32[]) + ME.RNG_TASK[] = Threads.nthreads() > 1 # threaded subtrees want task RNG + ME.InitRandom(ctrl.seed) + nnvtxs = 0 + graph = nothing + cptr = Int[] + cind = Int[] + if ctrl.compress + cptr = Vector{Int}(undef, n + 1) + cind = Vector{Int}(undef, n) + graph = ME.CompressGraph(n, xadj, adjncy, nothing, cptr, cind) + if graph === nothing + ctrl.compress = false + else + nnvtxs = graph.nvtxs + cfactor = 1.0 * n / nnvtxs + cfactor > 1.5 && ctrl.nseps == 1 && (ctrl.nseps = 2) + end + end + graph === nothing && (graph = ME.SetupGraph(n, xadj, adjncy, nothing)) + iperm = Vector{Int}(undef, n) + perm = Vector{Int}(undef, n) + dissect!(backend, ctrl, graph, iperm, graph.nvtxs) + if ctrl.compress + for i in 0:(nnvtxs - 1) + perm[iperm[i + 1] + 1] = i + end + l = 0 + for ii in 0:(nnvtxs - 1) + i = perm[ii + 1] + for j in cptr[i + 1]:(cptr[i + 2] - 1) + iperm[cind[j + 1] + 1] = l + l += 1 + end + end + end + for i in 0:(n - 1) + perm[iperm[i + 1] + 1] = i + end + return perm .+ 1, iperm .+ 1 +end + +end # module diff --git a/bench/ordering/MetisExact.jl b/bench/ordering/MetisExact.jl new file mode 100644 index 0000000..3d1b2b4 --- /dev/null +++ b/bench/ordering/MetisExact.jl @@ -0,0 +1,2411 @@ +module MetisExact + +using SparseArrays +using Libdl +using Random + +# Faithful Julia translation of METIS 5.1.0 (METIS_jll 5.1.4+0 = upstream +# 5.1.0 at amontoison/METIS e827ffe, build-system patches only), NodeND path +# with OMETIS defaults, idx_t = Int32 semantics, 0-based ids kept internally +# (C `a[i]` with value index i becomes `a[i + 1]`). Goal: bitwise-identical +# perm/iperm to libmetis before anything is moved to the GPU. +# +# Source of truth: ~/git/scratch/vixen-metis/metis-src/METIS/{libmetis,GKlib}. +# Every section names the C file and function it translates. + +export InitRandom, irand, irandInRange, irandArrayPermute! + +# --------------------------------------------------------------------------- +# GKlib/random.c — METIS_jll is built WITHOUT -DUSE_GKRAND (GKlib option +# default OFF, Yggdrasil does not set it), so the Mersenne-Twister code is +# compiled out and gk_randinit/gk_randint32 are libc srand()/rand(). We call +# the same libc, which is exact by construction on any host, and matches the +# shipped binary because both resolve to the process' C library. + +# Two RNG modes. :libc reproduces the shipped binary bitwise (process-global +# libc rand). RNG_TASK[] = true switches to Julia's task-local RNG so +# independent subtrees can run on threads deterministically (each task seeds +# its own stream); orderings are then METIS-quality but not bitwise-libmetis. +const RNG_TASK = Ref(false) + +function isrand(seed) + if RNG_TASK[] + Random.seed!(Random.TaskLocalRNG(), UInt32(seed % UInt32)) + else + ccall(:srand, Cvoid, (Cuint,), Cuint(seed % UInt32)) + end + return nothing +end +irand() = RNG_TASK[] ? Int(rand(Random.TaskLocalRNG(), UInt32) & 0x7fffffff) : + Int(ccall(:rand, Cint, ())) +irandInRange(max) = Int(mod(irand(), max)) + +# libmetis/util.c InitRandom +InitRandom(seed) = isrand(seed == -1 ? 4321 : seed) + +# gk_mkrandom.h randArrayPermute: p is 1-based storage of 0-based values when +# flag == 1; swaps are position swaps so the caller's base does not matter. +function irandArrayPermute!(n::Int, p::AbstractVector, nshuffles::Int, flag::Int) + if flag == 1 + for i in 1:n + p[i] = i - 1 + end + end + if n < 10 + for _ in 1:n + v = irandInRange(n) + 1 + u = irandInRange(n) + 1 + p[v], p[u] = p[u], p[v] + end + else + for _ in 1:nshuffles + v = irandInRange(n - 3) + 1 + u = irandInRange(n - 3) + 1 + # the template swaps (v+k, u+2), (v+1, u+3), (v+2, u+0), (v+3, u+1) + p[v], p[u + 2] = p[u + 2], p[v] + p[v + 1], p[u + 3] = p[u + 3], p[v + 1] + p[v + 2], p[u] = p[u], p[v + 2] + p[v + 3], p[u + 1] = p[u + 1], p[v + 3] + end + end + return p +end + + +# --------------------------------------------------------------------------- +# GKlib/gk_mksort.h GK_MKQSORT — the glibc-derived quicksort. Tie order is +# observable downstream (compression keys, matching orders), so the exact +# algorithm is translated, pointers becoming 1-based indices. + +const _QSORT_MAX_THRESH = 4 + +function gk_qsort!(v::AbstractVector, n::Int, lt::F) where {F} + n == 0 && return v + if n > _QSORT_MAX_THRESH + lo = 1 + hi = n + stack = Vector{Tuple{Int, Int}}() # (_lo, _hi); C pushes a sentinel + push!(stack, (0, 0)) + while true + mid = lo + ((hi - lo) >> 1) + if lt(v[mid], v[lo]) + v[mid], v[lo] = v[lo], v[mid] + end + jump = false + if lt(v[hi], v[mid]) + v[mid], v[hi] = v[hi], v[mid] + else + jump = true + end + if !jump && lt(v[mid], v[lo]) + v[mid], v[lo] = v[lo], v[mid] + end + left = lo + 1 + right = hi - 1 + while true + while lt(v[left], v[mid]) + left += 1 + end + while lt(v[mid], v[right]) + right -= 1 + end + if left < right + v[left], v[right] = v[right], v[left] + if mid == left + mid = right + elseif mid == right + mid = left + end + left += 1 + right -= 1 + elseif left == right + left += 1 + right -= 1 + break + end + left <= right || break + end + if right - lo <= _QSORT_MAX_THRESH + if hi - left <= _QSORT_MAX_THRESH + lo, hi = pop!(stack) + lo == 0 && break # popped the sentinel: done + else + lo = left + end + elseif hi - left <= _QSORT_MAX_THRESH + hi = right + elseif right - lo > hi - left + push!(stack, (lo, right)) + lo = left + else + push!(stack, (left, hi)) + hi = right + end + end + end + # insertion-sort tail + thresh = min(1 + _QSORT_MAX_THRESH, n) + tmp = 1 + for run in 2:thresh + lt(v[run], v[tmp]) && (tmp = run) + end + if tmp != 1 + v[tmp], v[1] = v[1], v[tmp] + end + run = 2 + while (run += 1) <= n + tmp = run - 1 + while lt(v[run], v[tmp]) + tmp -= 1 + end + tmp += 1 + if tmp != run + x = v[run] + for k in run:-1:(tmp + 1) + v[k] = v[k - 1] + end + v[tmp] = x + end + end + return v +end + +# libmetis/gklib.c — key/value pairs as tuples (key, val) +ikvsorti!(n::Int, v) = gk_qsort!(v, n, (a, b) -> a[1] < b[1]) +ikvsortii!(n::Int, v) = gk_qsort!(v, n, (a, b) -> a[1] < b[1] || (a[1] == b[1] && a[2] < b[2])) + +# libmetis/bucketsort.c BucketSortKeysInc (counting sort; keys in [0, max]) +function BucketSortKeysInc(n::Int, max::Int, keys, tperm, perm) + counts = zeros(Int, max + 2) + for i in 0:(n - 1) + counts[keys[i + 1] + 1] += 1 + end + for i in 1:max # MAKECSR over 0..max + counts[i + 1] += counts[i] + end + for i in (max + 1):-1:1 + counts[i + 1] = counts[i] + end + counts[1] = 0 + for ii in 0:(n - 1) + i = tperm[ii + 1] + perm[counts[keys[i + 1] + 1] + 1] = i + counts[keys[i + 1] + 1] += 1 + end + return perm +end + +# --------------------------------------------------------------------------- +# libmetis/struct.h graph_t (ncon = 1 fields only) — ids and offsets stay +# 0-based as in C; every array is a 1-based Julia Vector indexed [x + 1]. + +const NOVALUE = Int[] +const DBG = Ref(false) +const PROF = Dict{Symbol, Float64}() +const PROF_LOCK = ReentrantLock() +prof!(k, t) = (lock(() -> (PROF[k] = get(PROF, k, 0.0) + t), PROF_LOCK); nothing) + +mutable struct GraphT + nvtxs::Int + nedges::Int + xadj::Vector{Int} + vwgt::Vector{Int} + adjncy::Vector{Int} + adjwgt::Vector{Int} + tvwgt::Int + invtvwgt::Float32 # real_t = float32 in the shipped binary + label::Vector{Int} + cmap::Vector{Int} + # partition/refinement data + where_::Vector{Int} + pwgts::Vector{Int} + nbnd::Int + bndptr::Vector{Int} + bndind::Vector{Int} + id::Vector{Int} + ed::Vector{Int} + nrinfo::Vector{NTuple{2, Int}} # nrinfo_t.edegrees[2] + mincut::Int + coarser::Union{Nothing, GraphT} + finer::Union{Nothing, GraphT} +end + +CreateGraph() = GraphT(-1, -1, NOVALUE, NOVALUE, NOVALUE, NOVALUE, 0, 0.0f0, + NOVALUE, NOVALUE, NOVALUE, NOVALUE, -1, NOVALUE, NOVALUE, + NOVALUE, NOVALUE, NTuple{2, Int}[], -1, nothing, nothing) + +# libmetis/graph.c SetupGraph for OMETIS (ncon = 1, objtype = NODE) +function SetupGraph(nvtxs::Int, xadj::Vector{Int}, adjncy::Vector{Int}, + vwgt::Union{Nothing, Vector{Int}}) + graph = CreateGraph() + graph.nvtxs = nvtxs + graph.nedges = xadj[nvtxs + 1] + graph.xadj = xadj + graph.adjncy = adjncy + graph.vwgt = vwgt === nothing ? ones(Int, nvtxs) : vwgt + graph.adjwgt = ones(Int, graph.nedges) + SetupGraph_tvwgt!(graph) + SetupGraph_label!(graph) + return graph +end + +function SetupGraph_tvwgt!(graph::GraphT) + graph.tvwgt = sum(view(graph.vwgt, 1:graph.nvtxs); init = 0) + graph.invtvwgt = Float32(1.0 / (graph.tvwgt > 0 ? graph.tvwgt : 1)) + return nothing +end + +function SetupGraph_label!(graph::GraphT) + graph.label = collect(0:(graph.nvtxs - 1)) + return nothing +end + +function SetupSplitGraph(graph::GraphT, snvtxs::Int, snedges::Int) + sgraph = CreateGraph() + sgraph.nvtxs = snvtxs + sgraph.nedges = snedges + sgraph.xadj = Vector{Int}(undef, snvtxs + 1) + sgraph.vwgt = Vector{Int}(undef, snvtxs) + sgraph.adjncy = Vector{Int}(undef, snedges) + sgraph.adjwgt = Vector{Int}(undef, snedges) + sgraph.label = Vector{Int}(undef, snvtxs) + return sgraph +end + +# --------------------------------------------------------------------------- +# libmetis/compress.c CompressGraph — identical-adjacency vertices merge. +# Returns nothing when compression is below COMPRESSION_FRACTION. + +const COMPRESSION_FRACTION = 0.85 + +function CompressGraph(nvtxs::Int, xadj::Vector{Int}, adjncy::Vector{Int}, + vwgt::Union{Nothing, Vector{Int}}, cptr::Vector{Int}, cind::Vector{Int}) + mark = fill(-1, nvtxs) + map_ = fill(-1, nvtxs) + keys = Vector{Tuple{Int, Int}}(undef, nvtxs) + for i in 0:(nvtxs - 1) + k = 0 + for j in xadj[i + 1]:(xadj[i + 2] - 1) + k += adjncy[j + 1] + end + keys[i + 1] = (k + i, i) + end + ikvsorti!(nvtxs, keys) + l = 0 + cptr[1] = 0 + cnvtxs = 0 + for i in 0:(nvtxs - 1) + ii = keys[i + 1][2] + if map_[ii + 1] == -1 + mark[ii + 1] = i + for j in xadj[ii + 1]:(xadj[ii + 2] - 1) + mark[adjncy[j + 1] + 1] = i + end + map_[ii + 1] = cnvtxs + l += 1 + cind[l] = ii + for jj in (i + 1):(nvtxs - 1) + iii = keys[jj + 1][2] + if keys[i + 1][1] != keys[jj + 1][1] || + xadj[ii + 2] - xadj[ii + 1] != xadj[iii + 2] - xadj[iii + 1] + break + end + if map_[iii + 1] == -1 + ok = true + for p in xadj[iii + 1]:(xadj[iii + 2] - 1) + if mark[adjncy[p + 1] + 1] != i + ok = false + break + end + end + if ok + map_[iii + 1] = cnvtxs + l += 1 + cind[l] = iii + end + end + end + cnvtxs += 1 + cptr[cnvtxs + 1] = l + end + end + cnvtxs >= COMPRESSION_FRACTION * nvtxs && return nothing + graph = CreateGraph() + cnedges = 0 + for i in 0:(cnvtxs - 1) + ii = cind[cptr[i + 1] + 1] + cnedges += xadj[ii + 2] - xadj[ii + 1] + end + cxadj = Vector{Int}(undef, cnvtxs + 1) + cvwgt = zeros(Int, cnvtxs) + cadjncy = Vector{Int}(undef, cnedges) + fill!(mark, -1) + l = 0 + cxadj[1] = 0 + for i in 0:(cnvtxs - 1) + mark[i + 1] = i # remove diagonal in compressed graph + for j in cptr[i + 1]:(cptr[i + 2] - 1) + ii = cind[j + 1] + cvwgt[i + 1] += vwgt === nothing ? 1 : vwgt[ii + 1] + for jj in xadj[ii + 1]:(xadj[ii + 2] - 1) + k = map_[adjncy[jj + 1] + 1] + if mark[k + 1] != i + mark[k + 1] = i + l += 1 + cadjncy[l] = k + end + end + end + cxadj[i + 2] = l + end + graph.nvtxs = cnvtxs + graph.nedges = l + graph.xadj = cxadj + graph.vwgt = cvwgt + graph.adjncy = cadjncy + graph.adjwgt = ones(Int, l) + SetupGraph_tvwgt!(graph) + SetupGraph_label!(graph) + return graph +end + +# --------------------------------------------------------------------------- +# GKlib/gk_mkpqueue.h instantiated as rpq (libmetis/gklib.c: +# GK_MKPQUEUE(rpq, rpq_t, rkv_t, real_t, idx_t, ..., key_gt)): a MAX-heap +# with real_t = Float32 keys and idx_t vals. Heap layout (which of several +# equal-key items surfaces first) is observable through FM, so the sift +# logic is translated verbatim. KEY_LT(a, b) = a > b. + +mutable struct RPQ + nnodes::Int + heap_key::Vector{Float32} + heap_val::Vector{Int} + locator::Vector{Int} # by node id (0-based), -1 = absent +end + +rpqCreate(maxnodes::Int) = RPQ(0, Vector{Float32}(undef, maxnodes), + Vector{Int}(undef, maxnodes), fill(-1, maxnodes)) +rpqLength(q::RPQ) = q.nnodes + +@inline _rpq_lt(a::Float32, b::Float32) = a > b # KEY_LT = key_gt (max-heap) + +function rpqReset(q::RPQ) + for i in q.nnodes:-1:1 + q.locator[q.heap_val[i] + 1] = -1 + end + q.nnodes = 0 + return nothing +end + +function rpqInsert(q::RPQ, node::Int, key::Float32) + i = q.nnodes # 0-based heap position + q.nnodes += 1 + while i > 0 + j = (i - 1) >> 1 + if _rpq_lt(key, q.heap_key[j + 1]) + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + else + break + end + end + q.heap_key[i + 1] = key + q.heap_val[i + 1] = node + q.locator[node + 1] = i + return nothing +end + +function rpqDelete(q::RPQ, node::Int) + i = q.locator[node + 1] + q.locator[node + 1] = -1 + q.nnodes -= 1 + if q.nnodes > 0 && q.heap_val[q.nnodes + 1] != node + node = q.heap_val[q.nnodes + 1] + newkey = q.heap_key[q.nnodes + 1] + oldkey = q.heap_key[i + 1] + if _rpq_lt(newkey, oldkey) # filter up + while i > 0 + j = (i - 1) >> 1 + if _rpq_lt(newkey, q.heap_key[j + 1]) + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + else + break + end + end + else # filter down + nnodes = q.nnodes + while (j = (i << 1) + 1) < nnodes + if _rpq_lt(q.heap_key[j + 1], newkey) + if j + 1 < nnodes && _rpq_lt(q.heap_key[j + 2], q.heap_key[j + 1]) + j += 1 + end + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + elseif j + 1 < nnodes && _rpq_lt(q.heap_key[j + 2], newkey) + j += 1 + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + else + break + end + end + end + q.heap_key[i + 1] = newkey + q.heap_val[i + 1] = node + q.locator[node + 1] = i + end + return nothing +end + +function rpqUpdate(q::RPQ, node::Int, newkey::Float32) + i = q.locator[node + 1] + oldkey = q.heap_key[i + 1] + if _rpq_lt(newkey, oldkey) # filter up + while i > 0 + j = (i - 1) >> 1 + if _rpq_lt(newkey, q.heap_key[j + 1]) + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + else + break + end + end + else # filter down + nnodes = q.nnodes + while (j = (i << 1) + 1) < nnodes + if _rpq_lt(q.heap_key[j + 1], newkey) + if j + 1 < nnodes && _rpq_lt(q.heap_key[j + 2], q.heap_key[j + 1]) + j += 1 + end + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + elseif j + 1 < nnodes && _rpq_lt(q.heap_key[j + 2], newkey) + j += 1 + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + else + break + end + end + end + q.heap_key[i + 1] = newkey + q.heap_val[i + 1] = node + q.locator[node + 1] = i + return nothing +end + +function rpqGetTop(q::RPQ) + q.nnodes == 0 && return -1 + q.nnodes -= 1 + vtx = q.heap_val[1] + q.locator[vtx + 1] = -1 + i = q.nnodes + if i > 0 + key = q.heap_key[i + 1] + node = q.heap_val[i + 1] + i = 0 + while (j = 2 * i + 1) < q.nnodes + if _rpq_lt(q.heap_key[j + 1], key) + if j + 1 < q.nnodes && _rpq_lt(q.heap_key[j + 2], q.heap_key[j + 1]) + j += 1 + end + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + elseif j + 1 < q.nnodes && _rpq_lt(q.heap_key[j + 2], key) + j += 1 + q.heap_key[i + 1] = q.heap_key[j + 1] + q.heap_val[i + 1] = q.heap_val[j + 1] + q.locator[q.heap_val[i + 1] + 1] = i + i = j + else + break + end + end + q.heap_key[i + 1] = key + q.heap_val[i + 1] = node + q.locator[node + 1] = i + end + return vtx +end + + +# --------------------------------------------------------------------------- +# ctrl_t subset for METIS_OP_OMETIS (libmetis/options.c SetupCtrl). All +# values that the OMETIS path reads. idx arithmetic stays in Int64: every C +# idx_t intermediate on this path fits Int32 (checked per use), so Int64 +# arithmetic is value-identical. + +mutable struct CtrlT + ctype::Int # METIS_CTYPE_SHEM = 2 + rtype::Int # METIS_RTYPE_SEP1SIDED = 3 / SEP2SIDED = 4 + iptype::Int # METIS_IPTYPE_EDGE = 2 + nseps::Int + niter::Int + ufactor::Int + compress::Bool + ccorder::Bool + pfactor::Float64 + no2hop::Bool + seed::Int + CoarsenTo::Int + maxvwgt::Int + ubfactors::Float32 # real_t: I2RUBFACTOR(ufactor) + 0.0000499 + pijbm::Vector{Float32} # real_t: Setup2WayBalMultipliers (len 2) +end + +const METIS_CTYPE_RM = 0 +const METIS_CTYPE_SHEM = 1 +const METIS_RTYPE_FM = 0 +const METIS_RTYPE_SEP2SIDED = 3 +const METIS_RTYPE_SEP1SIDED = 2 +const METIS_IPTYPE_EDGE = 2 +const METIS_IPTYPE_NODE = 3 +const UNMATCHED = -1 +const COARSEN_FRACTION = 0.85 +const UNMATCHEDFOR2HOP = 0.10 +const HTLENGTH = (1 << 11) - 1 +const MMDSWITCH = 120 +const LARGENIPARTS = 7 +const SMALLNIPARTS = 5 +const IDX_MAX = 2147483647 + +# --------------------------------------------------------------------------- +# libmetis/coarsen.c + +function CoarsenGraph(ctrl::CtrlT, graph::GraphT) + eqewgts = true + for i in 1:(graph.nedges - 1) + if graph.adjwgt[1] != graph.adjwgt[i + 1] + eqewgts = false + break + end + end + ctrl.maxvwgt = trunc(Int, 1.5 * graph.tvwgt / ctrl.CoarsenTo) + while true + DBG[] && println(lpad(graph.nvtxs, 10), lpad(graph.nedges, 11), + lpad(sum(view(graph.adjwgt, 1:graph.nedges); init = 0), 11), + " [", ctrl.CoarsenTo, "] [", lpad(ctrl.maxvwgt, 9), ":", lpad(graph.tvwgt, 8), " ]") + isempty(graph.cmap) && (graph.cmap = Vector{Int}(undef, graph.nvtxs)) + if ctrl.ctype == METIS_CTYPE_RM || eqewgts || graph.nedges == 0 + Match_RM(ctrl, graph) + else + Match_SHEM(ctrl, graph) + end + graph = graph.coarser::GraphT + eqewgts = false + (graph.nvtxs > ctrl.CoarsenTo && + graph.nvtxs < COARSEN_FRACTION * (graph.finer::GraphT).nvtxs && + graph.nedges > graph.nvtxs ÷ 2) || break + end + DBG[] && println(lpad(graph.nvtxs, 10), lpad(graph.nedges, 11), + lpad(sum(view(graph.adjwgt, 1:graph.nedges); init = 0), 11), + " [", ctrl.CoarsenTo, "] [", lpad(ctrl.maxvwgt, 9), ":", lpad(graph.tvwgt, 8), " ]") + return graph +end + +function CoarsenGraphNlevels(ctrl::CtrlT, graph::GraphT, nlevels::Int) + eqewgts = true + for i in 1:(graph.nedges - 1) + if graph.adjwgt[1] != graph.adjwgt[i + 1] + eqewgts = false + break + end + end + ctrl.maxvwgt = trunc(Int, 1.5 * graph.tvwgt / ctrl.CoarsenTo) + for _ in 1:nlevels + isempty(graph.cmap) && (graph.cmap = Vector{Int}(undef, graph.nvtxs)) + if ctrl.ctype == METIS_CTYPE_RM || eqewgts || graph.nedges == 0 + Match_RM(ctrl, graph) + else + Match_SHEM(ctrl, graph) + end + graph = graph.coarser::GraphT + eqewgts = false + if graph.nvtxs < ctrl.CoarsenTo || + graph.nvtxs > COARSEN_FRACTION * (graph.finer::GraphT).nvtxs || + graph.nedges < graph.nvtxs ÷ 2 + break + end + end + return graph +end + +function Match_RM(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + cmap = graph.cmap + maxvwgt = ctrl.maxvwgt + nunmatched = 0 + match = fill(UNMATCHED, nvtxs) + perm = Vector{Int}(undef, nvtxs) + irandArrayPermute!(nvtxs, perm, nvtxs ÷ 8, 1) + cnvtxs = 0 + last_unmatched = 0 # 0-based position + for pi in 0:(nvtxs - 1) + i = perm[pi + 1] + match[i + 1] == UNMATCHED || continue + maxidx = i + if vwgt[i + 1] < maxvwgt + if xadj[i + 1] == xadj[i + 2] # island vertex + last_unmatched = max(pi, last_unmatched) + 1 + while last_unmatched < nvtxs + j = perm[last_unmatched + 1] + if match[j + 1] == UNMATCHED + maxidx = j + break + end + last_unmatched += 1 + end + else + for j in xadj[i + 1]:(xadj[i + 2] - 1) + k = adjncy[j + 1] + if match[k + 1] == UNMATCHED && vwgt[i + 1] + vwgt[k + 1] <= maxvwgt + maxidx = k + break + end + end + if maxidx == i && 3 * vwgt[i + 1] < maxvwgt + nunmatched += 1 + maxidx = UNMATCHED + end + end + end + if maxidx != UNMATCHED + cmap[i + 1] = cmap[maxidx + 1] = cnvtxs + cnvtxs += 1 + match[i + 1] = maxidx + match[maxidx + 1] = i + end + end + if !ctrl.no2hop && nunmatched > UNMATCHEDFOR2HOP * nvtxs + cnvtxs = Match_2Hop(ctrl, graph, perm, match, cnvtxs, nunmatched) + end + cnvtxs = 0 + for i in 0:(nvtxs - 1) + if match[i + 1] == UNMATCHED + match[i + 1] = i + cmap[i + 1] = cnvtxs + cnvtxs += 1 + elseif i <= match[i + 1] + cmap[i + 1] = cmap[match[i + 1] + 1] = cnvtxs + cnvtxs += 1 + end + end + CreateCoarseGraph(ctrl, graph, cnvtxs, match) + return cnvtxs +end + +function Match_SHEM(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + cmap = graph.cmap + maxvwgt = ctrl.maxvwgt + nunmatched = 0 + match = fill(UNMATCHED, nvtxs) + perm = Vector{Int}(undef, nvtxs) + tperm = Vector{Int}(undef, nvtxs) + degrees = Vector{Int}(undef, nvtxs) + irandArrayPermute!(nvtxs, tperm, nvtxs ÷ 8, 1) + avgdegree = trunc(Int, 0.7 * (xadj[nvtxs + 1] ÷ nvtxs)) + for i in 0:(nvtxs - 1) + d = xadj[i + 2] - xadj[i + 1] + degrees[i + 1] = d > avgdegree ? avgdegree : d + end + BucketSortKeysInc(nvtxs, avgdegree, degrees, tperm, perm) + cnvtxs = 0 + last_unmatched = 0 + for pi in 0:(nvtxs - 1) + i = perm[pi + 1] + match[i + 1] == UNMATCHED || continue + maxidx = i + maxwgt = -1 + if vwgt[i + 1] < maxvwgt + if xadj[i + 1] == xadj[i + 2] # island vertex + last_unmatched = max(pi, last_unmatched) + 1 + while last_unmatched < nvtxs + j = perm[last_unmatched + 1] + if match[j + 1] == UNMATCHED + maxidx = j + break + end + last_unmatched += 1 + end + else + for j in xadj[i + 1]:(xadj[i + 2] - 1) + k = adjncy[j + 1] + if match[k + 1] == UNMATCHED && maxwgt < adjwgt[j + 1] && + vwgt[i + 1] + vwgt[k + 1] <= maxvwgt + maxidx = k + maxwgt = adjwgt[j + 1] + end + end + if maxidx == i && 3 * vwgt[i + 1] < maxvwgt + nunmatched += 1 + maxidx = UNMATCHED + end + end + end + if maxidx != UNMATCHED + cmap[i + 1] = cmap[maxidx + 1] = cnvtxs + cnvtxs += 1 + match[i + 1] = maxidx + match[maxidx + 1] = i + end + end + if !ctrl.no2hop && nunmatched > UNMATCHEDFOR2HOP * nvtxs + cnvtxs = Match_2Hop(ctrl, graph, perm, match, cnvtxs, nunmatched) + end + cnvtxs = 0 + for i in 0:(nvtxs - 1) + if match[i + 1] == UNMATCHED + match[i + 1] = i + cmap[i + 1] = cnvtxs + cnvtxs += 1 + elseif i <= match[i + 1] + cmap[i + 1] = cmap[match[i + 1] + 1] = cnvtxs + cnvtxs += 1 + end + end + CreateCoarseGraph(ctrl, graph, cnvtxs, match) + return cnvtxs +end + +function Match_2Hop(ctrl::CtrlT, graph::GraphT, perm, match, cnvtxs::Int, nunmatched::Int) + cnvtxs, nunmatched = Match_2HopAny(ctrl, graph, perm, match, cnvtxs, nunmatched, 2) + cnvtxs, nunmatched = Match_2HopAll(ctrl, graph, perm, match, cnvtxs, nunmatched, 64) + if nunmatched > 1.5 * UNMATCHEDFOR2HOP * graph.nvtxs + cnvtxs, nunmatched = Match_2HopAny(ctrl, graph, perm, match, cnvtxs, nunmatched, 3) + end + if nunmatched > 2.0 * UNMATCHEDFOR2HOP * graph.nvtxs + cnvtxs, nunmatched = Match_2HopAny(ctrl, graph, perm, match, cnvtxs, nunmatched, graph.nvtxs) + end + return cnvtxs +end + +function Match_2HopAny(ctrl::CtrlT, graph::GraphT, perm, match, cnvtxs::Int, + nunmatched::Int, maxdegree::Int) + nvtxs = graph.nvtxs + xadj = graph.xadj + adjncy = graph.adjncy + cmap = graph.cmap + colptr = zeros(Int, nvtxs + 1) + for i in 0:(nvtxs - 1) + if match[i + 1] == UNMATCHED && xadj[i + 2] - xadj[i + 1] < maxdegree + for j in xadj[i + 1]:(xadj[i + 2] - 1) + colptr[adjncy[j + 1] + 1] += 1 + end + end + end + for i in 1:(nvtxs - 1) # MAKECSR + colptr[i + 1] += colptr[i] + end + for i in nvtxs:-1:1 + colptr[i + 1] = colptr[i] + end + colptr[1] = 0 + rowind = Vector{Int}(undef, colptr[nvtxs + 1]) + for pi in 0:(nvtxs - 1) + i = perm[pi + 1] + if match[i + 1] == UNMATCHED && xadj[i + 2] - xadj[i + 1] < maxdegree + for j in xadj[i + 1]:(xadj[i + 2] - 1) + c = adjncy[j + 1] + rowind[colptr[c + 1] + 1] = i + colptr[c + 1] += 1 + end + end + end + for i in nvtxs:-1:1 # SHIFTCSR + colptr[i + 1] = colptr[i] + end + colptr[1] = 0 + for pi in 0:(nvtxs - 1) + i = perm[pi + 1] + colptr[i + 2] - colptr[i + 1] < 2 && continue + jj = colptr[i + 2] + j = colptr[i + 1] + while j < jj + if match[rowind[j + 1] + 1] == UNMATCHED + jj -= 1 + while jj > j + if match[rowind[jj + 1] + 1] == UNMATCHED + cmap[rowind[j + 1] + 1] = cmap[rowind[jj + 1] + 1] = cnvtxs + cnvtxs += 1 + match[rowind[j + 1] + 1] = rowind[jj + 1] + match[rowind[jj + 1] + 1] = rowind[j + 1] + nunmatched -= 2 + break + end + jj -= 1 + end + end + j += 1 + end + end + return cnvtxs, nunmatched +end + +function Match_2HopAll(ctrl::CtrlT, graph::GraphT, perm, match, cnvtxs::Int, + nunmatched::Int, maxdegree::Int) + nvtxs = graph.nvtxs + xadj = graph.xadj + adjncy = graph.adjncy + cmap = graph.cmap + mask = IDX_MAX ÷ maxdegree + keys = Vector{Tuple{Int, Int}}(undef, nunmatched) + ncand = 0 + for pi in 0:(nvtxs - 1) + i = perm[pi + 1] + idegree = xadj[i + 2] - xadj[i + 1] + if match[i + 1] == UNMATCHED && idegree > 1 && idegree < maxdegree + k = 0 + for j in xadj[i + 1]:(xadj[i + 2] - 1) + k += adjncy[j + 1] % mask + end + ncand += 1 + keys[ncand] = ((k % mask) * maxdegree + idegree, i) + end + end + ikvsorti!(ncand, keys) + mark = zeros(Int, nvtxs) + for pi in 0:(ncand - 1) + i = keys[pi + 1][2] + match[i + 1] == UNMATCHED || continue + for j in xadj[i + 1]:(xadj[i + 2] - 1) + mark[adjncy[j + 1] + 1] = i + end + for pk in (pi + 1):(ncand - 1) + k = keys[pk + 1][2] + match[k + 1] == UNMATCHED || continue + keys[pi + 1][1] != keys[pk + 1][1] && break + xadj[i + 2] - xadj[i + 1] != xadj[k + 2] - xadj[k + 1] && break + allmatch = true + for jj in xadj[k + 1]:(xadj[k + 2] - 1) + if mark[adjncy[jj + 1] + 1] != i + allmatch = false + break + end + end + if allmatch + cmap[i + 1] = cmap[k + 1] = cnvtxs + cnvtxs += 1 + match[i + 1] = k + match[k + 1] = i + nunmatched -= 2 + break + end + end + end + return cnvtxs, nunmatched +end + +function SetupCoarseGraph(graph::GraphT, cnvtxs::Int) + cgraph = CreateGraph() + cgraph.nvtxs = cnvtxs + cgraph.finer = graph + graph.coarser = cgraph + cgraph.xadj = Vector{Int}(undef, cnvtxs + 1) + cgraph.adjncy = Vector{Int}(undef, graph.nedges) + cgraph.adjwgt = Vector{Int}(undef, graph.nedges) + cgraph.vwgt = Vector{Int}(undef, cnvtxs) + return cgraph +end + +function CreateCoarseGraph(ctrl::CtrlT, graph::GraphT, cnvtxs::Int, match) + mask = HTLENGTH + if cnvtxs < 2 * mask || graph.nedges ÷ graph.nvtxs > mask ÷ 20 + return CreateCoarseGraphNoMask(ctrl, graph, cnvtxs, match) + end + nvtxs = graph.nvtxs + xadj = graph.xadj + for v in 0:(nvtxs - 1) + if xadj[v + 2] - xadj[v + 1] > (mask >> 3) + return CreateCoarseGraphNoMask(ctrl, graph, cnvtxs, match) + end + end + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + cmap = graph.cmap + cgraph = SetupCoarseGraph(graph, cnvtxs) + cxadj = cgraph.xadj + cvwgt = cgraph.vwgt + cadjncy = cgraph.adjncy + cadjwgt = cgraph.adjwgt + htable = fill(-1, mask + 1) + cxadj[1] = 0 + cnvtxs = 0 + cnedges = 0 # base of the current coarse row + for v in 0:(nvtxs - 1) + u = match[v + 1] + u < v && continue + cvwgt[cnvtxs + 1] = vwgt[v + 1] + nedges = 0 + for j in xadj[v + 1]:(xadj[v + 2] - 1) + k = cmap[adjncy[j + 1] + 1] + kk = k & mask + m = htable[kk + 1] + if m == -1 + cadjncy[cnedges + nedges + 1] = k + cadjwgt[cnedges + nedges + 1] = adjwgt[j + 1] + htable[kk + 1] = nedges + nedges += 1 + elseif cadjncy[cnedges + m + 1] == k + cadjwgt[cnedges + m + 1] += adjwgt[j + 1] + else + jj = 0 + while jj < nedges + if cadjncy[cnedges + jj + 1] == k + cadjwgt[cnedges + jj + 1] += adjwgt[j + 1] + break + end + jj += 1 + end + if jj == nedges + cadjncy[cnedges + nedges + 1] = k + cadjwgt[cnedges + nedges + 1] = adjwgt[j + 1] + nedges += 1 + end + end + end + if v != u + cvwgt[cnvtxs + 1] += vwgt[u + 1] + for j in xadj[u + 1]:(xadj[u + 2] - 1) + k = cmap[adjncy[j + 1] + 1] + kk = k & mask + m = htable[kk + 1] + if m == -1 + cadjncy[cnedges + nedges + 1] = k + cadjwgt[cnedges + nedges + 1] = adjwgt[j + 1] + htable[kk + 1] = nedges + nedges += 1 + elseif cadjncy[cnedges + m + 1] == k + cadjwgt[cnedges + m + 1] += adjwgt[j + 1] + else + jj = 0 + while jj < nedges + if cadjncy[cnedges + jj + 1] == k + cadjwgt[cnedges + jj + 1] += adjwgt[j + 1] + break + end + jj += 1 + end + if jj == nedges + cadjncy[cnedges + nedges + 1] = k + cadjwgt[cnedges + nedges + 1] = adjwgt[j + 1] + nedges += 1 + end + end + end + # remove the contracted adjacency weight + jj = htable[(cnvtxs & mask) + 1] + if jj >= 0 && cadjncy[cnedges + jj + 1] != cnvtxs + jj = 0 + while jj < nedges + cadjncy[cnedges + jj + 1] == cnvtxs && break + jj += 1 + end + end + if jj >= 0 && jj < nedges && cadjncy[cnedges + jj + 1] == cnvtxs + nedges -= 1 + cadjncy[cnedges + jj + 1] = cadjncy[cnedges + nedges + 1] + cadjwgt[cnedges + jj + 1] = cadjwgt[cnedges + nedges + 1] + end + end + for j in 0:(nedges - 1) + htable[(cadjncy[cnedges + j + 1] & mask) + 1] = -1 + end + htable[(cnvtxs & mask) + 1] = -1 + cnedges += nedges + cnvtxs += 1 + cxadj[cnvtxs + 1] = cnedges + end + cgraph.nedges = cnedges + SetupGraph_tvwgt!(cgraph) + return nothing +end + +function CreateCoarseGraphNoMask(ctrl::CtrlT, graph::GraphT, cnvtxs::Int, match) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + cmap = graph.cmap + cgraph = SetupCoarseGraph(graph, cnvtxs) + cxadj = cgraph.xadj + cvwgt = cgraph.vwgt + cadjncy = cgraph.adjncy + cadjwgt = cgraph.adjwgt + htable = fill(-1, cnvtxs) + cxadj[1] = 0 + cnvtxs = 0 + cnedges = 0 + for v in 0:(nvtxs - 1) + u = match[v + 1] + u < v && continue + cvwgt[cnvtxs + 1] = vwgt[v + 1] + nedges = 0 + for j in xadj[v + 1]:(xadj[v + 2] - 1) + k = cmap[adjncy[j + 1] + 1] + m = htable[k + 1] + if m == -1 + cadjncy[cnedges + nedges + 1] = k + cadjwgt[cnedges + nedges + 1] = adjwgt[j + 1] + htable[k + 1] = nedges + nedges += 1 + else + cadjwgt[cnedges + m + 1] += adjwgt[j + 1] + end + end + if v != u + cvwgt[cnvtxs + 1] += vwgt[u + 1] + for j in xadj[u + 1]:(xadj[u + 2] - 1) + k = cmap[adjncy[j + 1] + 1] + m = htable[k + 1] + if m == -1 + cadjncy[cnedges + nedges + 1] = k + cadjwgt[cnedges + nedges + 1] = adjwgt[j + 1] + htable[k + 1] = nedges + nedges += 1 + else + cadjwgt[cnedges + m + 1] += adjwgt[j + 1] + end + end + j = htable[cnvtxs + 1] + if j != -1 + nedges -= 1 + cadjncy[cnedges + j + 1] = cadjncy[cnedges + nedges + 1] + cadjwgt[cnedges + j + 1] = cadjwgt[cnedges + nedges + 1] + htable[cnvtxs + 1] = -1 + end + end + for j in 0:(nedges - 1) + htable[cadjncy[cnedges + j + 1] + 1] = -1 + end + cnedges += nedges + cnvtxs += 1 + cxadj[cnvtxs + 1] = cnedges + end + cgraph.nedges = cnedges + SetupGraph_tvwgt!(cgraph) + return nothing +end + + +rpqSeeTopVal(q::RPQ) = q.nnodes == 0 ? -1 : q.heap_val[1] + +# --------------------------------------------------------------------------- +# libmetis/refine.c (2-way edge-cut params; used at the coarsest level only) + +function Allocate2WayPartitionMemory!(graph::GraphT) + graph.pwgts = Vector{Int}(undef, 2) + graph.where_ = Vector{Int}(undef, graph.nvtxs) + graph.bndptr = Vector{Int}(undef, graph.nvtxs) + graph.bndind = Vector{Int}(undef, graph.nvtxs) + graph.id = Vector{Int}(undef, graph.nvtxs) + graph.ed = Vector{Int}(undef, graph.nvtxs) + return nothing +end + +function Compute2WayPartitionParams!(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + where_ = graph.where_ + id = graph.id + ed = graph.ed + pwgts = graph.pwgts + fill!(pwgts, 0) + bndptr = graph.bndptr + fill!(bndptr, -1) + bndind = graph.bndind + for i in 0:(nvtxs - 1) + pwgts[where_[i + 1] + 1] += vwgt[i + 1] + end + nbnd = 0 + mincut = 0 + for i in 0:(nvtxs - 1) + istart = xadj[i + 1] + iend = xadj[i + 2] + me = where_[i + 1] + tid = 0 + ted = 0 + for j in istart:(iend - 1) + if me == where_[adjncy[j + 1] + 1] + tid += adjwgt[j + 1] + else + ted += adjwgt[j + 1] + end + end + id[i + 1] = tid + ed[i + 1] = ted + if ted > 0 || istart == iend + bndind[nbnd + 1] = i # BNDInsert + bndptr[i + 1] = nbnd + nbnd += 1 + mincut += ted + end + end + graph.mincut = mincut ÷ 2 + graph.nbnd = nbnd + return nothing +end + +# --------------------------------------------------------------------------- +# libmetis/fm.c FM_2WayRefine / FM_2WayCutRefine (ncon = 1) + +FM_2WayRefine!(ctrl::CtrlT, graph::GraphT, ntpwgts, niter::Int) = + FM_2WayCutRefine!(ctrl, graph, ntpwgts, niter) + +function FM_2WayCutRefine!(ctrl::CtrlT, graph::GraphT, ntpwgts, niter::Int) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + where_ = graph.where_ + id = graph.id + ed = graph.ed + pwgts = graph.pwgts + bndptr = graph.bndptr + bndind = graph.bndind + moved = fill(-1, nvtxs) + swaps = Vector{Int}(undef, nvtxs) + perm = Vector{Int}(undef, nvtxs) + tpwgts0 = trunc(Int, graph.tvwgt * ntpwgts[1]) + tpwgts = (tpwgts0, graph.tvwgt - tpwgts0) + limit = trunc(Int, min(max(0.01 * nvtxs, 15.0), 100.0)) + avgvwgt = min((pwgts[1] + pwgts[2]) ÷ 20, 2 * (pwgts[1] + pwgts[2]) ÷ nvtxs) + queues = (rpqCreate(nvtxs), rpqCreate(nvtxs)) + origdiff = abs(tpwgts[1] - pwgts[1]) + for _pass in 0:(niter - 1) + rpqReset(queues[1]) + rpqReset(queues[2]) + mincutorder = -1 + newcut = mincut = initcut = graph.mincut + mindiff = abs(tpwgts[1] - pwgts[1]) + nbnd = graph.nbnd + irandArrayPermute!(nbnd, perm, nbnd, 1) + for ii in 0:(nbnd - 1) + i = perm[ii + 1] + b = bndind[i + 1] + rpqInsert(queues[where_[b + 1] + 1], b, Float32(ed[b + 1] - id[b + 1])) + end + nswaps = 0 + while nswaps < nvtxs + from = tpwgts[1] - pwgts[1] < tpwgts[2] - pwgts[2] ? 0 : 1 + to = (from + 1) % 2 + higain = rpqGetTop(queues[from + 1]) + higain == -1 && break + newcut -= ed[higain + 1] - id[higain + 1] + pwgts[to + 1] += vwgt[higain + 1] + pwgts[from + 1] -= vwgt[higain + 1] + if (newcut < mincut && abs(tpwgts[1] - pwgts[1]) <= origdiff + avgvwgt) || + (newcut == mincut && abs(tpwgts[1] - pwgts[1]) < mindiff) + mincut = newcut + mindiff = abs(tpwgts[1] - pwgts[1]) + mincutorder = nswaps + elseif nswaps - mincutorder > limit + newcut += ed[higain + 1] - id[higain + 1] + pwgts[from + 1] += vwgt[higain + 1] + pwgts[to + 1] -= vwgt[higain + 1] + break + end + where_[higain + 1] = to + moved[higain + 1] = nswaps + swaps[nswaps + 1] = higain + id[higain + 1], ed[higain + 1] = ed[higain + 1], id[higain + 1] + if ed[higain + 1] == 0 && xadj[higain + 1] < xadj[higain + 2] + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] # BNDDelete + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + end + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + kwgt = to == where_[k + 1] ? adjwgt[j + 1] : -adjwgt[j + 1] + id[k + 1] += kwgt + ed[k + 1] -= kwgt + if bndptr[k + 1] != -1 + if ed[k + 1] == 0 + bndind[bndptr[k + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[k + 1] + nbnd -= 1 + bndptr[k + 1] = -1 + moved[k + 1] == -1 && rpqDelete(queues[where_[k + 1] + 1], k) + elseif moved[k + 1] == -1 + rpqUpdate(queues[where_[k + 1] + 1], k, Float32(ed[k + 1] - id[k + 1])) + end + elseif ed[k + 1] > 0 + bndind[nbnd + 1] = k + bndptr[k + 1] = nbnd + nbnd += 1 + moved[k + 1] == -1 && rpqInsert(queues[where_[k + 1] + 1], k, Float32(ed[k + 1] - id[k + 1])) + end + end + nswaps += 1 + end + for i in 0:(nswaps - 1) + moved[swaps[i + 1] + 1] = -1 + end + nswaps -= 1 + while nswaps > mincutorder + higain = swaps[nswaps + 1] + to = where_[higain + 1] = (where_[higain + 1] + 1) % 2 + id[higain + 1], ed[higain + 1] = ed[higain + 1], id[higain + 1] + if ed[higain + 1] == 0 && bndptr[higain + 1] != -1 && xadj[higain + 1] < xadj[higain + 2] + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + elseif ed[higain + 1] > 0 && bndptr[higain + 1] == -1 + bndind[nbnd + 1] = higain + bndptr[higain + 1] = nbnd + nbnd += 1 + end + pwgts[to + 1] += vwgt[higain + 1] + pwgts[((to + 1) % 2) + 1] -= vwgt[higain + 1] + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + kwgt = to == where_[k + 1] ? adjwgt[j + 1] : -adjwgt[j + 1] + id[k + 1] += kwgt + ed[k + 1] -= kwgt + if bndptr[k + 1] != -1 && ed[k + 1] == 0 + bndind[bndptr[k + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[k + 1] + nbnd -= 1 + bndptr[k + 1] = -1 + end + if bndptr[k + 1] == -1 && ed[k + 1] > 0 + bndind[nbnd + 1] = k + bndptr[k + 1] = nbnd + nbnd += 1 + end + end + nswaps -= 1 + end + graph.mincut = mincut + graph.nbnd = nbnd + (mincutorder <= 0 || mincut == initcut) && break + end + return nothing +end + +# --------------------------------------------------------------------------- +# libmetis/balance.c (ncon = 1) + +# mcutil.c ComputeLoadImbalanceDiff (ncon = 1, nparts = 2), real_t = Float32 +function ComputeLoadImbalanceDiff(graph::GraphT, pijbm::Vector{Float32}, ub::Float32) + mx = -1.0f0 + for j in 1:2 + cur = Float32(graph.pwgts[j]) * pijbm[j] - ub + cur > mx && (mx = cur) + end + return mx +end + +function Balance2Way!(ctrl::CtrlT, graph::GraphT, ntpwgts) + ComputeLoadImbalanceDiff(graph, ctrl.pijbm, ctrl.ubfactors) <= 0 && return nothing + abs(trunc(Int, ntpwgts[1] * Float32(graph.tvwgt) - graph.pwgts[1])) < + 3 * graph.tvwgt ÷ graph.nvtxs && return nothing + if graph.nbnd > 0 + Bnd2WayBalance!(ctrl, graph, ntpwgts) + else + General2WayBalance!(ctrl, graph, ntpwgts) + end + return nothing +end + +function Bnd2WayBalance!(ctrl::CtrlT, graph::GraphT, ntpwgts) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + where_ = graph.where_ + id = graph.id + ed = graph.ed + pwgts = graph.pwgts + bndptr = graph.bndptr + bndind = graph.bndind + moved = fill(-1, nvtxs) + perm = Vector{Int}(undef, nvtxs) + tpwgts0 = trunc(Int, graph.tvwgt * ntpwgts[1]) + tpwgts = (tpwgts0, graph.tvwgt - tpwgts0) + mindiff = abs(tpwgts[1] - pwgts[1]) + from = pwgts[1] < tpwgts[1] ? 1 : 0 + to = (from + 1) % 2 + queue = rpqCreate(nvtxs) + nbnd = graph.nbnd + irandArrayPermute!(nbnd, perm, nbnd ÷ 5, 1) + for ii in 0:(nbnd - 1) + i = perm[ii + 1] + b = bndind[i + 1] + if where_[b + 1] == from && vwgt[b + 1] <= mindiff + rpqInsert(queue, b, Float32(ed[b + 1] - id[b + 1])) + end + end + mincut = graph.mincut + for _ in 0:(nvtxs - 1) + higain = rpqGetTop(queue) + higain == -1 && break + pwgts[to + 1] + vwgt[higain + 1] > tpwgts[to + 1] && break + mincut -= ed[higain + 1] - id[higain + 1] + pwgts[to + 1] += vwgt[higain + 1] + pwgts[from + 1] -= vwgt[higain + 1] + where_[higain + 1] = to + moved[higain + 1] = 1 + id[higain + 1], ed[higain + 1] = ed[higain + 1], id[higain + 1] + if ed[higain + 1] == 0 && xadj[higain + 1] < xadj[higain + 2] + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + end + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + kwgt = to == where_[k + 1] ? adjwgt[j + 1] : -adjwgt[j + 1] + id[k + 1] += kwgt + ed[k + 1] -= kwgt + if bndptr[k + 1] != -1 + if ed[k + 1] == 0 + bndind[bndptr[k + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[k + 1] + nbnd -= 1 + bndptr[k + 1] = -1 + if moved[k + 1] == -1 && where_[k + 1] == from && vwgt[k + 1] <= mindiff + rpqDelete(queue, k) + end + elseif moved[k + 1] == -1 && where_[k + 1] == from && vwgt[k + 1] <= mindiff + rpqUpdate(queue, k, Float32(ed[k + 1] - id[k + 1])) + end + elseif ed[k + 1] > 0 + bndind[nbnd + 1] = k + bndptr[k + 1] = nbnd + nbnd += 1 + if moved[k + 1] == -1 && where_[k + 1] == from && vwgt[k + 1] <= mindiff + rpqInsert(queue, k, Float32(ed[k + 1] - id[k + 1])) + end + end + end + end + graph.mincut = mincut + graph.nbnd = nbnd + return nothing +end + +function General2WayBalance!(ctrl::CtrlT, graph::GraphT, ntpwgts) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + adjwgt = graph.adjwgt + where_ = graph.where_ + id = graph.id + ed = graph.ed + pwgts = graph.pwgts + bndptr = graph.bndptr + bndind = graph.bndind + moved = fill(-1, nvtxs) + perm = Vector{Int}(undef, nvtxs) + tpwgts0 = trunc(Int, graph.tvwgt * ntpwgts[1]) + tpwgts = (tpwgts0, graph.tvwgt - tpwgts0) + mindiff = abs(tpwgts[1] - pwgts[1]) + from = pwgts[1] < tpwgts[1] ? 1 : 0 + to = (from + 1) % 2 + queue = rpqCreate(nvtxs) + irandArrayPermute!(nvtxs, perm, nvtxs ÷ 5, 1) + for ii in 0:(nvtxs - 1) + i = perm[ii + 1] + if where_[i + 1] == from && vwgt[i + 1] <= mindiff + rpqInsert(queue, i, Float32(ed[i + 1] - id[i + 1])) + end + end + mincut = graph.mincut + nbnd = graph.nbnd + for _ in 0:(nvtxs - 1) + higain = rpqGetTop(queue) + higain == -1 && break + pwgts[to + 1] + vwgt[higain + 1] > tpwgts[to + 1] && break + mincut -= ed[higain + 1] - id[higain + 1] + pwgts[to + 1] += vwgt[higain + 1] + pwgts[from + 1] -= vwgt[higain + 1] + where_[higain + 1] = to + moved[higain + 1] = 1 + id[higain + 1], ed[higain + 1] = ed[higain + 1], id[higain + 1] + if ed[higain + 1] == 0 && bndptr[higain + 1] != -1 && xadj[higain + 1] < xadj[higain + 2] + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + end + if ed[higain + 1] > 0 && bndptr[higain + 1] == -1 + bndind[nbnd + 1] = higain + bndptr[higain + 1] = nbnd + nbnd += 1 + end + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + kwgt = to == where_[k + 1] ? adjwgt[j + 1] : -adjwgt[j + 1] + id[k + 1] += kwgt + ed[k + 1] -= kwgt + if moved[k + 1] == -1 && where_[k + 1] == from && vwgt[k + 1] <= mindiff + rpqUpdate(queue, k, Float32(ed[k + 1] - id[k + 1])) + end + if ed[k + 1] == 0 && bndptr[k + 1] != -1 + bndind[bndptr[k + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[k + 1] + nbnd -= 1 + bndptr[k + 1] = -1 + elseif ed[k + 1] > 0 && bndptr[k + 1] == -1 + bndind[nbnd + 1] = k + bndptr[k + 1] = nbnd + nbnd += 1 + end + end + end + graph.mincut = mincut + graph.nbnd = nbnd + return nothing +end + + +# --------------------------------------------------------------------------- +# libmetis/initpart.c (iptype = EDGE path) + +function RandomBisection!(ctrl::CtrlT, graph::GraphT, ntpwgts, niparts::Int) + nvtxs = graph.nvtxs + vwgt = graph.vwgt + Allocate2WayPartitionMemory!(graph) + where_ = graph.where_ + bestwhere = Vector{Int}(undef, nvtxs) + perm = Vector{Int}(undef, nvtxs) + zeromaxpwgt = trunc(Int, ctrl.ubfactors * Float32(graph.tvwgt) * ntpwgts[1]) + bestcut = 0 + for inbfs in 0:(niparts - 1) + fill!(where_, 1) + if inbfs > 0 + irandArrayPermute!(nvtxs, perm, nvtxs ÷ 2, 1) + pw0 = 0 + for ii in 0:(nvtxs - 1) + i = perm[ii + 1] + if pw0 + vwgt[i + 1] < zeromaxpwgt + where_[i + 1] = 0 + pw0 += vwgt[i + 1] + pw0 > zeromaxpwgt && break + end + end + end + Compute2WayPartitionParams!(ctrl, graph) + Balance2Way!(ctrl, graph, ntpwgts) + FM_2WayRefine!(ctrl, graph, ntpwgts, 4) + if inbfs == 0 || bestcut > graph.mincut + bestcut = graph.mincut + copyto!(bestwhere, where_) + bestcut == 0 && break + end + end + graph.mincut = bestcut + copyto!(where_, bestwhere) + return nothing +end + +function GrowBisection!(ctrl::CtrlT, graph::GraphT, ntpwgts, niparts::Int) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + Allocate2WayPartitionMemory!(graph) + where_ = graph.where_ + bestwhere = Vector{Int}(undef, nvtxs) + queue = Vector{Int}(undef, nvtxs) + touched = Vector{Int}(undef, nvtxs) + onemaxpwgt = trunc(Int, ctrl.ubfactors * Float32(graph.tvwgt) * ntpwgts[2]) + oneminpwgt = trunc(Int, (1.0 / ctrl.ubfactors) * graph.tvwgt * ntpwgts[2]) + bestcut = 0 + for inbfs in 0:(niparts - 1) + fill!(where_, 1) + fill!(touched, 0) + pw0 = 0 + pw1 = graph.tvwgt + queue[1] = irandInRange(nvtxs) + touched[queue[1] + 1] = 1 + first = 0 + last = 1 + nleft = nvtxs - 1 + drain = false + while true + if first == last # empty: disconnected + (nleft == 0 || drain) && break + k = irandInRange(nleft) + i = 0 + while i < nvtxs + if touched[i + 1] == 0 + k == 0 && break + k -= 1 + end + i += 1 + end + queue[1] = i + touched[i + 1] = 1 + first = 0 + last = 1 + nleft -= 1 + end + i = queue[first + 1] + first += 1 + if pw0 > 0 && pw1 - vwgt[i + 1] < oneminpwgt + drain = true + continue + end + where_[i + 1] = 0 + pw0 += vwgt[i + 1] + pw1 -= vwgt[i + 1] + pw1 <= onemaxpwgt && break + drain = false + for j in xadj[i + 1]:(xadj[i + 2] - 1) + k = adjncy[j + 1] + if touched[k + 1] == 0 + queue[last + 1] = k + last += 1 + touched[k + 1] = 1 + nleft -= 1 + end + end + end + pw1 == 0 && (where_[irandInRange(nvtxs) + 1] = 1) + pw0 == 0 && (where_[irandInRange(nvtxs) + 1] = 0) + Compute2WayPartitionParams!(ctrl, graph) + Balance2Way!(ctrl, graph, ntpwgts) + FM_2WayRefine!(ctrl, graph, ntpwgts, ctrl.niter) + if inbfs == 0 || bestcut > graph.mincut + bestcut = graph.mincut + copyto!(bestwhere, where_) + bestcut == 0 && break + end + end + graph.mincut = bestcut + copyto!(where_, bestwhere) + return nothing +end + +function InitSeparator!(ctrl::CtrlT, graph::GraphT, niparts::Int) + ntpwgts = (0.5f0, 0.5f0) + # Setup2WayBalMultipliers + ctrl.pijbm = Float32[graph.invtvwgt / ntpwgts[1], graph.invtvwgt / ntpwgts[2]] + if graph.nedges == 0 + RandomBisection!(ctrl, graph, ntpwgts, niparts) + else + GrowBisection!(ctrl, graph, ntpwgts, niparts) + end + Compute2WayPartitionParams!(ctrl, graph) + ConstructSeparator!(ctrl, graph) + DBG[] && println("Initial Sep: ", graph.mincut) + return nothing +end + +# --------------------------------------------------------------------------- +# libmetis/srefine.c + +function Allocate2WayNodePartitionMemory!(graph::GraphT) + graph.pwgts = Vector{Int}(undef, 3) + graph.where_ = Vector{Int}(undef, graph.nvtxs) + graph.bndptr = Vector{Int}(undef, graph.nvtxs) + graph.bndind = Vector{Int}(undef, graph.nvtxs) + graph.nrinfo = Vector{NTuple{2, Int}}(undef, graph.nvtxs) + return nothing +end + +function Compute2WayNodePartitionParams!(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + where_ = graph.where_ + rinfo = graph.nrinfo + pwgts = graph.pwgts + fill!(pwgts, 0) + bndind = graph.bndind + bndptr = graph.bndptr + fill!(bndptr, -1) + nbnd = 0 + for i in 0:(nvtxs - 1) + me = where_[i + 1] + pwgts[me + 1] += vwgt[i + 1] + if me == 2 + bndind[nbnd + 1] = i + bndptr[i + 1] = nbnd + nbnd += 1 + e0 = 0 + e1 = 0 + for j in xadj[i + 1]:(xadj[i + 2] - 1) + other = where_[adjncy[j + 1] + 1] + if other != 2 + other == 0 ? (e0 += vwgt[adjncy[j + 1] + 1]) : (e1 += vwgt[adjncy[j + 1] + 1]) + end + end + rinfo[i + 1] = (e0, e1) + end + end + graph.mincut = pwgts[3] + graph.nbnd = nbnd + return nothing +end + +function Project2WayNodePartition!(ctrl::CtrlT, graph::GraphT) + cgraph = graph.coarser::GraphT + cwhere = cgraph.where_ + nvtxs = graph.nvtxs + cmap = graph.cmap + Allocate2WayNodePartitionMemory!(graph) + where_ = graph.where_ + for i in 0:(nvtxs - 1) + where_[i + 1] = cwhere[cmap[i + 1] + 1] + end + graph.coarser = nothing + Compute2WayNodePartitionParams!(ctrl, graph) + return nothing +end + +function Refine2WayNode!(ctrl::CtrlT, orggraph::GraphT, graph::GraphT) + if graph === orggraph + Compute2WayNodePartitionParams!(ctrl, graph) + else + while true + graph = graph.finer::GraphT + Project2WayNodePartition!(ctrl, graph) + FM_2WayNodeBalance!(ctrl, graph) + if ctrl.rtype == METIS_RTYPE_SEP2SIDED + FM_2WayNodeRefine2Sided!(ctrl, graph, ctrl.niter) + else + FM_2WayNodeRefine1Sided!(ctrl, graph, ctrl.niter) + end + graph === orggraph && break + end + end + return nothing +end + +# --------------------------------------------------------------------------- +# libmetis/separator.c ConstructSeparator (iptype = EDGE) + +function ConstructSeparator!(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + nbnd = graph.nbnd + bndind = graph.bndind + where_ = copy(graph.where_) + for i in 0:(nbnd - 1) + j = bndind[i + 1] + xadj[j + 2] - xadj[j + 1] > 0 && (where_[j + 1] = 2) + end + Allocate2WayNodePartitionMemory!(graph) + copyto!(graph.where_, where_) + Compute2WayNodePartitionParams!(ctrl, graph) + FM_2WayNodeRefine2Sided!(ctrl, graph, 1) + FM_2WayNodeRefine1Sided!(ctrl, graph, 4) + return nothing +end + + +# --------------------------------------------------------------------------- +# libmetis/sfm.c — node-separator FM + +function FM_2WayNodeRefine2Sided!(ctrl::CtrlT, graph::GraphT, niter::Int) + nvtxs = graph.nvtxs + xadj = graph.xadj + adjncy = graph.adjncy + vwgt = graph.vwgt + bndind = graph.bndind + bndptr = graph.bndptr + where_ = graph.where_ + pwgts = graph.pwgts + rinfo = graph.nrinfo + queues = (rpqCreate(nvtxs), rpqCreate(nvtxs)) + moved = Vector{Int}(undef, nvtxs) + swaps = Vector{Int}(undef, nvtxs) + mptr = Vector{Int}(undef, nvtxs + 1) + mind = Vector{Int}(undef, 2 * nvtxs) + mult = Float32(0.5 * ctrl.ubfactors) + badmaxpwgt = trunc(Int, mult * (pwgts[1] + pwgts[2] + pwgts[3])) + for pass in 0:(niter - 1) + fill!(moved, -1) + rpqReset(queues[1]) + rpqReset(queues[2]) + mincutorder = -1 + initcut = mincut = graph.mincut + nbnd = graph.nbnd + irandArrayPermute!(nbnd, swaps, nbnd, 1) + for ii in 0:(nbnd - 1) + i = bndind[swaps[ii + 1] + 1] + rpqInsert(queues[1], i, Float32(vwgt[i + 1] - rinfo[i + 1][2])) + rpqInsert(queues[2], i, Float32(vwgt[i + 1] - rinfo[i + 1][1])) + end + limit = ctrl.compress ? min(5 * nbnd, 400) : min(2 * nbnd, 300) + mptr[1] = 0 + nmind = 0 + mindiff = abs(pwgts[1] - pwgts[2]) + to = pwgts[1] < pwgts[2] ? 0 : 1 + nswaps = 0 + while nswaps < nvtxs + u0 = rpqSeeTopVal(queues[1]) + u1 = rpqSeeTopVal(queues[2]) + if u0 != -1 && u1 != -1 + g0 = vwgt[u0 + 1] - rinfo[u0 + 1][2] + g1 = vwgt[u1 + 1] - rinfo[u1 + 1][1] + to = g0 > g1 ? 0 : (g0 < g1 ? 1 : pass % 2) + u = (u0, u1) + if pwgts[to + 1] + vwgt[u[to + 1] + 1] > badmaxpwgt + to = (to + 1) % 2 + end + elseif u0 == -1 && u1 == -1 + break + elseif u0 != -1 && pwgts[1] + vwgt[u0 + 1] <= badmaxpwgt + to = 0 + elseif u1 != -1 && pwgts[2] + vwgt[u1 + 1] <= badmaxpwgt + to = 1 + else + break + end + other = (to + 1) % 2 + higain = rpqGetTop(queues[to + 1]) + moved[higain + 1] == -1 && rpqDelete(queues[other + 1], higain) + if nmind + xadj[higain + 2] - xadj[higain + 1] >= 2 * nvtxs - 1 + break + end + eo = rinfo[higain + 1][other + 1] + pwgts[3] -= vwgt[higain + 1] - eo + newdiff = abs(pwgts[to + 1] + vwgt[higain + 1] - (pwgts[other + 1] - eo)) + if pwgts[3] < mincut || (pwgts[3] == mincut && newdiff < mindiff) + mincut = pwgts[3] + mincutorder = nswaps + mindiff = newdiff + elseif nswaps - mincutorder > 2 * limit || + (nswaps - mincutorder > limit && pwgts[3] > 1.10 * mincut) + pwgts[3] += vwgt[higain + 1] - eo + break + end + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] # BNDDelete + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + pwgts[to + 1] += vwgt[higain + 1] + where_[higain + 1] = to + moved[higain + 1] = nswaps + swaps[nswaps + 1] = higain + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + if where_[k + 1] == 2 + oldgain = vwgt[k + 1] - rinfo[k + 1][to + 1] + rinfo[k + 1] = to == 0 ? (rinfo[k + 1][1] + vwgt[higain + 1], rinfo[k + 1][2]) : + (rinfo[k + 1][1], rinfo[k + 1][2] + vwgt[higain + 1]) + if moved[k + 1] == -1 || moved[k + 1] == -(2 + other) + rpqUpdate(queues[other + 1], k, Float32(oldgain - vwgt[higain + 1])) + end + elseif where_[k + 1] == other + bndind[nbnd + 1] = k # BNDInsert + bndptr[k + 1] = nbnd + nbnd += 1 + nmind += 1 + mind[nmind] = k + where_[k + 1] = 2 + pwgts[other + 1] -= vwgt[k + 1] + e0 = 0 + e1 = 0 + for jj in xadj[k + 1]:(xadj[k + 2] - 1) + kk = adjncy[jj + 1] + if where_[kk + 1] != 2 + where_[kk + 1] == 0 ? (e0 += vwgt[kk + 1]) : (e1 += vwgt[kk + 1]) + else + oldgain = vwgt[kk + 1] - rinfo[kk + 1][other + 1] + rinfo[kk + 1] = other == 0 ? (rinfo[kk + 1][1] - vwgt[k + 1], rinfo[kk + 1][2]) : + (rinfo[kk + 1][1], rinfo[kk + 1][2] - vwgt[k + 1]) + if moved[kk + 1] == -1 || moved[kk + 1] == -(2 + to) + rpqUpdate(queues[to + 1], kk, Float32(oldgain + vwgt[k + 1])) + end + end + end + rinfo[k + 1] = (e0, e1) + if moved[k + 1] == -1 + rpqInsert(queues[to + 1], k, Float32(vwgt[k + 1] - (other == 0 ? e0 : e1))) + moved[k + 1] = -(2 + to) + end + end + end + mptr[nswaps + 2] = nmind + nswaps += 1 + end + nswaps -= 1 + while nswaps > mincutorder + higain = swaps[nswaps + 1] + to = where_[higain + 1] + other = (to + 1) % 2 + pwgts[3] += vwgt[higain + 1] + pwgts[to + 1] -= vwgt[higain + 1] + where_[higain + 1] = 2 + bndind[nbnd + 1] = higain + bndptr[higain + 1] = nbnd + nbnd += 1 + e0 = 0 + e1 = 0 + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + if where_[k + 1] == 2 + rinfo[k + 1] = to == 0 ? (rinfo[k + 1][1] - vwgt[higain + 1], rinfo[k + 1][2]) : + (rinfo[k + 1][1], rinfo[k + 1][2] - vwgt[higain + 1]) + else + where_[k + 1] == 0 ? (e0 += vwgt[k + 1]) : (e1 += vwgt[k + 1]) + end + end + rinfo[higain + 1] = (e0, e1) + for j in mptr[nswaps + 1]:(mptr[nswaps + 2] - 1) + k = mind[j + 1] + where_[k + 1] = other + pwgts[other + 1] += vwgt[k + 1] + pwgts[3] -= vwgt[k + 1] + bndind[bndptr[k + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[k + 1] + nbnd -= 1 + bndptr[k + 1] = -1 + for jj in xadj[k + 1]:(xadj[k + 2] - 1) + kk = adjncy[jj + 1] + if where_[kk + 1] == 2 + rinfo[kk + 1] = other == 0 ? (rinfo[kk + 1][1] + vwgt[k + 1], rinfo[kk + 1][2]) : + (rinfo[kk + 1][1], rinfo[kk + 1][2] + vwgt[k + 1]) + end + end + end + nswaps -= 1 + end + graph.mincut = mincut + graph.nbnd = nbnd + (mincutorder == -1 || mincut >= initcut) && break + end + return nothing +end + +function FM_2WayNodeRefine1Sided!(ctrl::CtrlT, graph::GraphT, niter::Int) + nvtxs = graph.nvtxs + xadj = graph.xadj + adjncy = graph.adjncy + vwgt = graph.vwgt + bndind = graph.bndind + bndptr = graph.bndptr + where_ = graph.where_ + pwgts = graph.pwgts + rinfo = graph.nrinfo + queue = rpqCreate(nvtxs) + swaps = Vector{Int}(undef, nvtxs) + mptr = Vector{Int}(undef, nvtxs + 1) + mind = Vector{Int}(undef, 2 * nvtxs) + mult = Float32(0.5 * ctrl.ubfactors) + badmaxpwgt = trunc(Int, mult * (pwgts[1] + pwgts[2] + pwgts[3])) + to = pwgts[1] < pwgts[2] ? 1 : 0 + for pass in 0:(2 * niter - 1) + other = to + to = (to + 1) % 2 + rpqReset(queue) + mincutorder = -1 + initcut = mincut = graph.mincut + nbnd = graph.nbnd + irandArrayPermute!(nbnd, swaps, nbnd, 1) + for ii in 0:(nbnd - 1) + i = bndind[swaps[ii + 1] + 1] + rpqInsert(queue, i, Float32(vwgt[i + 1] - rinfo[i + 1][other + 1])) + end + limit = ctrl.compress ? min(5 * nbnd, 500) : min(3 * nbnd, 300) + mptr[1] = 0 + nmind = 0 + mindiff = abs(pwgts[1] - pwgts[2]) + nswaps = 0 + while nswaps < nvtxs + higain = rpqGetTop(queue) + higain == -1 && break + if nmind + xadj[higain + 2] - xadj[higain + 1] >= 2 * nvtxs - 1 + break + end + pwgts[to + 1] + vwgt[higain + 1] > badmaxpwgt && break + eo = rinfo[higain + 1][other + 1] + pwgts[3] -= vwgt[higain + 1] - eo + newdiff = abs(pwgts[to + 1] + vwgt[higain + 1] - (pwgts[other + 1] - eo)) + if pwgts[3] < mincut || (pwgts[3] == mincut && newdiff < mindiff) + mincut = pwgts[3] + mincutorder = nswaps + mindiff = newdiff + elseif nswaps - mincutorder > 3 * limit || + (nswaps - mincutorder > limit && pwgts[3] > 1.10 * mincut) + pwgts[3] += vwgt[higain + 1] - eo + break + end + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + pwgts[to + 1] += vwgt[higain + 1] + where_[higain + 1] = to + swaps[nswaps + 1] = higain + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + if where_[k + 1] == 2 + rinfo[k + 1] = to == 0 ? (rinfo[k + 1][1] + vwgt[higain + 1], rinfo[k + 1][2]) : + (rinfo[k + 1][1], rinfo[k + 1][2] + vwgt[higain + 1]) + elseif where_[k + 1] == other + bndind[nbnd + 1] = k + bndptr[k + 1] = nbnd + nbnd += 1 + nmind += 1 + mind[nmind] = k + where_[k + 1] = 2 + pwgts[other + 1] -= vwgt[k + 1] + e0 = 0 + e1 = 0 + for jj in xadj[k + 1]:(xadj[k + 2] - 1) + kk = adjncy[jj + 1] + if where_[kk + 1] != 2 + where_[kk + 1] == 0 ? (e0 += vwgt[kk + 1]) : (e1 += vwgt[kk + 1]) + else + rinfo[kk + 1] = other == 0 ? (rinfo[kk + 1][1] - vwgt[k + 1], rinfo[kk + 1][2]) : + (rinfo[kk + 1][1], rinfo[kk + 1][2] - vwgt[k + 1]) + rpqUpdate(queue, kk, Float32(vwgt[kk + 1] - rinfo[kk + 1][other + 1])) + end + end + rinfo[k + 1] = (e0, e1) + rpqInsert(queue, k, Float32(vwgt[k + 1] - (other == 0 ? e0 : e1))) + end + end + mptr[nswaps + 2] = nmind + nswaps += 1 + end + nswaps -= 1 + while nswaps > mincutorder + higain = swaps[nswaps + 1] + pwgts[3] += vwgt[higain + 1] + pwgts[to + 1] -= vwgt[higain + 1] + where_[higain + 1] = 2 + bndind[nbnd + 1] = higain + bndptr[higain + 1] = nbnd + nbnd += 1 + e0 = 0 + e1 = 0 + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + if where_[k + 1] == 2 + rinfo[k + 1] = to == 0 ? (rinfo[k + 1][1] - vwgt[higain + 1], rinfo[k + 1][2]) : + (rinfo[k + 1][1], rinfo[k + 1][2] - vwgt[higain + 1]) + else + where_[k + 1] == 0 ? (e0 += vwgt[k + 1]) : (e1 += vwgt[k + 1]) + end + end + rinfo[higain + 1] = (e0, e1) + for j in mptr[nswaps + 1]:(mptr[nswaps + 2] - 1) + k = mind[j + 1] + where_[k + 1] = other + pwgts[other + 1] += vwgt[k + 1] + pwgts[3] -= vwgt[k + 1] + bndind[bndptr[k + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[k + 1] + nbnd -= 1 + bndptr[k + 1] = -1 + for jj in xadj[k + 1]:(xadj[k + 2] - 1) + kk = adjncy[jj + 1] + if where_[kk + 1] == 2 + rinfo[kk + 1] = other == 0 ? (rinfo[kk + 1][1] + vwgt[k + 1], rinfo[kk + 1][2]) : + (rinfo[kk + 1][1], rinfo[kk + 1][2] + vwgt[k + 1]) + end + end + end + nswaps -= 1 + end + graph.mincut = mincut + graph.nbnd = nbnd + pass % 2 == 1 && (mincutorder == -1 || mincut >= initcut) && break + end + return nothing +end + +function FM_2WayNodeBalance!(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + adjncy = graph.adjncy + vwgt = graph.vwgt + bndind = graph.bndind + bndptr = graph.bndptr + where_ = graph.where_ + pwgts = graph.pwgts + rinfo = graph.nrinfo + mult = Float32(0.5 * ctrl.ubfactors) + badmaxpwgt = trunc(Int, mult * (pwgts[1] + pwgts[2])) + max(pwgts[1], pwgts[2]) < badmaxpwgt && return nothing + abs(pwgts[1] - pwgts[2]) < 3 * graph.tvwgt ÷ nvtxs && return nothing + to = pwgts[1] < pwgts[2] ? 0 : 1 + other = (to + 1) % 2 + queue = rpqCreate(nvtxs) + perm = Vector{Int}(undef, nvtxs) + moved = fill(-1, nvtxs) + nbnd = graph.nbnd + irandArrayPermute!(nbnd, perm, nbnd, 1) + for ii in 0:(nbnd - 1) + i = bndind[perm[ii + 1] + 1] + rpqInsert(queue, i, Float32(vwgt[i + 1] - rinfo[i + 1][other + 1])) + end + for _ in 0:(nvtxs - 1) + higain = rpqGetTop(queue) + higain == -1 && break + moved[higain + 1] = 1 + gain = vwgt[higain + 1] - rinfo[higain + 1][other + 1] + badmaxpwgt = trunc(Int, mult * (pwgts[1] + pwgts[2])) + pwgts[to + 1] > pwgts[other + 1] && break + gain < 0 && pwgts[other + 1] < badmaxpwgt && break + pwgts[to + 1] + vwgt[higain + 1] > badmaxpwgt && continue + pwgts[3] -= gain + bndind[bndptr[higain + 1] + 1] = bndind[nbnd] + bndptr[bndind[nbnd] + 1] = bndptr[higain + 1] + nbnd -= 1 + bndptr[higain + 1] = -1 + pwgts[to + 1] += vwgt[higain + 1] + where_[higain + 1] = to + for j in xadj[higain + 1]:(xadj[higain + 2] - 1) + k = adjncy[j + 1] + if where_[k + 1] == 2 + rinfo[k + 1] = to == 0 ? (rinfo[k + 1][1] + vwgt[higain + 1], rinfo[k + 1][2]) : + (rinfo[k + 1][1], rinfo[k + 1][2] + vwgt[higain + 1]) + elseif where_[k + 1] == other + bndind[nbnd + 1] = k + bndptr[k + 1] = nbnd + nbnd += 1 + where_[k + 1] = 2 + pwgts[other + 1] -= vwgt[k + 1] + e0 = 0 + e1 = 0 + for jj in xadj[k + 1]:(xadj[k + 2] - 1) + kk = adjncy[jj + 1] + if where_[kk + 1] != 2 + where_[kk + 1] == 0 ? (e0 += vwgt[kk + 1]) : (e1 += vwgt[kk + 1]) + else + oldgain = vwgt[kk + 1] - rinfo[kk + 1][other + 1] + rinfo[kk + 1] = other == 0 ? (rinfo[kk + 1][1] - vwgt[k + 1], rinfo[kk + 1][2]) : + (rinfo[kk + 1][1], rinfo[kk + 1][2] - vwgt[k + 1]) + moved[kk + 1] == -1 && rpqUpdate(queue, kk, Float32(oldgain + vwgt[k + 1])) + end + end + rinfo[k + 1] = (e0, e1) + rpqInsert(queue, k, Float32(vwgt[k + 1] - (other == 0 ? e0 : e1))) + end + end + end + graph.mincut = pwgts[3] + graph.nbnd = nbnd + return nothing +end + + +# --------------------------------------------------------------------------- +# libmetis/ometis.c — bisection drivers, split, MMD leaves, NodeND + +function MlevelNodeBisectionMultiple!(ctrl::CtrlT, graph::GraphT) + if ctrl.nseps == 1 || graph.nvtxs < (ctrl.compress ? 1000 : 2000) + MlevelNodeBisectionL2!(ctrl, graph, LARGENIPARTS) + return nothing + end + bestwhere = Vector{Int}(undef, graph.nvtxs) + mincut = graph.tvwgt + for i in 0:(ctrl.nseps - 1) + MlevelNodeBisectionL2!(ctrl, graph, LARGENIPARTS) + if i == 0 || graph.mincut < mincut + mincut = graph.mincut + i < ctrl.nseps - 1 && copyto!(bestwhere, graph.where_) + end + mincut == 0 && break + end + if mincut != graph.mincut + copyto!(graph.where_, bestwhere) + Compute2WayNodePartitionParams!(ctrl, graph) + end + return nothing +end + +function MlevelNodeBisectionL2!(ctrl::CtrlT, graph::GraphT, niparts::Int) + if graph.nvtxs < 5000 + MlevelNodeBisectionL1!(ctrl, graph, niparts) + return nothing + end + nruns = 5 + ctrl.CoarsenTo = max(100, graph.nvtxs ÷ 30) + t = @elapsed (cgraph = CoarsenGraphNlevels(ctrl, graph, 4)) + prof!(:coarsen, t) + bestwhere = Vector{Int}(undef, cgraph.nvtxs) + mincut = graph.tvwgt + for i in 0:(nruns - 1) + MlevelNodeBisectionL1!(ctrl, cgraph, trunc(Int, 0.7 * niparts)) + if i == 0 || cgraph.mincut < mincut + mincut = cgraph.mincut + i < nruns - 1 && copyto!(bestwhere, cgraph.where_) + end + mincut == 0 && break + end + if mincut != cgraph.mincut + copyto!(cgraph.where_, bestwhere) + Compute2WayNodePartitionParams!(ctrl, cgraph) + end + t = @elapsed Refine2WayNode!(ctrl, graph, cgraph) + prof!(:refine_l2, t) + return nothing +end + +function MlevelNodeBisectionL1!(ctrl::CtrlT, graph::GraphT, niparts::Int) + ctrl.CoarsenTo = graph.nvtxs ÷ 8 + ctrl.CoarsenTo > 100 && (ctrl.CoarsenTo = 100) + ctrl.CoarsenTo < 40 && (ctrl.CoarsenTo = 40) + t = @elapsed (cgraph = CoarsenGraph(ctrl, graph)) + prof!(:coarsen, t) + niparts = max(1, cgraph.nvtxs <= ctrl.CoarsenTo ? niparts ÷ 2 : niparts) + t = @elapsed InitSeparator!(ctrl, cgraph, niparts) + prof!(:initsep, t) + t = @elapsed Refine2WayNode!(ctrl, graph, cgraph) + prof!(:refine, t) + return nothing +end + +# SplitGraphOrder: tri-section into the left/right subgraphs (separator out) +function SplitGraphOrder(ctrl::CtrlT, graph::GraphT) + nvtxs = graph.nvtxs + xadj = graph.xadj + vwgt = graph.vwgt + adjncy = graph.adjncy + label = graph.label + where_ = graph.where_ + bndptr = graph.bndptr + bndind = graph.bndind + rename = Vector{Int}(undef, nvtxs) + snvtxs = [0, 0, 0] + snedges = [0, 0, 0] + for i in 0:(nvtxs - 1) + k = where_[i + 1] + rename[i + 1] = snvtxs[k + 1] + snvtxs[k + 1] += 1 + snedges[k + 1] += xadj[i + 2] - xadj[i + 1] + end + lgraph = SetupSplitGraph(graph, snvtxs[1], snedges[1]) + rgraph = SetupSplitGraph(graph, snvtxs[2], snedges[2]) + # mark the boundary vertices of the two partitions through bndptr + for ii in 0:(graph.nbnd - 1) + i = bndind[ii + 1] + for j in xadj[i + 1]:(xadj[i + 2] - 1) + bndptr[adjncy[j + 1] + 1] = 1 + end + end + sn = [0, 0] + se = [0, 0] + lgraph.xadj[1] = 0 + rgraph.xadj[1] = 0 + sg = (lgraph, rgraph) + for i in 0:(nvtxs - 1) + mypart = where_[i + 1] + mypart == 2 && continue + g = sg[mypart + 1] + istart = xadj[i + 1] + iend = xadj[i + 2] + if bndptr[i + 1] == -1 # interior vertex + l = se[mypart + 1] + for j in istart:(iend - 1) + l += 1 + g.adjncy[l] = adjncy[j + 1] + end + se[mypart + 1] = l + else + l = se[mypart + 1] + for j in istart:(iend - 1) + k = adjncy[j + 1] + if where_[k + 1] == mypart + l += 1 + g.adjncy[l] = k + end + end + se[mypart + 1] = l + end + g.vwgt[sn[mypart + 1] + 1] = vwgt[i + 1] + g.label[sn[mypart + 1] + 1] = label[i + 1] + sn[mypart + 1] += 1 + g.xadj[sn[mypart + 1] + 1] = se[mypart + 1] + end + for mypart in 0:1 + g = sg[mypart + 1] + for i in 1:se[mypart + 1] + g.adjwgt[i] = 1 + g.adjncy[i] = rename[g.adjncy[i] + 1] + end + g.nvtxs = sn[mypart + 1] + g.nedges = se[mypart + 1] + SetupGraph_tvwgt!(g) + end + return lgraph, rgraph +end + +# MMDOrder: multiple minimum degree on a small piece, through the exported +# libmetis__genmmd of the shipped binary (pure arrays; exact by construction) +const GENMMD_PTR = Ref{Ptr{Cvoid}}(C_NULL) # set via init_libmetis!(path) + +function init_libmetis!(path::AbstractString) + h = Libdl.dlopen(path) + GENMMD_PTR[] = Libdl.dlsym(h, :libmetis__genmmd) + return nothing +end + +function MMDOrder!(ctrl::CtrlT, graph::GraphT, order::Vector{Int}, lastvtx::Int) + t0 = time() + nvtxs = graph.nvtxs + # 1-based CSR copies for the Fortran-style genmmd + xadj1 = Vector{Int32}(undef, nvtxs + 1) + for i in 1:(nvtxs + 1) + xadj1[i] = Int32(graph.xadj[i] + 1) + end + adjncy1 = Vector{Int32}(undef, max(graph.nedges, 1)) + for i in 1:graph.nedges + adjncy1[i] = Int32(graph.adjncy[i] + 1) + end + perm = Vector{Int32}(undef, nvtxs + 5) + iperm = Vector{Int32}(undef, nvtxs + 5) + head = Vector{Int32}(undef, nvtxs + 5) + qsize = Vector{Int32}(undef, nvtxs + 5) + list = Vector{Int32}(undef, nvtxs + 5) + marker = Vector{Int32}(undef, nvtxs + 5) + nofsub = Ref{Int32}(0) + ccall(GENMMD_PTR[], Cvoid, + (Int32, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Int32, + Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Int32, Ref{Int32}), + Int32(nvtxs), xadj1, adjncy1, iperm, perm, Int32(1), + head, qsize, list, marker, Int32(IDX_MAX), nofsub) + label = graph.label + firstvtx = lastvtx - nvtxs + for i in 0:(nvtxs - 1) + order[label[i + 1] + 1] = firstvtx + Int(iperm[i + 1]) - 1 + end + prof!(:mmd, time() - t0) + return nothing +end + +function MlevelNestedDissection!(ctrl::CtrlT, graph::GraphT, order::Vector{Int}, lastvtx::Int) + MlevelNodeBisectionMultiple!(ctrl, graph) + DBG[] && println("Nvtxs: ", lpad(graph.nvtxs, 6), ", [", lpad(graph.pwgts[1], 6), + lpad(graph.pwgts[2], 7), lpad(graph.pwgts[3], 7), "]") + nbnd = graph.nbnd + bndind = graph.bndind + label = graph.label + for i in 0:(nbnd - 1) + lastvtx -= 1 + order[label[bndind[i + 1] + 1] + 1] = lastvtx + end + t = @elapsed ((lgraph, rgraph) = SplitGraphOrder(ctrl, graph)) + prof!(:split, t) + if lgraph.nvtxs > MMDSWITCH && lgraph.nedges > 0 + MlevelNestedDissection!(ctrl, lgraph, order, lastvtx - rgraph.nvtxs) + else + MMDOrder!(ctrl, lgraph, order, lastvtx - rgraph.nvtxs) + end + if rgraph.nvtxs > MMDSWITCH && rgraph.nedges > 0 + MlevelNestedDissection!(ctrl, rgraph, order, lastvtx) + else + MMDOrder!(ctrl, rgraph, order, lastvtx) + end + return nothing +end + +""" + NodeND(nvtxs, xadj, adjncy; vwgt = nothing, seed = -1) -> (perm, iperm) + +libmetis METIS_NodeND with OMETIS defaults, 0-based CSR in, 0-based perm and +iperm out. ccorder and pfactor are not translated (defaults are off). +""" +function NodeND(nvtxs::Int, xadj::Vector{Int}, adjncy::Vector{Int}; + vwgt::Union{Nothing, Vector{Int}} = nothing, seed::Int = -1) + ctrl = CtrlT(METIS_CTYPE_SHEM, METIS_RTYPE_SEP1SIDED, METIS_IPTYPE_EDGE, + 1, 10, 200, true, false, 0.0, false, seed, 100, 0, + Float32(Float64(Float32(1.0 + 0.001 * 200)) + 0.0000499), Float32[]) + InitRandom(ctrl.seed) + nnvtxs = 0 + graph = nothing + cptr = Int[] + cind = Int[] + if ctrl.compress + cptr = Vector{Int}(undef, nvtxs + 1) + cind = Vector{Int}(undef, nvtxs) + graph = CompressGraph(nvtxs, xadj, adjncy, vwgt, cptr, cind) + if graph === nothing + ctrl.compress = false + else + nnvtxs = graph.nvtxs + cfactor = 1.0 * nvtxs / nnvtxs + cfactor > 1.5 && ctrl.nseps == 1 && (ctrl.nseps = 2) + end + end + graph === nothing && (graph = SetupGraph(nvtxs, xadj, adjncy, vwgt)) + iperm = Vector{Int}(undef, nvtxs) + perm = Vector{Int}(undef, nvtxs) + MlevelNestedDissection!(ctrl, graph, iperm, graph.nvtxs) + if ctrl.compress # uncompress the ordering + for i in 0:(nnvtxs - 1) + perm[iperm[i + 1] + 1] = i + end + l = 0 + for ii in 0:(nnvtxs - 1) + i = perm[ii + 1] + for j in cptr[i + 1]:(cptr[i + 2] - 1) + iperm[cind[j + 1] + 1] = l + l += 1 + end + end + end + for i in 0:(nvtxs - 1) + perm[iperm[i + 1] + 1] = i + end + return perm, iperm +end + +""" + nodend_exact(A::SparseMatrixCSC; seed = -1) -> (perm, iperm) 1-based + +The `Metis.permutation(A)` replacement: symmetric pattern of `A`, diagonal +dropped, vwgt = NULL, all-default options. +""" +function nodend_exact(A; seed::Int = -1) + n = size(A, 1) + xadj = Vector{Int}(undef, n + 1) + adjncy = Int[] + rows = SparseArrays.rowvals(A) + xadj[1] = 0 + for j in 1:n + for p in SparseArrays.nzrange(A, j) + i = rows[p] + i != j && push!(adjncy, i - 1) + end + xadj[j + 1] = length(adjncy) + end + perm, iperm = NodeND(n, xadj, adjncy; seed) + return perm .+ 1, iperm .+ 1 +end + +end # module diff --git a/bench/ordering/gpu_nodend.jl b/bench/ordering/gpu_nodend.jl new file mode 100644 index 0000000..a969a7a --- /dev/null +++ b/bench/ordering/gpu_nodend.jl @@ -0,0 +1,31 @@ +# GPU-hybrid METIS NodeND — loader for the ordering prototype (bench-level, +# stock src untouched). +# +# MetisExact.jl bitwise-faithful Julia translation of METIS 5.1.0's NodeND +# path (all-default options, the PR #107 configuration); +# validated perm- and iperm-identical to libmetis on lap2d/ +# lap3d and the case78484 KKT dumps up to n = 2.09M. +# GpuND.jl KernelAbstractions device kernels for the pipeline stages +# (hashed-handshake HEM matching incl. 2-hop, device-resident +# contraction), stage-validated against the exact port. +# Hybrid.jl nodend_hybrid(backend, A): device coarsening chain + exact +# host initial separator / FM refinement / MMD, threaded +# sub-threshold subtrees. Fill quality 1.003–1.015 nnz_L vs +# exact METIS on the benchmark set. +# +# Usage (any KA backend; CUDABackend for the GPU, CPU() works for smoke): +# include("bench/ordering/gpu_nodend.jl") +# perm, iperm = Hybrid.nodend_hybrid(backend, A) # A symmetric sparse +# +# MetisExact ccalls libc rand() and the binary's genmmd, so libmetis must be +# initialized once per process (Metis.jl's jll carries it): + +using Metis + +include(joinpath(@__DIR__, "MetisExact.jl")) +include(joinpath(@__DIR__, "GpuND.jl")) +include(joinpath(@__DIR__, "GpuRefine.jl")) +include(joinpath(@__DIR__, "Hybrid.jl")) +include(joinpath(@__DIR__, "GpuBatch.jl")) + +MetisExact.init_libmetis!(Metis.libmetis) diff --git a/bench/ordering/nodend_smoke.jl b/bench/ordering/nodend_smoke.jl new file mode 100644 index 0000000..25377d4 --- /dev/null +++ b/bench/ordering/nodend_smoke.jl @@ -0,0 +1,51 @@ +# Smoke + quality gate for the GPU-hybrid NodeND prototype: a 2-D Laplacian +# ordered by the exact port and by the hybrid, scored by nnz_L through SDS's +# own symbolic analysis. Passes = both perms valid and the hybrid within 5% +# fill of the exact port (measured 1.003–1.015 on the benchmark set; 5% is +# the alarm threshold, not the expectation). +# julia +1.13 --project=bench bench/ordering/nodend_smoke.jl [grid=60] +using SparseDirectSolver, SparseArrays, LinearAlgebra, Printf +const SDS = SparseDirectSolver +const KA = SDS.KernelAbstractions + +include(joinpath(@__DIR__, "gpu_nodend.jl")) + +laplacian2d(nx, ny) = begin + Tx = spdiagm(-1 => -ones(nx - 1), 0 => 4ones(nx), 1 => -ones(nx - 1)) + Ty = spdiagm(-1 => -ones(ny - 1), 1 => -ones(ny - 1)) + Ix = sparse(1.0I, nx, nx); Iy = sparse(1.0I, ny, ny) + kron(Iy, Tx) + kron(Ty, Ix) +end + +function nnzL(A, perm) + n = size(A, 1) + Lt = sparse(transpose(tril(A))) + P = SDS.SymmetricPattern(Vector{Int}(Lt.colptr), Vector{Int}(Lt.rowval), n, + SDS.STRUCTURE_SPD; view = SDS.VIEW_LOWER, index = SDS.INDEX_ONE) + Int(SDS.evaluate_ordering(P, perm).nnz_L) +end + +g = isempty(ARGS) ? 60 : parse(Int, ARGS[1]) +A = laplacian2d(g, g) +n = size(A, 1) +println("lap2d $(g)x$(g) n = $n") + +pe, _ = MetisExact.nodend_exact(A) +@assert isperm(pe) +ph, _ = Hybrid.nodend_hybrid(KA.CPU(), A) +@assert isperm(ph) + +le = nnzL(A, pe) +lh = nnzL(A, ph) +@printf "nnz_L exact %d hybrid %d ratio %.4f\n" le lh lh / le +@assert lh / le < 1.05 "hybrid fill regressed beyond the 5% alarm threshold" + +# device-refinement path (only engages above Hybrid.GPU_MIN_NVTXS) +Hybrid.GPU_REFINE[] = true +pg, _ = Hybrid.nodend_hybrid(KA.CPU(), A) +Hybrid.GPU_REFINE[] = false +@assert isperm(pg) +lg = nnzL(A, pg) +@printf "nnz_L gpu-refine %d ratio vs exact %.4f (engaged: %s)\n" lg lg / le (n >= Hybrid.GPU_MIN_NVTXS[] ? "yes" : "no — below threshold") +n >= Hybrid.GPU_MIN_NVTXS[] && @assert lg / le < 1.08 "gpu-refine fill beyond alarm threshold" +println("SMOKE OK")