Vendor BQM and add paper benchmarks

Project.toml

@@ -4,14 +4,14 @@ version = "0.2.0"
 authors = ["François Pacaud <francoispacaud8@gmail.com>", "Alexis Montoison <amontoison@anl.gov>"]
 
 [deps]
-BatchQuadraticModels = "412afcf0-3ec6-4826-be17-8d792afa05a8"
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseMatricesCOO = "fa32481b-f100-4b48-8dc8-c62f61b13870"
 
 [weakdeps]
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
@@ -21,16 +21,13 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 MadNLPGPU = "d72a61cc-809d-412f-99be-fd81f4b8a598"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 
-[sources]
-BatchQuadraticModels = {rev = "main", url = "https://github.com/klamike/BatchQuadraticModels.jl.git"}
-
 [extensions]
 MadIPMCUDAExt = ["Atomix", "CUDA", "CUDSS", "KernelAbstractions", "MadNLPGPU"]
 MadIPMMathOptInterfaceExt = "MathOptInterface"
 
 [compat]
+Adapt = "4"
 Atomix = "1"
-BatchQuadraticModels = "0.1"
 CUDA = "5.4.0"
 CUDSS = "0.6"
 KernelAbstractions = "0.9"
@@ -40,10 +37,10 @@ MadNLP = "0.9.1"
 MadNLPGPU = "0.8"
 MadNLPTests = "0.6"
 MathOptInterface = "1.48"
-NLPModels = "0.21.10"
+NLPModels = "0.21.12"
 Printf = "1.10"
-QuadraticModels = "0.9.14"
 SparseArrays = "1.10"
+SparseMatricesCOO = "0.2.4"
 Test = "1.10"
 julia = "1.10"
 

benchmark/1-scalability.jl

@@ -0,0 +1,120 @@
+
+include("common.jl")
+
+const NETLIB_PATH = fetch_netlib()
+const NNZJ_THRESHOLD = 5_000
+const WARMUP_INSTANCE = "ADLITTLE.SIF"
+
+@memoize function load_instance(case)
+    qpdat = readqps(joinpath(NETLIB_PATH, case))
+    return qps_model(qpdat)
+end
+
+function warmup(instance)
+    qp = load_instance(instance)
+    _warmup(qp; linear_solver=Ma57Solver)
+    return
+end
+
+function select_netlib()
+    cases = filter(x -> endswith(x, ".SIF"), readdir(NETLIB_PATH))
+    selected = String[]
+    for case in cases
+        try
+            qp = load_instance(case)
+            # Select only medium-sized instances
+            if NLPModels.get_nnzj(qp) <= NNZJ_THRESHOLD
+                push!(selected, case)
+            end
+        catch ex
+            println("Fail to load $(case)")
+        end
+    end
+    return selected
+end
+
+# Columns: nvar, ncon, nnzj of the loaded instance; then, per batch size b,
+#   CPU on instances 1:b: converged count, mean iter, summed init time, summed
+#       solve time (one CPU, sequentially), max init time, max solve time
+#       (b CPUs with perfect scaling);
+#   GPU solving the same instances 1:b as one batch: converged count, mean iter,
+#       init time, solve time.
+# Both sides solve the same presolved, scaled, standard-form instances; all
+# times are wall clock.
+function benchmark_scalability(cases, batches; options...)
+    shift = 3
+    m = shift + 10*length(batches)
+    results = zeros(length(cases), m)
+
+    for (k, case) in enumerate(cases)
+        @info case
+        refresh_memory()
+        qp = load_instance(case)
+        results[k, 1] = NLPModels.get_nvar(qp)
+        results[k, 2] = NLPModels.get_ncon(qp)
+        results[k, 3] = NLPModels.get_nnzj(qp)
+        # Test pure scalability, do not change cost vector here.
+        qps = build_qps(qp, batches[end]; shift_c=false)
+        # CPU: every instance of the largest batch, sequentially
+        cpu = timed_cpu_sequential(qps; linear_solver=Ma57Solver, options...)
+        for (l, batch) in enumerate(batches)
+            results[k, shift+10*(l-1) .+ (1:6)] .= cpu_summary(cpu..., batch)
+            # GPU: the same instances as one batch
+            refresh_memory()
+            gpu_bnlp = to_gpu(ObjRHSBatchQuadraticModel(qps[1:batch]))
+            _, stats, t_init, t_solve = timed_gpu_solve(gpu_bnlp; cudss_pivot_epsilon=1e-8, options...)
+            results[k, shift+10*(l-1) .+ (7:10)] .= gpu_summary(stats, t_init, t_solve)
+        end
+    end
+
+    return [cases results]
+end
+
+function parse_args(args::Vector{String})
+    # Default options
+    max_batch = 12
+    device = nothing
+    tol = 1e-6
+    for arg in args
+        if startswith(arg, "--tol=")
+            tol = parse(Float64, split(arg, "=")[2])
+        elseif startswith(arg, "--max-batch=")
+            max_batch = parse(Int, split(arg, "=")[2])
+        elseif startswith(arg, "--device=")
+            device = parse(Int, split(arg, "=")[2])
+        end
+    end
+    return (
+        max_batch=max_batch,
+        tol=tol,
+        device=device,
+    )
+end
+
+function @main(args::Vector{String})
+    pargs = parse_args(args)
+
+    # Set-up device
+    if !isnothing(pargs.device)
+        CUDA.device!(pargs.device)
+    end
+
+    @info "Warmup"
+    warmup(WARMUP_INSTANCE)
+
+    batches = [2^i for i in 0:pargs.max_batch]
+    cases = select_netlib()
+    @info "#instances: $(length(cases))"
+    results = benchmark_scalability(
+        cases,
+        batches;
+        print_level=MadNLP.ERROR,
+        tol=pargs.tol,
+        max_iter=500,
+        regularization = MadIPM.FixedRegularization(1e-10, -1e-10),
+    )
+    mkpath("results")
+    writedlm(joinpath("results", "1-scalability-netlib.csv"), results)
+end
+
+

benchmark/2-performance.jl

@@ -0,0 +1,206 @@
+
+include("common.jl")
+
+using MIPLIB
+
+
+const NETLIB_PATH = fetch_netlib()
+const MIPLIB_INSTANCES = "miplib_problems.txt"
+const WARMUP_INSTANCE = "ADLITTLE.SIF"
+
+function select_netlib_instance()
+    cases = filter(x -> endswith(x, ".SIF"), readdir(NETLIB_PATH))
+    selected = String[]
+    for case in cases
+        try
+            qp = readqps(joinpath(NETLIB_PATH, case))
+            push!(selected, case)
+        catch ex
+            println("Fail to load $(case)")
+        end
+    end
+    return selected
+end
+
+function select_miplib_instance()
+    return [
+        "30n20b8.mps.gz",
+        "aflow40b.mps.gz",
+        "ash608gpia-3col.mps.gz",
+        "biella1.mps.gz",
+        "binkar10_1.mps.gz",
+        "bnatt350.mps.gz",
+        "core2536-691.mps.gz",
+        "cov1075.mps.gz",
+        "eil33-2.mps.gz",
+        "eilB101.mps.gz",
+        "enlight13.mps.gz",
+        "enlight14.mps.gz",
+        "glass4.mps.gz",
+        "gmu-35-40.mps.gz",
+        "iis-100-0-cov.mps.gz",
+        "iis-bupa-cov.mps.gz",
+        "iis-pima-cov.mps.gz",
+        "lectsched-4-obj.mps.gz",
+        "m100n500k4r1.mps.gz",
+        "macrophage.mps.gz",
+        "map18.mps.gz",               # slow 4
+        "map20.mps.gz",               # slow 4
+        "mik-250-1-100-1.mps.gz",
+        "mine-166-5.mps.gz",
+        "mine-90-10.mps.gz",
+        "n3div36.mps.gz",
+        "neos-1109824.mps.gz",
+        "neos13.mps.gz",
+        "neos18.mps.gz",
+        "neos-934278.mps.gz",
+        "noswot.mps.gz",
+        "pg5_34.mps.gz",
+        "pw-myciel4.mps.gz",
+        "qiu.mps.gz",
+        "rail507.mps.gz",
+        "reblock67.mps.gz",
+        "rmatr100-p10.mps.gz",
+        "rmatr100-p5.mps.gz",
+        "rmine6.mps.gz",
+        "sp98ic.mps.gz",
+        "tanglegram1.mps.gz",
+        "tanglegram2.mps.gz",
+        "timtab1.mps.gz",
+    ]
+end
+
+@memoize function load_netlib_instance(case)
+    qpdat = readqps(joinpath(NETLIB_PATH, case))
+    return qps_model(qpdat)
+end
+
+@memoize function load_miplib_instance(case)
+    return qps_model(MIPLIB.miplib2010_data(case))
+end
+
+# Columns: nvar, ncon, nnzj of the loaded instance; then, per batch size b,
+#   CPU on instances 1:b: converged count, mean iter, summed init time, summed
+#       solve time (one CPU, sequentially), max init time, max solve time
+#       (b CPUs with perfect scaling);
+#   GPU solving the same instances 1:b as one batch: converged count, mean iter,
+#       init time, solve time.
+# Both sides solve the same presolved, scaled, standard-form instances, cost
+# perturbations included; all times are wall clock. Failures are marked -1.
+function benchmark_lps(cases, batches, load_instance; options...)
+    shift = 3
+    m = shift + 10*length(batches)
+    results = zeros(length(cases), m)
+
+    for (k, case) in enumerate(cases)
+        @info case
+        refresh_memory()
+        # Load instance
+        qp = load_instance(case)
+        results[k, 1] = NLPModels.get_nvar(qp)
+        results[k, 2] = NLPModels.get_ncon(qp)
+        results[k, 3] = NLPModels.get_nnzj(qp)
+        qps = try
+            build_qps(qp, batches[end])
+        catch ex
+            println("$(case) fails in presolve with message $(ex)")
+            results[k, shift+1:end] .= -1
+            continue
+        end
+        # CPU: every instance of the largest batch, sequentially
+        cpu = try
+            timed_cpu_sequential(qps; linear_solver=Ma27Solver, options...)
+        catch ex
+            println("$(case) fails on CPU with message $(ex)")
+            nothing
+        end
+        for (l, batch) in enumerate(batches)
+            cpu_cols = shift+10*(l-1) .+ (1:6)
+            gpu_cols = shift+10*(l-1) .+ (7:10)
+            if cpu === nothing
+                results[k, cpu_cols] .= -1
+            else
+                results[k, cpu_cols] .= cpu_summary(cpu..., batch)
+            end
+            # GPU: the same instances as one batch
+            try
+                refresh_memory()
+                gpu_bnlp = to_gpu(ObjRHSBatchQuadraticModel(qps[1:batch]))
+                _, stats, t_init, t_solve = timed_gpu_solve(gpu_bnlp; cudss_pivot_epsilon=1e-8, options...)
+                results[k, gpu_cols] .= gpu_summary(stats, t_init, t_solve)
+            catch ex
+                println("$(case) fails on GPU with message $(ex)")
+                results[k, gpu_cols] .= -1
+            end
+        end
+    end
+    return [cases results]
+end
+
+function parse_args(args::Vector{String})
+    # Default options
+    max_batch = 7
+    device = nothing
+    tol = 1e-6
+    benchmark = :netlib
+    for arg in args
+        if startswith(arg, "--tol=")
+            tol = parse(Float64, split(arg, "=")[2])
+        elseif startswith(arg, "--max-batch=")
+            max_batch = parse(Int, split(arg, "=")[2])
+        elseif startswith(arg, "--device=")
+            device = parse(Int, split(arg, "=")[2])
+        elseif startswith(arg, "--benchmark=")
+            benchmark = Symbol(split(arg, "=")[2])
+        end
+    end
+    return (
+        max_batch=max_batch,
+        tol=tol,
+        device=device,
+        benchmark=benchmark,
+    )
+end
+
+function @main(args::Vector{String})
+    pargs = parse_args(args)
+
+    # Set-up device
+    if !isnothing(pargs.device)
+        CUDA.device!(pargs.device)
+    end
+
+    @info "Warmup"
+    _warmup(load_netlib_instance(WARMUP_INSTANCE); linear_solver=Ma27Solver)
+
+    batches = [2^i for i in 0:pargs.max_batch]
+    if pargs.benchmark == :netlib
+        cases = select_netlib_instance()
+        mkpath("results")
+        results = benchmark_lps(
+            cases,
+            batches,
+            load_netlib_instance;
+            print_level=MadNLP.ERROR,
+            tol=pargs.tol,
+            max_iter=300,
+            regularization = MadIPM.FixedRegularization(1e-8, -1e-8),
+        )
+        writedlm(joinpath("results", "2-benchmark-netlib.csv"), results)
+    elseif pargs.benchmark == :miplib
+        cases = select_miplib_instance()
+        mkpath("results")
+        results = benchmark_lps(
+            cases,
+            batches,
+            load_miplib_instance;
+            print_level=MadNLP.ERROR,
+            tol=pargs.tol,
+            max_iter=300,
+            regularization = MadIPM.FixedRegularization(1e-8, -1e-8),
+        )
+        writedlm(joinpath("results", "2-benchmark-miplib.csv"), results)
+    end
+    return
+end
+