From 79f98b4e14eb63968fc08f56f20f725b59b7fb70 Mon Sep 17 00:00:00 2001 From: William R Sweeney Date: Fri, 31 May 2019 15:46:35 -0400 Subject: [PATCH] removed gc_haven, added tests --- src/MUMPS3.jl | 15 ++- src/convenience_wrappers.jl | 70 ++++++------- src/icntl_alibis.jl | 8 +- src/interface.jl | 193 +++++++++++++++++------------------- src/mumps_struc.jl | 58 ++--------- src/printing.jl | 92 ++++++----------- test/advanced_solve.jl | 19 ++++ test/basic_solve.jl | 17 ++++ test/runtests.jl | 6 +- test/schur_complement.jl | 21 ++++ test/simple_test.jl | 10 -- 11 files changed, 242 insertions(+), 267 deletions(-) create mode 100644 test/advanced_solve.jl create mode 100644 test/basic_solve.jl create mode 100644 test/schur_complement.jl delete mode 100644 test/simple_test.jl diff --git a/src/MUMPS3.jl b/src/MUMPS3.jl index 12b1b97..08a6308 100644 --- a/src/MUMPS3.jl +++ b/src/MUMPS3.jl @@ -38,11 +38,20 @@ SparseArrays function __init__() if haskey(ENV,"MUMPS_PREFIX") - global MUMPS_LIB = joinpath(ENV["MUMPS_PREFIX"],"lib/libmumps_simple.dylib") + global MUMPS_LIB_S = joinpath(ENV["MUMPS_PREFIX"],"lib/libsmumps.dylib") + global MUMPS_LIB_D = joinpath(ENV["MUMPS_PREFIX"],"lib/libdmumps.dylib") + global MUMPS_LIB_C = joinpath(ENV["MUMPS_PREFIX"],"lib/libcmumps.dylib") + global MUMPS_LIB_Z = joinpath(ENV["MUMPS_PREFIX"],"lib/libzmumps.dylib") else - global MUMPS_LIB = "/usr/local/opt/brewsci-mumps/lib/libmumps_simple.dylib" + global MUMPS_LIB_S = "/usr/local/opt/brewsci-mumps/lib/libsmumps.dylib" + global MUMPS_LIB_D = "/usr/local/opt/brewsci-mumps/lib/libdmumps.dylib" + global MUMPS_LIB_C = "/usr/local/opt/brewsci-mumps/lib/libcmumps.dylib" + global MUMPS_LIB_Z = "/usr/local/opt/brewsci-mumps/lib/libzmumps.dylib" end - global LIB = dlopen(MUMPS_LIB) + global LIB_S = dlopen(MUMPS_LIB_S) + global LIB_D = dlopen(MUMPS_LIB_D) + global LIB_C = dlopen(MUMPS_LIB_C) + global LIB_Z = dlopen(MUMPS_LIB_Z) end include("mumps_types.jl") diff --git a/src/convenience_wrappers.jl b/src/convenience_wrappers.jl index 5e76211..4cb115e 100644 --- a/src/convenience_wrappers.jl +++ b/src/convenience_wrappers.jl @@ -81,14 +81,14 @@ See also: [`mumps_solve`](@ref), [`get_sol!`](@ref), [`get_sol`](@ref) function mumps_solve!(mumps::Mumps) @assert has_matrix(mumps) "matrix not yet provided to mumps object" @assert has_rhs(mumps) "rhs not yet provided to mumps object" - if mumps.mumpsc.job ∈ [2,4] # if already factored, just solve - mumps.mumpsc.job = 3 - elseif mumps.mumpsc.job ∈ [1] # if analyzed only, factorize and solve - mumps.mumpsc.job=5 - elseif mumps.mumpsc.job ∈ [3,5,6] # is solved already, retrieve solution + if mumps.job ∈ [2,4] # if already factored, just solve + mumps.job = 3 + elseif mumps.job ∈ [1] # if analyzed only, factorize and solve + mumps.job=5 + elseif mumps.job ∈ [3,5,6] # is solved already, retrieve solution return nothing else # else analyze, factor, solve - mumps.mumpsc.job=6 + mumps.job=6 end invoke_mumps!(mumps) return nothing @@ -106,13 +106,13 @@ function mumps_solve!(x::AbstractArray,A::AbstractArray,rhs::AbstractArray; kwar end function mumps_solve!(x::AbstractArray,mumps::Mumps,rhs::AbstractArray) provide_rhs!(mumps,rhs) - if mumps.mumpsc.job ∈ [2,4] # if already factored, just solve - elseif mumps.mumpsc.job ∈ [1] # if analyzed only, factorize and solve - mumps.mumpsc.job=5 - elseif mumps.mumpsc.job ∈ [3,5,6] # is solved already, reset to solve only - mumps.mumpsc.job=2 + if mumps.job ∈ [2,4] # if already factored, just solve + elseif mumps.job ∈ [1] # if analyzed only, factorize and solve + mumps.job=5 + elseif mumps.job ∈ [3,5,6] # is solved already, reset to solve only + mumps.job=2 else # else analyze, factor, solve - mumps.mumpsc.job=6 + mumps.job=6 end mumps_solve!(x,mumps) end @@ -128,7 +128,7 @@ If only input is `mumps` must also have been provided `y`. See also: [`mumps_solve!`](@ref) """ function mumps_solve(mumps::Mumps) - if mumps.mumpsc.job ∈ [3,5,6] + if mumps.job ∈ [3,5,6] x = get_sol(mumps) else x = get_rhs(mumps) @@ -160,13 +160,13 @@ See also: [`mumps_factorize`](@ref) """ function mumps_factorize!(mumps::Mumps) @assert has_matrix(mumps) "matrix not yet provided to mumps object" - if mumps.mumpsc.job ∈ [2,3,4,5,6] # already factored + if mumps.job ∈ [2,3,4,5,6] # already factored @warn "already factored" return nothing - elseif mumps.mumpsc.job ∈ [1] # if analyzed only, factorize - mumps.mumpsc.job=2 + elseif mumps.job ∈ [1] # if analyzed only, factorize + mumps.job=2 else # else analyze, factor - mumps.mumpsc.job=4 + mumps.job=4 end invoke_mumps!(mumps) end @@ -205,12 +205,12 @@ See also: [`mumps_det`](@ref) """ function mumps_det!(mumps::Mumps; discard=true) @assert has_matrix(mumps) "matrix not yet provided to mumps object" - if has_det(mumps) && mumps.mumpsc.job>1 + if has_det(mumps) && mumps.job>1 return nothing end set_icntl!(mumps,31,Int(discard)) set_icntl!(mumps,33,1) - mumps.mumpsc.job>0 ? mumps.mumpsc.job=1 : nothing + mumps.job>0 ? mumps.job=1 : nothing mumps_factorize!(mumps) end """ @@ -242,14 +242,14 @@ See also: [`mumps_schur_complement`](@ref), [`get_schur_complement!`](@ref), [`g function mumps_schur_complement!(mumps::Mumps, schur_inds::AbstractArray{Int,1}) @assert has_matrix(mumps) "matrix not yet provided to mumps object" set_schur_centralized_by_column!(mumps, schur_inds) - if mumps.mumpsc.job ∈ [1] # if analyzed only, factorize - mumps.mumpsc.job=2 + if mumps.job ∈ [1] # if analyzed only, factorize + mumps.job=2 else # else analyze, factor - mumps.mumpsc.job=4 + mumps.job=4 end invoke_mumps!(mumps) end -mumps_schur_complement!(mumps::Mumps, x::SparseMatrixCSC) = mumps_schur_complement!(mumps,unique!(sort!(x.rowval))) +mumps_schur_complement!(mumps::Mumps, x::SparseMatrixCSC) = mumps_schur_complement!(mumps,unique!(sort(x.rowval))) mumps_schur_complement!(mumps::Mumps, x::SparseVector) = mumps_schur_complement!(mumps,x.nzind) """ mumps_schur_complement(A,schur_inds) -> S @@ -287,12 +287,12 @@ function mumps_select_inv!(x::AbstractSparseArray,mumps::Mumps) set_icntl!(mumps,30,1) set_icntl!(mumps,20,3) provide_rhs!(mumps,x) - if mumps.mumpsc.job ∈ [2,4] # if already factored, just solve - mumps.mumpsc.job = 3 - elseif mumps.mumpsc.job ∈ [1] # if analyzed only, factorize and solve - mumps.mumpsc.job=5 + if mumps.job ∈ [2,4] # if already factored, just solve + mumps.job = 3 + elseif mumps.job ∈ [1] # if analyzed only, factorize and solve + mumps.job=5 else # else analyze, factor, solve - mumps.mumpsc.job=6 + mumps.job=6 end invoke_mumps!(mumps) get_rhs!(x,mumps) @@ -338,7 +338,7 @@ See also: [`finalize`](@ref) function initialize!(mumps::Mumps) suppress_display!(mumps) mumps.finalized=false - mumps.mumpsc.job=-1 + mumps.job=-1 invoke_mumps!(mumps) end @@ -357,14 +357,14 @@ function finalize!(mumps::Mumps) end function finalize_unsafe!(mumps::Mumps) mumps.finalized=true - mumps.mumpsc.job=-2 + mumps.job=-2 invoke_mumps_unsafe!(mumps) end function Base.:\(mumps::Mumps,y) suppress_display!(mumps) - if mumps.mumpsc.job < 2 + if mumps.job < 2 mumps_factorize!(mumps) end x = mumps_solve(mumps,y) @@ -372,24 +372,24 @@ function Base.:\(mumps::Mumps,y) end function LinearAlgebra.ldiv!(mumps::Mumps,y) suppress_display!(mumps) - if mumps.mumpsc.job < 2 + if mumps.job < 2 mumps_factorize!(mumps) else - mumps.mumpsc.job=2 + mumps.job=2 end provide_rhs!(mumps,y) mumps_solve!(y,mumps) end function LinearAlgebra.ldiv!(x,mumps::Mumps,y) suppress_display!(mumps) - if mumps.mumpsc.job < 2 + if mumps.job < 2 mumps_factorize!(mumps) end mumps_solve!(x,mumps,y) end function LinearAlgebra.inv(mumps::Mumps) suppress_display!(mumps) - y = sparse(1:mumps.mumpsc.n,1:mumps.mumpsc.n,1:mumps.mumpsc.n,mumps.mumpsc.n,mumps.mumpsc.n) + y = sparse(1:mumps.n,1:mumps.n,1:mumps.n,mumps.n,mumps.n) mumps_select_inv!(y,mumps) finalize!(mumps) return y diff --git a/src/icntl_alibis.jl b/src/icntl_alibis.jl index b5e3432..a261b31 100644 --- a/src/icntl_alibis.jl +++ b/src/icntl_alibis.jl @@ -14,7 +14,7 @@ const ICNTL_DEFAULT = (6, 0, 6, 2, 0, 7, 7, 77, 1, 0, 0, 1, 0, 20, 0, 0, 0, 0, 0 reset ICNTL to its default """ function default_icntl!(mumps::Mumps) - mumps.mumpsc.icntl = ICNTL_DEFAULT + mumps.icntl = ICNTL_DEFAULT return nothing end @@ -24,16 +24,16 @@ set_info_stream!(mumps::Mumps,i) = set_icntl!(mumps,3,i; displaylevel=0) set_print_level!(mumps::Mumps,i) = set_icntl!(mumps,4,i; displaylevel=0) suppress_printing!(mumps::Mumps) = set_print_level!(mumps,1) -toggle_printing!(mumps::Mumps) = set_print_level!(mumps,mod1(mumps.mumpsc.icntl[4]+1,2)) +toggle_printing!(mumps::Mumps) = set_print_level!(mumps,mod1(mumps.icntl[4]+1,2)) suppress_display! = suppress_printing! toggle_display! = toggle_printing! sparse_matrix!(mumps::Mumps) = set_icntl!(mumps,5,0; displaylevel=0) dense_matrix!(mumps::Mumps) = set_icntl!(mumps,5,1; displaylevel=0) -LinearAlgebra.transpose!(mumps::Mumps) = set_icntl!(mumps,9,mod(mumps.mumpsc.icntl[9]+1,2); displaylevel=0) +LinearAlgebra.transpose!(mumps::Mumps) = set_icntl!(mumps,9,mod(mumps.icntl[9]+1,2); displaylevel=0) sparse_rhs!(mumps::Mumps) = set_icntl!(mumps, 20, 1; displaylevel=0) dense_rhs!(mumps::Mumps) = set_icntl!(mumps, 20, 0; displaylevel=0) -toggle_null_pivot!(mumps::Mumps) = set_icntl!(mumps,24,mod(mumps.mumpsc.icntl[24]+1,2); displaylevel=0) +toggle_null_pivot!(mumps::Mumps) = set_icntl!(mumps,24,mod(mumps.icntl[24]+1,2); displaylevel=0) diff --git a/src/interface.jl b/src/interface.jl index 5aad4a4..4dfc732 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -26,20 +26,23 @@ initialized. See also: [`invoke_mumps!`](@ref) """ -invoke_mumps_unsafe!(mumps::Mumps) = invoke_mumps_unsafe!(mumps.mumpsc) -@inline function invoke_mumps_unsafe!(mumpsc::MumpsC{TC,TR}) where {TC,TR} +@inline function invoke_mumps_unsafe!(mumps::Mumps{TC,TR}) where {TC,TR} @assert MPI.Initialized() "must call MPI.Init() exactly once before calling mumps" if TC==Float32 cfun = :smumps_c + LIB = LIB_S elseif TC==Float64 cfun = :dmumps_c + LIB = LIB_D elseif TC==ComplexF32 cfun = :cmumps_c + LIB = LIB_C elseif TC==ComplexF64 cfun = :zmumps_c + LIB = LIB_Z end sym = dlsym(LIB,cfun) - ccall(sym,Cvoid,(Ref{MumpsC{TC,TR}},), mumpsc) + ccall(sym,Cvoid,(Ref{Mumps{TC,TR}},), mumps) return nothing end """ @@ -60,7 +63,11 @@ function invoke_mumps!(mumps::Mumps) check_finalized(mumps) invoke_mumps_unsafe!(mumps) end -check_finalized(mumps::Mumps) = @assert !mumps.finalized "Mumps object already finalized" +function check_finalized(mumps::Mumps) + if mumps.finalized + throw(ErrorException("Mumps object already finalized")) + end +end """ @@ -70,10 +77,10 @@ Set the integer control parameters according to ICNTL[i]=val See also: [`display_icntl`](@ref) """ -function set_icntl!(mumps::Mumps,i::Int,val::Int; displaylevel=mumps.mumpsc.icntl[4]-1) - icntl = mumps.mumpsc.icntl - mumps.mumpsc.icntl = (icntl[1:i-1]...,convert(MUMPS_INT,val),icntl[i+1:end]...) - displaylevel>0 ? display_icntl(stdout,mumps.mumpsc.icntl,i,val) : nothing +function set_icntl!(mumps::Mumps,i::Int,val::Int; displaylevel=mumps.icntl[4]-1) + icntl = mumps.icntl + mumps.icntl = (icntl[1:i-1]...,convert(MUMPS_INT,val),icntl[i+1:end]...) + displaylevel>0 ? display_icntl(stdout,mumps.icntl,i,val) : nothing return nothing end @@ -85,10 +92,10 @@ Set the real/complex control parameters according to CNTL[i]=val See also: [`display_cntl`](@ref) """ -function set_cntl!(mumps::Mumps{TC,TR},i::Int,val::Float64; displaylevel=mumps.mumpsc.icntl[4]-1) where {TC,TR} - cntl = mumps.mumpsc.cntl - mumps.mumpsc.cntl = (cntl[1:i-1]...,convert(TR,val),cntl[i+1:end]...) - displaylevel>0 ? display_cntl(stdout,mumps.mumpsc.cntl,i,val) : nothing +function set_cntl!(mumps::Mumps{TC,TR},i::Int,val::Float64; displaylevel=mumps.icntl[4]-1) where {TC,TR} + cntl = mumps.cntl + mumps.cntl = (cntl[1:i-1]...,convert(TR,val),cntl[i+1:end]...) + displaylevel>0 ? display_cntl(stdout,mumps.cntl,i,val) : nothing return nothing end @@ -99,7 +106,7 @@ end Set the phase to `job`. See MUMPS manual for options. """ function set_job!(mumps::Mumps,i) - mumps.mumpsc.job=i + mumps.job=i return nothing end @@ -112,8 +119,8 @@ set name of directory in which to store out-of-core files. function set_save_dir!(mumps,dir::String) @assert length(dir≤255) "directory name has $(length(dir)) characters, must be ≤255" i = length(dir+1) - save_dir = mumps.mumpsc.save_dir - mumps.mumpsc.save_dir = (dir...,'\0',save_dir[i+2:end]...) + save_dir = mumps.save_dir + mumps.save_dir = (dir...,'\0',save_dir[i+2:end]...) return nothing end @@ -126,8 +133,8 @@ prefix for out-of-core files. function set_save_prefix!(mumps,prefix::String) @assert length(prefix) "prefix name has $(length(prefix)) characters, must be ≤255" i = length(prefix+1) - save_prefix = mumps.mumpsc.save_prefix - mumps.mumpsc.save_prefix = (prefix...,'\0',save_prefix[i+2:end]...) + save_prefix = mumps.save_prefix + mumps.save_prefix = (prefix...,'\0',save_prefix[i+2:end]...) return nothing end @@ -149,10 +156,10 @@ function provide_matrix!(mumps::Mumps{T},A::AbstractArray{TA}) where {T,TA} @warn "matrix with element type $TA will attempt be converted to Mumps type $T" end if is_matrix_assembled(mumps) - typeof(A)<:SparseMatrixCSC ? nothing : (@warn "matrix is dense, but ICNTL[5]=$(mumps.mumpsc.icntl[5]) indicates assembled. attempting to convert matrix to sparse") + typeof(A)<:SparseMatrixCSC ? nothing : (@warn "matrix is dense, but ICNTL[5]=$(mumps.icntl[5]) indicates assembled. attempting to convert matrix to sparse") _provide_matrix_assembled!(mumps,convert(SparseMatrixCSC,A)) else - !(typeof(A)<:SparseMatrixCSC) ? nothing : (@warn "matrix is sparse, but ICNTL[5]=$(mumps.mumpsc.icntl[5]) indicates elemental. attempting to convert matrix to dense") + !(typeof(A)<:SparseMatrixCSC) ? nothing : (@warn "matrix is sparse, but ICNTL[5]=$(mumps.icntl[5]) indicates elemental. attempting to convert matrix to dense") _provide_matrix_elemental!(mumps,convert(Matrix,A)) end end @@ -164,17 +171,16 @@ function _provide_matrix_assembled!(mumps::Mumps,A::SparseMatrixCSC{T}) where T end end function _provide_matrix_assembled_centralized!(mumps::Mumps{T},A::SparseMatrixCSC) where T - mumpsc, gc_haven = mumps.mumpsc, mumps.gc_haven if is_symmetric(mumps) I,J,V = findnz(triu(A)) else I,J,V = findnz(A) end irn, jcn, a = convert.((Array{MUMPS_INT},Array{MUMPS_INT},Array{T}),(I,J,V)) - gc_haven.irn, gc_haven.jcn, gc_haven.a = irn, jcn, a - mumpsc.irn, mumpsc.jcn, mumpsc.a = pointer.((irn,jcn,a)) - mumpsc.n = A.n - mumpsc.nnz = length(V) + mumps.irn, mumps.jcn, mumps.a = pointer.((irn,jcn,a)) + mumps.n = A.n + mumps.nnz = length(V) + append!(mumps.gc_haven,[Ref(irn),Ref(jcn),Ref(a)]) return nothing end function _provide_matrix_assembled_distributed!(mumps::Mumps{T},A::SparseMatrixCSC) where T @@ -182,20 +188,18 @@ function _provide_matrix_assembled_distributed!(mumps::Mumps{T},A::SparseMatrixC return nothing end function _provide_matrix_elemental!(mumps::Mumps{T},A::Array) where T - mumpsc, gc_haven = mumps.mumpsc, mumps.gc_haven - mumpsc.n = size(A,1) - mumpsc.nelt = 1 - eltptr = convert.(MUMPS_INT,[1,mumpsc.n+1]) - eltvar = convert.(MUMPS_INT,collect(1:mumpsc.n)) - gc_haven.eltptr, gc_haven.eltvar = eltptr, eltvar - mumpsc.eltptr, mumpsc.eltvar = pointer(eltptr), pointer(eltvar) + mumps.n = size(A,1) + mumps.nelt = 1 + eltptr = convert.(MUMPS_INT,[1,mumps.n+1]) + eltvar = convert.(MUMPS_INT,collect(1:mumps.n)) + mumps.eltptr, mumps.eltvar = pointer.(eltptr, eltvar) if is_symmetric(mumps) - a_elt = convert.(T,[A[i,j] for i ∈ 1:mumpsc.n for j ∈ 1:i]) + a_elt = convert.(T,[A[i,j] for i ∈ 1:mumps.n for j ∈ 1:i]) else a_elt = convert.(T,A[:]) end - gc_haven.a_elt = a_elt - mumpsc.a_elt = pointer(a_elt) + mumps.a_elt = pointer(a_elt) + append!(mumps.gc_haven,[Ref(eltptr),Ref(eltvar),Ref(a_elt)]) return nothing end @@ -220,38 +224,33 @@ function provide_rhs!(mumps::Mumps,rhs::AbstractMatrix) end end function provide_rhs_sparse!(mumps::Mumps{T},rhs::AbstractMatrix) where T - gc_haven, mumpsc = mumps.gc_haven, mumps.mumpsc rhs = convert(SparseMatrixCSC,rhs) - mumpsc.nz_rhs = length(rhs.nzval) - mumpsc.nrhs = size(rhs,2) + mumps.nz_rhs = length(rhs.nzval) + mumps.nrhs = size(rhs,2) rhs_sparse = convert.(T,rhs.nzval) irhs_sparse = convert.(MUMPS_INT,rhs.rowval) irhs_ptr = convert.(MUMPS_INT,rhs.colptr) - gc_haven.rhs_sparse = rhs_sparse - gc_haven.irhs_sparse = irhs_sparse - gc_haven.irhs_ptr = irhs_ptr - - mumpsc.rhs_sparse = pointer(rhs_sparse) - mumpsc.irhs_sparse = pointer(irhs_sparse) - mumpsc.irhs_ptr = pointer(irhs_ptr) + mumps.rhs_sparse = pointer(rhs_sparse) + mumps.irhs_sparse = pointer(irhs_sparse) + mumps.irhs_ptr = pointer(irhs_ptr) if is_sol_central(mumps) y = fill(convert(T,NaN),prod(size(rhs))) - gc_haven.rhs = y - mumpsc.rhs = pointer(y) - mumpsc.lrhs = size(rhs,1) + mumps.rhs = pointer(y) + mumps.lrhs = size(rhs,1) + push!(mumps.gc_haven,Ref(y)) end + append!(mumps.gc_haven,[Ref(rhs_sparse),Ref(irhs_sparse),Ref(irhs_ptr)]) return nothing end function provide_rhs_dense!(mumps::Mumps{T},rhs::AbstractMatrix) where T - gc_haven, mumpsc = mumps.gc_haven, mumps.mumpsc y = convert(Matrix{T},rhs)[:] - gc_haven.rhs = y - mumpsc.rhs = pointer(y) - mumpsc.lrhs = size(rhs,1) - mumpsc.nrhs = size(rhs,2) + mumps.rhs = pointer(y) + mumps.lrhs = size(rhs,1) + mumps.nrhs = size(rhs,2) + push!(mumps.gc_haven,Ref(y)) return nothing end provide_rhs!(mumps::Mumps,rhs::AbstractVector) = provide_rhs!(mumps,repeat(rhs,1,1)) @@ -276,18 +275,18 @@ end function get_rhs_unsafe!(x::SparseMatrixCSC,mumps::Mumps) @assert !is_rhs_dense(mumps) "rhs is dense, target is sparse. try with dense target" for i ∈ LinearIndices(x.colptr) - x.colptr[i] = unsafe_load(mumps.mumpsc.irhs_ptr,i) + x.colptr[i] = unsafe_load(mumps.irhs_ptr,i) end for i ∈ LinearIndices(x.rowval) - x.rowval[i] = unsafe_load(mumps.mumpsc.irhs_sparse,i) - x.nzval[i] = unsafe_load(mumps.mumpsc.rhs_sparse,i) + x.rowval[i] = unsafe_load(mumps.irhs_sparse,i) + x.nzval[i] = unsafe_load(mumps.rhs_sparse,i) end return nothing end function get_rhs_unsafe!(x::Union{SubArray,Array},mumps::Mumps) @assert is_rhs_dense(mumps) "rhs is sparse, target is dense. try with sparse target" for i ∈ LinearIndices(x) - x[i] = unsafe_load(mumps.mumpsc.rhs,i) + x[i] = unsafe_load(mumps.s,i) end return nothing end @@ -299,15 +298,15 @@ Retrieve right hand side from `mumps` See also: [`get_rhs!`](@ref), [`get_sol!`](@ref), [`get_sol`](@ref) """ function get_rhs(mumps::Mumps{T}) where T - n = mumps.mumpsc.nrhs + n = mumps.nrhs if !is_rhs_dense(mumps) - m = mumps.mumpsc.n - colptr = ones(MUMPS_INT,mumps.mumpsc.nrhs+1) - rowval = ones(MUMPS_INT,mumps.mumpsc.nz_rhs) - nzval = Array{T}(undef,mumps.mumpsc.nz_rhs) + m = mumps.n + colptr = ones(MUMPS_INT,mumps.nrhs+1) + rowval = ones(MUMPS_INT,mumps.nz_rhs) + nzval = Array{T}(undef,mumps.nz_rhs) x = SparseMatrixCSC(m,n,colptr,rowval,nzval) else - m = mumps.mumpsc.lrhs + m = mumps.lrhs x = Array{T}(undef,m,n) end get_rhs!(x,mumps) @@ -324,7 +323,7 @@ See also: [`get_rhs!`](@ref), [`get_rhs`](@ref), [`get_sol`](@ref) """ function get_sol!(x::Union{SubArray,Array},mumps::Mumps) check_finalized(mumps) - if mumps.mumpsc.job ∉ [3,5,6] + if mumps.job ∉ [3,5,6] @warn "mumps has not passed through a solution phase" end if has_rhs(mumps) @@ -336,7 +335,7 @@ function get_sol!(x::Union{SubArray,Array},mumps::Mumps) end function get_sol_unsafe!(x::Union{SubArray,Array},mumps::Mumps) for i ∈ LinearIndices(x) - x[i] = unsafe_load(mumps.mumpsc.rhs,i) + x[i] = unsafe_load(mumps.rhs,i) end return nothing end @@ -348,10 +347,10 @@ Retrieve solution from `mumps` See also: [`get_rhs!`](@ref), [`get_rhs`](@ref), [`get_sol!`](@ref) """ function get_sol(mumps::Mumps{T}) where T - if mumps.mumpsc.job ∉ [3,5,6] + if mumps.job ∉ [3,5,6] @warn "mumps has not passed through a solution phase" end - x = Array{T}(undef,mumps.mumpsc.lrhs,mumps.mumpsc.nrhs) + x = Array{T}(undef,mumps.lrhs,mumps.nrhs) get_sol!(x,mumps) return x end @@ -370,7 +369,7 @@ function get_schur_complement!(S,mumps::Mumps) end function get_schur_complement_unsafe!(S,mumps::Mumps) for i ∈ LinearIndices(S) - S[i] = unsafe_load(mumps.mumpsc.schur,i) + S[i] = unsafe_load(mumps.schur,i) end return nothing end @@ -382,7 +381,7 @@ Retrieve Schur complement matrix `S` from `mumps` See also: [`get_schur_complement!`](@ref), [`mumps_schur!`](@ref), [`mumps_schur`](@ref) """ function get_schur_complement(mumps::Mumps{T}) where T - S = Array{T}(undef,mumps.mumpsc.size_schur,mumps.mumpsc.size_schur) + S = Array{T}(undef,mumps.size_schur,mumps.size_schur) get_schur_complement!(S,mumps) return S end @@ -397,20 +396,18 @@ method suggested in the MUMPS manual See also: [`mumps_schur!`](@ref), [`mumps_schur`](@ref) """ function set_schur_centralized_by_column!(mumps::Mumps{T},schur_inds::AbstractArray{Int}) where T - gc_haven, mumpsc = mumps.gc_haven, mumps.mumpsc - mumpsc.size_schur = length(schur_inds) + mumps.size_schur = length(schur_inds) listvar_schur = convert.(MUMPS_INT,schur_inds) - gc_haven.listvar_schur = listvar_schur - mumpsc.listvar_schur = pointer(listvar_schur) - mumpsc.nprow = 1 - mumpsc.npcol = 1 - mumpsc.mblock = 100 - mumpsc.nblock = 100 - mumpsc.schur_lld = mumpsc.size_schur - schur = Array{T}(undef,mumps.mumpsc.size_schur^2) - gc_haven.schur = schur - mumpsc.schur = pointer(schur) + mumps.listvar_schur = pointer(listvar_schur) + mumps.nprow = 1 + mumps.npcol = 1 + mumps.mblock = 100 + mumps.nblock = 100 + mumps.schur_lld = mumps.size_schur + schur = Array{T}(undef,mumps.size_schur^2) + mumps.schur = pointer(schur) set_icntl!(mumps,19,3) + append!(mumps.gc_haven,[Ref(listvar_schur),Ref(schur)]) end @@ -420,9 +417,9 @@ function LinearAlgebra.det(mumps::Mumps{T}) where T d = T<:Complex ? complex(NaN,NaN) : NaN else if T<:Complex - d = complex(mumps.mumpsc.rinfog[12],mumps.mumpsc.rinfog[13])*2^mumps.mumpsc.infog[34] + d = complex(mumps.rinfog[12],mumps.rinfog[13])*2^mumps.infog[34] else - d = mumps.mumpsc.rinfog[12]*2^mumps.mumpsc.infog[34] + d = mumps.rinfog[12]*2^mumps.infog[34] end end return convert(T,d) @@ -432,37 +429,25 @@ end #################################################################### ### auxilliary functions #################################################################### -is_matrix_assembled(mumps::Mumps) = is_matrix_assembled(mumps.mumpsc) -is_matrix_assembled(mumpsc::MumpsC) = !(mumpsc.icntl[5] ∈ [1]) +is_matrix_assembled(mumps::Mumps) = !(mumps.icntl[5] ∈ [1]) -is_matrix_distributed(mumps::Mumps) = is_matrix_distributed(mumps.mumpsc) -is_matrix_distributed(mumpsc::MumpsC) = mumpsc.icntl[18] ∈ [1,2,3] +is_matrix_distributed(mumps::Mumps) = mumps.icntl[18] ∈ [1,2,3] -is_rhs_dense(mumps::Mumps) = is_rhs_dense(mumps.mumpsc) -is_rhs_dense(mumpsc::MumpsC) = mumpsc.icntl[20] ∉ [1,2,3] +is_rhs_dense(mumps::Mumps) = mumps.icntl[20] ∉ [1,2,3] -is_sol_central(mumps::Mumps) = is_sol_central(mumps.mumpsc) -is_sol_central(mumpsc::MumpsC) = mumpsc.icntl[21] ∉ [1] +is_sol_central(mumps::Mumps) = mumps.icntl[21] ∉ [1] -has_det(mumps::Mumps) = has_det(mumps.mumpsc) -has_det(mumpsc::MumpsC) = mumpsc.icntl[33] ∉ [0] +has_det(mumps::Mumps) = mumps.icntl[33] ∉ [0] -is_symmetric(mumps::Mumps) = is_symmetric(mumps.mumpsc) -is_symmetric(mumpsc::MumpsC) = mumpsc.sym ∈ [1,2] +is_symmetric(mumps::Mumps) = mumps.sym ∈ [1,2] -is_posdef(mumps::Mumps) = is_posdef(mumps.mumpsc) -is_posdef(mumpsc::MumpsC) = mumpsc.sym ∈ [1] +is_posdef(mumps::Mumps) = mumps.sym ∈ [1] -has_matrix(mumps::Mumps) = has_matrix(mumps.mumpsc) -has_matrix(mumpsc::MumpsC) = mumpsc.n > 0 +has_matrix(mumps::Mumps) = mumps.n > 0 -has_rhs(mumps::Mumps) = has_rhs(mumps.mumpsc) -has_rhs(mumpsc::MumpsC) = mumpsc.nrhs*mumpsc.lrhs>0 || mumpsc.nz_rhs>0 +has_rhs(mumps::Mumps) = mumps.nrhs*mumps.lrhs>0 || mumps.nz_rhs>0 -has_schur(mumps::Mumps) = has_schur(mumps.mumpsc) -has_schur(mumpsc::MumpsC) = mumpsc.size_schur > 0 +has_schur(mumps::Mumps) = mumps.size_schur > 0 LinearAlgebra.issymmetric(mumps::Mumps) = is_symmetric(mumps) -LinearAlgebra.issymmetric(mumpsc::MumpsC) = is_symmetric(mumpsc) LinearAlgebra.isposdef(mumps::Mumps) = is_posdef(mumps) -LinearAlgebra.isposdef(mumpsc::MumpsC) = is_posdef(mumpsc) diff --git a/src/mumps_struc.jl b/src/mumps_struc.jl index a6d5bf6..1a5a4ae 100644 --- a/src/mumps_struc.jl +++ b/src/mumps_struc.jl @@ -9,7 +9,7 @@ const MUMPS_VERSION = "5.2.0" const MUMPS_VERSION_MAX_LEN = 30 # mirror of structre in [sdcz]mumps_c.h -mutable struct MumpsC{TC,TR} +mutable struct Mumps{TC,TR} sym::MUMPS_INT # MANDATORY 0 for unsymmetric, 1 for symmetric and posdef, 2 for general symmetric. All others treated as 0 par::MUMPS_INT # MANDATORY 0 host not involved in parallel factorization and solve, 1 host is involved job::MUMPS_INT # MANDATORY -1 initializes package, must come first, -2 terminates, 1 analysis, 2 factorization, 3 solve, 4=1&2, 5=2&3, 6=1&2&3 @@ -106,57 +106,15 @@ mutable struct MumpsC{TC,TR} metis_options::NTuple{40,MUMPS_INT} - MumpsC{T}(sym::Int,par::Int,comm) where T = new{T,real(T)}(sym,par,-1,comm) -end - - -# structure to store chucnks of data to keep it safe from Julia gc -mutable struct GC_haven{TC,TR} - irn::Vector{MUMPS_INT} - jcn::Vector{MUMPS_INT} - a::Vector{TC} - irn_loc::Vector{MUMPS_INT} - jcn_loc::Vector{MUMPS_INT} - a_loc::Vector{TC} - eltptr::Vector{MUMPS_INT} - eltvar::Vector{MUMPS_INT} - a_elt::Vector{TC} - perm_in::Vector{MUMPS_INT} - sym_in::Vector{MUMPS_INT} - uns_in::Vector{MUMPS_INT} - colsca::Vector{TR} - rowsca::Vector{TR} - rhs::Vector{TC} - redrhs::Vector{TC} - rhs_sparse::Vector{TC} - sol_loc::Vector{TC} - rhs_loc::Vector{TC} - irhs_sparse::Vector{MUMPS_INT} - irhs_ptr::Vector{MUMPS_INT} - isol_loc::Vector{MUMPS_INT} - irhs_loc::Vector{MUMPS_INT} - pivnul_list::Vector{MUMPS_INT} - mapping::Vector{MUMPS_INT} - listvar_schur::Vector{MUMPS_INT} - schur::Vector{TC} - wk_user::Vector{TC} - - GC_haven{T}() where T = new{T,real(T)}() -end - - -# structure used througout package, which contains the mumps struct as well as -# a struct to hold data to keep it safe from garbage collection (due to necessity) -# of passing pointers to MUMPS that are not safe from Julia gc -mutable struct Mumps{TC,TR} - mumpsc::MumpsC{TC,TR} - gc_haven::GC_haven{TC,TR} + gc_haven::Array{Ref,1} finalized::Bool function Mumps{T}(sym::Int,par::Int,comm) where T - mumpsc = MumpsC{T}(sym,par,comm) - invoke_mumps_unsafe!(mumpsc) - new{T,real(T)}(mumpsc,GC_haven{T}(),false) + mumps = new{T,real(T)}(sym,par,-1,comm) + invoke_mumps_unsafe!(mumps) + mumps.gc_haven = Array{Ref,1}(undef,0) + mumps.finalized = false + return mumps end end @@ -166,6 +124,6 @@ function Base.real(T::TypeVar) if T<:Number return real(T) else - throw(ArgumentError("real not defined for type $T")) + throw(DomainError("real not defined for type $T")) end end diff --git a/src/printing.jl b/src/printing.jl index fdd03a4..ecd7b1f 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -1,8 +1,6 @@ export display_icntl -Base.show(io::IO,mumps::Mumps) = show(io,mumps.mumpsc) - -function Base.show(io::IO,mumpsc::MumpsC{TC,TR}) where {TC,TR} +function Base.show(io::IO,mumps::Mumps{TC,TR}) where {TC,TR} print(io,"Mumps{$TC,$TR}: ") if TC<:Float32 println(io,"single precision real") @@ -18,61 +16,61 @@ function Base.show(io::IO,mumpsc::MumpsC{TC,TR}) where {TC,TR} lib = "zmumps" end println(io,"lib: ", lib) - print(io,"job: ", mumpsc.job, " ") - if mumpsc.job==-3 + print(io,"job: ", mumps.job, " ") + if mumps.job==-3 println(io,"save/restore") - elseif mumpsc.job==-2 + elseif mumps.job==-2 println(io,"terminate") - elseif mumpsc.job==-1 + elseif mumps.job==-1 println(io,"initialize") - elseif mumpsc.job==1 + elseif mumps.job==1 println(io,"analyze") - elseif mumpsc.job==2 + elseif mumps.job==2 println(io,"factorize") - elseif mumpsc.job==3 + elseif mumps.job==3 println(io,"solve") - elseif mumpsc.job==4 + elseif mumps.job==4 println(io,"analyze + factorize") - elseif mumpsc.job==5 + elseif mumps.job==5 println(io,"factorize + solve") - elseif mumpsc.job==6 + elseif mumps.job==6 println(io,"analyze + factorize + solve") - elseif mumpsc.job==7 + elseif mumps.job==7 println(io,"save") - elseif mumpsc.job==8 + elseif mumps.job==8 println(io,"restore") else println(io,"unrecognized") end - print(io,"sym: ", mumpsc.sym) - if mumpsc.sym==1 + print(io,"sym: ", mumps.sym) + if mumps.sym==1 println(io," symmetric pos def") - elseif mumpsc.sym==2 + elseif mumps.sym==2 println(io," symmetric") else println(io," unsymmetric") end - print(io,"par: ", mumpsc.par) - mumpsc.par==0 ? println(io," host not among workers") : println(io," host among workers ") + print(io,"par: ", mumps.par) + mumps.par==0 ? println(io," host not among workers") : println(io," host among workers ") print(io,"matrix A: ") - if has_matrix(mumpsc) - print(io,"$(mumpsc.n)×$(mumpsc.n) ") - if is_matrix_assembled(mumpsc) - println(io,"sparse matrix, with $(mumpsc.nnz) nonzero elements") + if has_matrix(mumps) + print(io,"$(mumps.n)×$(mumps.n) ") + if is_matrix_assembled(mumps) + println(io,"sparse matrix, with $(mumps.nnz) nonzero elements") else - print(io,"elemental matrix with $(mumpsc.nelt) element") - mumpsc.nelt>1 ? println(io,"s") : println() + print(io,"elemental matrix with $(mumps.nelt) element") + mumps.nelt>1 ? println(io,"s") : println() end else println(io,"uninitialized") end print(io,"rhs B:") - rhs_type = is_rhs_dense(mumpsc) ? "dense" : "sparse" - nz_rhs = is_rhs_dense(mumpsc) ? "" : string(",with ",mumpsc.nz_rhs," nonzero elements") - if has_rhs(mumpsc) - lrhs = is_rhs_dense(mumpsc) ? mumpsc.lrhs : mumpsc.n - nrhs = mumpsc.nrhs + rhs_type = is_rhs_dense(mumps) ? "dense" : "sparse" + nz_rhs = is_rhs_dense(mumps) ? "" : string(",with ",mumps.nz_rhs," nonzero elements") + if has_rhs(mumps) + lrhs = is_rhs_dense(mumps) ? mumps.lrhs : mumps.n + nrhs = mumps.nrhs println(io," $lrhs×$nrhs ",rhs_type," matrix", nz_rhs) else println(io," uninitialized") @@ -82,33 +80,7 @@ function Base.show(io::IO,mumpsc::MumpsC{TC,TR}) where {TC,TR} icntl_inds = [4,9,13,19,30,33] for i ∈ eachindex(icntl_inds) print(io,"\t") - display_icntl(io,mumpsc.icntl,icntl_inds[i],mumpsc.icntl[icntl_inds[i]]) - end -end - -function Base.show(io::IO,gc_haven::GC_haven) - vars = (:irn, :jcn, :a, - :irn_loc, :jcn_loc, :a_loc, - :eltptr, :eltvar, :a_elt, - :perm_in, :sym_in, :uns_in, - :colsca, :rowsca, - :rhs, :redrhs, :rhs_sparse, - :sol_loc, :rhs_loc, :irhs_sparse, - :irhs_ptr, :isol_loc, :irhs_loc, - :pivnul_list, :mapping, - :listvar_schur, :schur, :wk_user) - pad = (":\t",":\t",":\t\t", - ": ","\t",":\t", - ":\t","\t",":\t", - ": ","\t",":\t", - ":\t",":\t", - ":\t",":\t",":\t", - ": ",":",": ", - ": ",": ",": ", - ": ",": ", - ": ",":\t",": ") - for i ∈ eachindex(vars) - println(io,vars[i],pad[i],isdefined(gc_haven,vars[i])) + display_icntl(io,mumps.icntl,icntl_inds[i],mumps.icntl[icntl_inds[i]]) end end @@ -119,7 +91,7 @@ Show the complete ICNTL integer array of `mumps`, with descriptions See also: [`set_icntl!`](@ref) """ -display_icntl(io::IO,mumps::Mumps) = display_icntl(io,mumps.mumpsc.icntl) +display_icntl(io::IO,mumps::Mumps) = display_icntl(io,mumps.mumps.icntl) function display_icntl(io::IO,icntl) for i ∈ eachindex(icntl) display_icntl(io,icntl,i,icntl[i]) @@ -451,7 +423,7 @@ Show the complete CNTL real array of `mumps`, with descriptions See also: [`set_cntl!`](@ref) """ -display_cntl(io::IO,mumps::Mumps) = display_cntl(io,mumps.mumpsc.cntl) +display_cntl(io::IO,mumps::Mumps) = display_cntl(io,mumps.mumps.cntl) function display_cntl(io::IO,cntl) for i ∈ eachindex(cntl) display_icntl(io,cntl,i,cntl[i]) diff --git a/test/advanced_solve.jl b/test/advanced_solve.jl new file mode 100644 index 0000000..57e2f88 --- /dev/null +++ b/test/advanced_solve.jl @@ -0,0 +1,19 @@ +using MUMPS3 +using MPI, LinearAlgebra, SparseArrays + +MPI.Initialized() ? nothing : MPI.Init() + +N, M = 2000, 20 + +@testset "Advanced solve: " begin + for i ∈ 1:2000 + A = sparse(1.0im*I,N,N) + sprand(N,N,1/N) + y = sprand(N,M,1/sqrt(N*M)) + + m = Mumps(A,y) + toggle_display!(m) + + x = mumps_solve(m) + @test norm(A*x-y) ≤ sqrt(eps()) + end +end diff --git a/test/basic_solve.jl b/test/basic_solve.jl new file mode 100644 index 0000000..ae92bc4 --- /dev/null +++ b/test/basic_solve.jl @@ -0,0 +1,17 @@ +using MUMPS3 +using MPI, LinearAlgebra, SparseArrays + +MPI.Initialized() ? nothing : MPI.Init() + +N, M = 2000, 20 + +@testset "Basic solve: " begin +for i ∈ 1:2000 + A = sparse(I,N,N) + sprand(N,N,1/N) + y = sprand(N,M,1/sqrt(N*M)) + + x = mumps_solve(A,y) + + @test norm(A*x-y) ≤ sqrt(eps()) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index ed59072..702983e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1 +1,5 @@ -include("simple_test.jl") +using Test + +include("basic_solve.jl") +include("advanced_solve.jl") +include("schur_complement.jl") diff --git a/test/schur_complement.jl b/test/schur_complement.jl new file mode 100644 index 0000000..5ec8e77 --- /dev/null +++ b/test/schur_complement.jl @@ -0,0 +1,21 @@ +using MUMPS3 +using MPI, LinearAlgebra, SparseArrays + +MPI.Initialized() ? nothing : MPI.Init() + +N, M = 2000, 20 +A = sparse(I,N,N) + sprand(N,N,1/N) +invA = inv(Matrix(A)) + +@testset "Schur complement: " begin + for i ∈ 1:2000 + + y = sprand(N,M,1/sqrt(N*M)) + + S = mumps_schur_complement(A,y) + + schur_inds = unique!(sort(y.rowval)) + + @test norm(inv(S) - invA[schur_inds,schur_inds]) ≤ sqrt(eps()) + end +end diff --git a/test/simple_test.jl b/test/simple_test.jl deleted file mode 100644 index f06b784..0000000 --- a/test/simple_test.jl +++ /dev/null @@ -1,10 +0,0 @@ -using MUMPS3,MPI,LinearAlgebra,SparseArrays - -MPI.Initialized() ? nothing : MPI.Init() - -N, M = 1000, 10 -A = sparse(I,N,N) + sprand(N,N,1/N) -y = sprand(N,M,1/sqrt(N*M)) - -x = mumps_solve(A,y) -norm(A*x-y)