Skip to content

Commit

Permalink
Use the new generid version of MA27
Browse files Browse the repository at this point in the history
  • Loading branch information
@amontoison
amontoison committed Dec 24, 2024
1 parent f3d4824 commit 93dcf5f
Showing 1 changed file with 45 additions and 75 deletions.
120 changes: 45 additions & 75 deletions lib/MadNLPHSL/src/ma27.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -123,22 +93,22 @@ 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)<M.maxfrt[1] && resize!(M.w,M.maxfrt[1])
length(M.iw1)<M.nsteps[1] && resize!(M.iw1,M.nsteps[1])
ma27c!(Int32(M.csc.n),M.a,M.la,M.iw,M.liw,M.w,M.maxfrt,rhs,
M.iw1,M.nsteps,M.icntl,M.info)
HSL.ma27cr(T, INT, M.csc.n,M.a,M.la,M.iw,M.liw,M.w,M.maxfrt,rhs,
M.iw1,M.nsteps,M.icntl,M.info)
M.info[1]<0 && throw(SolveException())
return rhs
end

is_inertia(::Ma27Solver) = true
function inertia(M::Ma27Solver)
dim = M.csc.n
rank = (Int(M.info[1])==3) ? Int(M.info[2]) : dim
rank = (Int(M.info[1]) == 3) ? Int(M.info[2]) : dim
neg = Int(M.info[15])

return (rank-neg,dim-rank,neg)
return (rank-neg, dim-rank,neg)
end

function improve!(M::Ma27Solver)
Expand Down

0 comments on commit 93dcf5f

Please sign in to comment.