diff --git a/lib/MadNLPHSL/src/ma27.jl b/lib/MadNLPHSL/src/ma27.jl index ddac4631..455a66ba 100644 --- a/lib/MadNLPHSL/src/ma27.jl +++ b/lib/MadNLPHSL/src/ma27.jl @@ -1,117 +1,87 @@ -ma27_default_icntl() = Int32[ +ma27_default_icntl(INT) = INT[ 6,6,0,2139062143,1,32639,32639,32639,32639,14,9,8,8,9,10,32639,32639,32639,32689,24,11,9,8,9,10,0,0,0,0,0] -ma27_default_cntl(T) = T[.1,1.0,0.,0.,0.] +ma27_default_cntl(T) = T[.1,1.0,0.,0.,0.] @kwdef mutable struct Ma27Options <: AbstractOptions ma27_pivtol::Float64 = 1e-8 ma27_pivtolmax::Float64 = 1e-4 ma27_liw_init_factor::Float64 = 5. - ma27_la_init_factor::Float64 =5. - ma27_meminc_factor::Float64 =2. + ma27_la_init_factor::Float64 = 5. + ma27_meminc_factor::Float64 = 2. end -mutable struct Ma27Solver{T} <: AbstractLinearSolver{T} - csc::SparseMatrixCSC{T,Int32} - I::Vector{Int32} - J::Vector{Int32} +mutable struct Ma27Solver{T,INT} <: AbstractLinearSolver{T} + csc::SparseMatrixCSC{T,INT} + I::Vector{INT} + J::Vector{INT} - icntl::Vector{Int32} + icntl::Vector{INT} cntl::Vector{T} - info::Vector{Int32} + info::Vector{INT} a::Vector{T} a_view::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true} - la::Int32 - ikeep::Vector{Int32} + la::INT + ikeep::Vector{INT} - iw::Vector{Int32} - liw::Int32 - iw1::Vector{Int32} - nsteps::Vector{Int32} + iw::Vector{INT} + liw::INT + iw1::Vector{INT} + nsteps::Vector{INT} w::Vector{T} - maxfrt::Vector{Int32} + maxfrt::Vector{INT} opt::Ma27Options logger::MadNLPLogger end - -for (fa, fb, fc, typ) in [ - (:ma27aq, :ma27bq, :ma27cq, Float128), - (:ma27ad, :ma27bd, :ma27cd, Float64), - (:ma27a , :ma27b , :ma27c , Float32) - ] - @eval begin - ma27a!( - n::Cint,nz::Cint,I::Vector{Cint},J::Vector{Cint}, - iw::Vector{Cint},liw::Cint,ikeep::Vector{Cint},iw1::Vector{Cint}, - nsteps::Vector{Cint},iflag::Cint,icntl::Vector{Cint},cntl::Vector{$typ}, - info::Vector{Cint},ops::$typ - ) = HSL.$fa(n,nz,I,J,iw,liw,ikeep,iw1,nsteps,iflag,icntl,cntl,info,ops) - - ma27b!( - n::Cint,nz::Cint,I::Vector{Cint},J::Vector{Cint}, - a::Vector{$typ},la::Cint,iw::Vector{Cint},liw::Cint, - ikeep::Vector{Cint},nsteps::Vector{Cint},maxfrt::Vector{Cint},iw1::Vector{Cint}, - icntl::Vector{Cint},cntl::Vector{$typ},info::Vector{Cint} - ) = HSL.$fb(n,nz,I,J,a,la,iw,liw,ikeep,nsteps,maxfrt,iw1,icntl,cntl,info) - - ma27c!( - n::Cint,a::Vector{$typ},la::Cint,iw::Vector{Cint}, - liw::Cint,w::Vector{$typ},maxfrt::Vector{Cint},rhs::Vector{$typ}, - iw1::Vector{Cint},nsteps::Vector{Cint},icntl::Vector{Cint}, - info::Vector{Cint} - ) = HSL.$fc(n,a,la,iw,liw,w,maxfrt,rhs,iw1,nsteps,icntl,info) - end -end - -function Ma27Solver(csc::SparseMatrixCSC{T}; +function Ma27Solver(csc::SparseMatrixCSC{T,INT}; opt=Ma27Options(),logger=MadNLPLogger(), -) where T +) where {T,INT} I,J = findIJ(csc) - nz=Int32(nnz(csc)) + nz = nnz(csc) |> INT - liw= Int32(2*(2*nz+3*csc.n+1)) - iw = Vector{Int32}(undef,liw) - ikeep= Vector{Int32}(undef,3*csc.n) - iw1 = Vector{Int32}(undef,2*csc.n) - nsteps=Int32[1] - iflag =Int32(0) + liw = INT(2*(2*nz+3*csc.n+1)) + iw = Vector{INT}(undef,liw) + ikeep = Vector{INT}(undef,3*csc.n) + iw1 = Vector{INT}(undef,2*csc.n) + nsteps = INT[1] + iflag = INT(0) - icntl= ma27_default_icntl() + icntl = ma27_default_icntl(INT) cntl = ma27_default_cntl(T) icntl[1:2] .= 0 cntl[1] = opt.ma27_pivtol - info = Vector{Int32}(undef,20) - ma27a!(Int32(csc.n),nz,I,J,iw,liw,ikeep,iw1,nsteps,Int32(0),icntl,cntl,info,zero(T)) + info = Vector{INT}(undef,20) + HSL.ma27ar(T, INT, csc.n,nz,I,J,iw,liw,ikeep,iw1,nsteps,iflag,icntl,cntl,info,zero(T)) info[1]<0 && throw(SymbolicException()) - la = ceil(Int32,max(nz,opt.ma27_la_init_factor*info[5])) + la = ceil(INT,max(nz,opt.ma27_la_init_factor*info[5])) a = Vector{T}(undef,la) a_view = view(a,1:nnz(csc)) # _madnlp_unsafe_wrap is not used because we may resize a - liw= ceil(Int32,opt.ma27_liw_init_factor*info[6]) + liw= ceil(INT,opt.ma27_liw_init_factor*info[6]) resize!(iw,liw) - maxfrt=Int32[1] + maxfrt=INT[1] - return Ma27Solver{T}(csc,I,J,icntl,cntl,info,a,a_view,la,ikeep,iw,liw, + return Ma27Solver{T,INT}(csc,I,J,icntl,cntl,info,a,a_view,la,ikeep,iw,liw, iw1,nsteps,Vector{T}(),maxfrt,opt,logger) end -function factorize!(M::Ma27Solver) - M.a_view.=M.csc.nzval +function factorize!(M::Ma27Solver{T,INT}) where {T, INT} + M.a_view .= M.csc.nzval while true - ma27b!(Int32(M.csc.n),Int32(nnz(M.csc)),M.I,M.J,M.a,M.la, - M.iw,M.liw,M.ikeep,M.nsteps,M.maxfrt, - M.iw1,M.icntl,M.cntl,M.info) + HSL.ma27br(T, INT, M.csc.n,nnz(M.csc),M.I,M.J,M.a,M.la, + M.iw,M.liw,M.ikeep,M.nsteps,M.maxfrt, + M.iw1,M.icntl,M.cntl,M.info) if M.info[1] == -3 - M.liw = ceil(Int32,M.opt.ma27_meminc_factor*M.liw) + M.liw = ceil(INT,M.opt.ma27_meminc_factor*M.liw) resize!(M.iw, M.liw) @debug(M.logger,"Reallocating memory: liw ($(M.liw))") elseif M.info[1] == -4 - M.la = ceil(Int32,M.opt.ma27_meminc_factor*M.la) + M.la = ceil(INT,M.opt.ma27_meminc_factor*M.la) resize!(M.a,M.la) @debug(M.logger,"Reallocating memory: la ($(M.la))") elseif M.info[1] < 0 @@ -123,11 +93,11 @@ function factorize!(M::Ma27Solver) return M end -function solve!(M::Ma27Solver{T},rhs::Vector{T}) where T +function solve!(M::Ma27Solver{T,INT}, rhs::Vector{T}) where {T,INT} length(M.w)