From cc505a549ba1f66dadbd8a2a41ae6e8f0673cfbd Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 1 May 2023 18:23:37 -0400 Subject: [PATCH 1/2] [HSL_MC66] Add a Julia interface --- src/hsl_mc66.jl | 98 +++++++++++++++++++++++++++++++++++++++++++ test/test_hsl_mc66.jl | 16 +++++++ 2 files changed, 114 insertions(+) create mode 100644 src/hsl_mc66.jl create mode 100644 test/test_hsl_mc66.jl diff --git a/src/hsl_mc66.jl b/src/hsl_mc66.jl new file mode 100644 index 0000000..ad9c6a8 --- /dev/null +++ b/src/hsl_mc66.jl @@ -0,0 +1,98 @@ +export mc66, mc66_control + +mutable struct fa14_seed + ix::Int32 +end + +fa14_seed() = fa14_seed(65535) + +mutable struct mc66_control{T <: BLAS.BlasReal} + lp::Int32 + wp::Int32 + mp::Int32 + print_level::Int32 + kl_aggressive::T + max_imbalance::T + mglevel::Int32 + grid_rdc_fac::T + coarsest_size::Int32 + coarsen_scheme::Int32 + num_coarsest_kl::Int32 +end + +function mc66_control{T}(; lp=6, wp=6, mp=6, print_level=-1, kl_aggressive=-1, max_imbalance=0.01, + mglevel=typemax(Int32), grid_rdc_fac=0.75, coarsest_size=100, + coarsen_scheme=1, num_coarsest_kl=4) where {T <: BLAS.BlasReal} + + control = mc66_control{T}(lp, wp, mp, print_level, kl_aggressive, max_imbalance, mglevel, + grid_rdc_fac, coarsest_size, coarsen_scheme, num_coarsest_kl) + return control +end + +# Thanks to `nm -D libhsl_mc66.so`! +function mc66s(m, n, nz, irn, jcn, nblocks, control, seed, row_order, info, rowptr, column_order, colptr, netcut, rowdiff, kblocks) + ccall((:__hsl_mc66_simple_MOD_monet, libhsl_mc66), + Cvoid, + (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32}, Ref{Int32}, Ref{mc66_control{Float32}}, Ref{fa14_seed}, Ptr{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32} , Ptr{Int32}, Ref{Int32}, Ref{Float32}, Ref{Int32}), + m , n , nz , irn , jcn , nblocks , control , seed , row_order , info , rowptr , column_order, colptr , netcut , rowdiff , kblocks ) +end + +function mc66s_print_message(info, control) + context = "" + ccall((:__hsl_mc66_single_MOD_monet_print_message, libhsl_mc66), + Cvoid, + (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{UInt8}), + info , control.lp, control.wp, context ) +end + +function mc66d(m, n, nz, irn, jcn, nblocks, control, seed, row_order, info, rowptr, column_order, colptr, netcut, rowdiff, kblocks) + ccall((:__hsl_mc66_double_MOD_monet, libhsl_mc66), + Cvoid, + (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32}, Ref{Int32}, Ref{mc66_control{Float64}}, Ref{fa14_seed}, Ptr{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32} , Ptr{Int32}, Ref{Int32}, Ref{Float64}, Ref{Int32}), + m , n , nz , irn , jcn , nblocks , control , seed , row_order , info , rowptr , column_order, colptr , netcut , rowdiff , kblocks ) +end + +function mc66d_print_message(info, control) + context = "" + ccall((:__hsl_mc66_double_MOD_monet_print_message, libhsl_mc66), + Cvoid, + (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{UInt8}), + info , control.lp, control.wp, context ) +end + +""" + rp, cp = mc66(A::SparseMatrixCSC, k::Integer) + +Order an unsymmetric matrix A into singly bordered blocked diagonal (SBBD) form. + + A[rp, cp] = [ A₁₁ S₁] + [ A₂₂ S₂] + [ • • ] + [ • • ] + [ • • ] + [ Aₖₖ Sₖ] +""" +function mc66(A::SparseMatrixCSC, nblocks::Integer) + m, n = size(A) + nblocks ≤ min(m,n) || error("1 ≤ nblocks ≤ min(m,n)") + m = Ref{Int32}(m) + n = Ref{Int32}(n) + nz = Ref{Int32}(nnz(A)) + irn, jcn, val = findnz(A) + irn = convert(Vector{Cint}, irn) + jcn = convert(Vector{Cint}, jcn) + nblocks = Ref{Int32}(nblocks) + control = mc66_control{Float64}() + seed = fa14_seed() + row_order = zeros(Cint, m[]) + info = Ref{Cint}() + rowptr = zeros(Cint, nblocks[]+1) + column_order = zeros(Cint, n[]) + colptr = zeros(Cint, nblocks[]+1) + netcut = Ref{Cint}() + rowdiff = Ref{Float64}() + kblocks = Ref{Cint}() + mc66d(m, n, nz, irn, jcn, nblocks, control, seed, row_order, info, rowptr, column_order, colptr, netcut, rowdiff, kblocks) + mc66d_print_message(info, control) + return row_order, column_order +end diff --git a/test/test_hsl_mc66.jl b/test/test_hsl_mc66.jl new file mode 100644 index 0000000..bef76f6 --- /dev/null +++ b/test/test_hsl_mc66.jl @@ -0,0 +1,16 @@ +@testset "hsl_mc66" begin + A = [1 1 0 0; + 1 0 1 1; + 1 1 0 0; + 1 0 1 1] + + B = [1 0 0 1; + 1 0 0 1; + 0 1 1 1; + 0 1 1 1] + C = SparseMatrixCSC(A) + rp, cp = mc66(C, 2) + @test isperm(rp) + @test isperm(cp) + @test A[rp,cp] == B +end From 87c4c94d79e6411fcb3471add86f9903309d085f Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 1 May 2023 18:45:25 -0400 Subject: [PATCH 2/2] Update the interface of hsl_mc66 --- src/HSL.jl | 1 + src/hsl_mc66.jl | 30 ++++++++++++------------------ test/runtests.jl | 1 + 3 files changed, 14 insertions(+), 18 deletions(-) diff --git a/src/HSL.jl b/src/HSL.jl index 399192a..475c2c5 100644 --- a/src/HSL.jl +++ b/src/HSL.jl @@ -36,6 +36,7 @@ include("wrappers.jl") # Interfaces include("hsl_ma57.jl") include("hsl_ma97.jl") +include("hsl_mc66.jl") include("kb07.jl") include("mc21.jl") include("mc77.jl") diff --git a/src/hsl_mc66.jl b/src/hsl_mc66.jl index ad9c6a8..9d60d3a 100644 --- a/src/hsl_mc66.jl +++ b/src/hsl_mc66.jl @@ -29,35 +29,29 @@ function mc66_control{T}(; lp=6, wp=6, mp=6, print_level=-1, kl_aggressive=-1, m return control end -# Thanks to `nm -D libhsl_mc66.so`! +# Thanks to `nm -D libhsl.so`! function mc66s(m, n, nz, irn, jcn, nblocks, control, seed, row_order, info, rowptr, column_order, colptr, netcut, rowdiff, kblocks) - ccall((:__hsl_mc66_simple_MOD_monet, libhsl_mc66), - Cvoid, - (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32}, Ref{Int32}, Ref{mc66_control{Float32}}, Ref{fa14_seed}, Ptr{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32} , Ptr{Int32}, Ref{Int32}, Ref{Float32}, Ref{Int32}), - m , n , nz , irn , jcn , nblocks , control , seed , row_order , info , rowptr , column_order, colptr , netcut , rowdiff , kblocks ) + @ccall libhsl.__hsl_mc66_simple_MOD_monet(m::Ref{Int32}, n::Ref{Int32}, nz::Ptr{Int32}, irn::Ptr{Int32}, jcn::Ptr{Int32}, + nblocks::Ref{Int32}, control::Ref{mc66_control{Float32}}, seed::Ref{fa14_seed}, + row_order::Ptr{Int32}, info::Ref{Int32}, rowptr::Ptr{Int32}, column_order::Ptr{Int32}, + colptr::Ptr{Int32}, netcut::Ref{Int32}, rowdiff::Ref{Float32}, kblocks::Ref{Int32})::Cvoid end function mc66s_print_message(info, control) context = "" - ccall((:__hsl_mc66_single_MOD_monet_print_message, libhsl_mc66), - Cvoid, - (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{UInt8}), - info , control.lp, control.wp, context ) + @ccall libhsl.__hsl_mc66_single_MOD_monet_print_message(info::Ref{Int32}, control.lp::Ref{Int32}, control.wp::Ref{Int32}, context::Ptr{UInt8})::Cvoid end function mc66d(m, n, nz, irn, jcn, nblocks, control, seed, row_order, info, rowptr, column_order, colptr, netcut, rowdiff, kblocks) - ccall((:__hsl_mc66_double_MOD_monet, libhsl_mc66), - Cvoid, - (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32}, Ref{Int32}, Ref{mc66_control{Float64}}, Ref{fa14_seed}, Ptr{Int32}, Ref{Int32}, Ptr{Int32}, Ptr{Int32} , Ptr{Int32}, Ref{Int32}, Ref{Float64}, Ref{Int32}), - m , n , nz , irn , jcn , nblocks , control , seed , row_order , info , rowptr , column_order, colptr , netcut , rowdiff , kblocks ) + @ccall libhsl.__hsl_mc66_double_MOD_monet(m::Ref{Int32}, n::Ref{Int32}, nz::Ptr{Int32}, irn::Ptr{Int32}, jcn::Ptr{Int32}, + nblocks::Ref{Int32}, control::Ref{mc66_control{Float64}}, seed::Ref{fa14_seed}, + row_order::Ptr{Int32}, info::Ref{Int32}, rowptr::Ptr{Int32}, column_order::Ptr{Int32}, + colptr::Ptr{Int32}, netcut::Ref{Int32}, rowdiff::Ref{Float64}, kblocks::Ref{Int32})::Cvoid end function mc66d_print_message(info, control) context = "" - ccall((:__hsl_mc66_double_MOD_monet_print_message, libhsl_mc66), - Cvoid, - (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ptr{UInt8}), - info , control.lp, control.wp, context ) + @ccall libhsl.__hsl_mc66_double_MOD_monet_print_message(info::Ref{Int32}, control.lp::Ref{Int32}, control.wp::Ref{Int32}, context::Ptr{UInt8})::Cvoid end """ @@ -93,6 +87,6 @@ function mc66(A::SparseMatrixCSC, nblocks::Integer) rowdiff = Ref{Float64}() kblocks = Ref{Cint}() mc66d(m, n, nz, irn, jcn, nblocks, control, seed, row_order, info, rowptr, column_order, colptr, netcut, rowdiff, kblocks) - mc66d_print_message(info, control) + # mc66d_print_message(info, control) return row_order, column_order end diff --git a/test/runtests.jl b/test/runtests.jl index 48f7c18..30459e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ if JULIAHSL_isfunctional() include("test_hsl_ma57.jl") # include("test_hsl_ma57_patch.jl") include("test_hsl_ma97.jl") + include("test_hsl_mc66.jl") include("test_kb07.jl") include("test_mc21.jl") include("test_mc77.jl")