From 149fe0a2f8c885c43a27abd3475dd941a2ee495e Mon Sep 17 00:00:00 2001 From: Abel Soares Siqueira Date: Fri, 23 Sep 2022 18:22:02 +0200 Subject: [PATCH 1/2] [WIP] Create specific operator for hprod --- src/NLPModels.jl | 1 + src/nlp/api.jl | 11 +---- src/nlp/operator.jl | 111 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 113 insertions(+), 10 deletions(-) create mode 100644 src/nlp/operator.jl diff --git a/src/NLPModels.jl b/src/NLPModels.jl index 17299f9b..d5e2ea31 100644 --- a/src/NLPModels.jl +++ b/src/NLPModels.jl @@ -37,6 +37,7 @@ Base type for a nonlinear least-squares model. """ abstract type AbstractNLSModel{T, S} <: AbstractNLPModel{T, S} end +include("nlp/operator.jl") for f in ["utils", "api", "counters", "meta", "show", "tools"] include("nlp/$f.jl") include("nls/$f.jl") diff --git a/src/nlp/api.jl b/src/nlp/api.jl index 0b03ccf0..b78adb5b 100644 --- a/src/nlp/api.jl +++ b/src/nlp/api.jl @@ -1225,16 +1225,7 @@ function hess_op!( obj_weight::Real = one(eltype(x)), ) @lencheck nlp.meta.nvar x Hv - prod! = @closure (res, v, α, β) -> begin - hprod!(nlp, x, v, Hv; obj_weight = obj_weight) - if β == 0 - @. res = α * Hv - else - @. res = α * Hv + β * res - end - return res - end - return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar, true, true, prod!, prod!, prod!) + return HprodOperator!(nlp, copy(x), Hv, obj_weight=obj_weight) end """ diff --git a/src/nlp/operator.jl b/src/nlp/operator.jl new file mode 100644 index 00000000..b7b8f194 --- /dev/null +++ b/src/nlp/operator.jl @@ -0,0 +1,111 @@ +using FastClosures, LinearOperators +import LinearOperators.AbstractLinearOperator + +mutable struct ModelOperator{T, I <: Integer, F, Ft, Fct, S, M <: AbstractNLPModel{T, S}} <: + AbstractLinearOperator{T} + x::S + nlp::M + nrow::I + ncol::I + symmetric::Bool + hermitian::Bool + prod!::F + tprod!::Ft + ctprod!::Fct + nprod::I + ntprod::I + nctprod::I + args5::Bool + use_prod5!::Bool # true for 5-args mul! and for composite operators created with operators that use the 3-args mul! + Mv5::S + Mtu5::S + allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated +end + +function ModelOperator( + x::S, + nlp::M, + nrow::I, + ncol::I, + symmetric::Bool, + hermitian::Bool, + prod!::F, + tprod!::Ft, + ctprod!::Fct, + nprod::I, + ntprod::I, + nctprod::I; +) where {T, I <: Integer, F, Ft, Fct, S, M <: AbstractNLPModel{T, S}} + Mv5, Mtu5 = S(undef, 0), S(undef, 0) + nargs = first(methods(prod!)).nargs - 1 + args5 = (nargs == 4) + (args5 == false) || (nargs != 2) || throw(LinearOperatorException("Invalid number of arguments")) + allocated5 = args5 ? true : false + use_prod5! = args5 ? true : false + return ModelOperator{T, I, F, Ft, Fct, S, M}( + x, + nlp, + nrow, + ncol, + symmetric, + hermitian, + prod!, + tprod!, + ctprod!, + nprod, + ntprod, + nctprod, + args5, + use_prod5!, + Mv5, + Mtu5, + allocated5, + ) +end + +ModelOperator( + x::S, + nlp::M, + nrow::I, + ncol::I, + symmetric::Bool, + hermitian::Bool, + prod!, + tprod!, + ctprod!, +) where {T, I <: Integer, S, M <: AbstractNLPModel{T, S}} = + ModelOperator(x, nlp, nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, 0, 0, 0) + +function HprodOperator!( + nlp::M, + x::S, + Hv::S; + obj_weight::Real = one(T), +) where {T, S, M <: AbstractNLPModel{T, S}} + prod! = @closure (res, v, α, β) -> begin + hprod!(nlp, x, v, Hv; obj_weight = obj_weight) + if β == 0 + @. res = α * Hv + else + @. res = α * Hv + β * res + end + return res + end + + return ModelOperator(x, nlp, nlp.meta.nvar, nlp.meta.nvar, true, true, prod!, prod!, prod!) +end + +function HprodOperator!( + nlp::M, + Hv::S; + obj_weight::Real = one(T), +) where {T, S, M <: AbstractNLPModel{T, S}} + x = copy(nlp.meta.x0) + HprodOperator!(nlp, x, Hv; obj_weight = obj_weight) +end + +function HprodOperator(nlp::M; obj_weight::Real = one(T)) where {T, S, M <: AbstractNLPModel{T, S}} + x = copy(nlp.meta.x0) + Hv = S(undef, nlp.meta.nvar) + HprodOperator!(nlp, x, Hv; obj_weight = obj_weight) +end From 7ceda7611d770c618b9da48cee300d14d0965972 Mon Sep 17 00:00:00 2001 From: Abel Soares Siqueira Date: Wed, 5 Oct 2022 00:35:00 +0200 Subject: [PATCH 2/2] Export ModelOperator and create update! for ModelOperator --- src/nlp/operator.jl | 6 ++++++ test/nlp/api.jl | 2 ++ test/nlp/dummy-model.jl | 2 +- 3 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/nlp/operator.jl b/src/nlp/operator.jl index b7b8f194..571ccee1 100644 --- a/src/nlp/operator.jl +++ b/src/nlp/operator.jl @@ -1,3 +1,5 @@ +export ModelOperator, update! + using FastClosures, LinearOperators import LinearOperators.AbstractLinearOperator @@ -76,6 +78,10 @@ ModelOperator( ) where {T, I <: Integer, S, M <: AbstractNLPModel{T, S}} = ModelOperator(x, nlp, nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, 0, 0, 0) +function update!(op::ModelOperator, x) + op.x .= x +end + function HprodOperator!( nlp::M, x::S, diff --git a/test/nlp/api.jl b/test/nlp/api.jl index 6d7d209e..c4ba855a 100644 --- a/test/nlp/api.jl +++ b/test/nlp/api.jl @@ -121,6 +121,8 @@ @test hprod!(nlp, hess_structure(nlp)..., hess_coord(nlp, x), v, Hv) ≈ H(x) * v Hop = hess_op(nlp, x) @test Hop * v ≈ H(x) * v + update!(Hop, -x) + @test Hop * v ≈ H(-x) * v Hop = hess_op!(nlp, x, Hv) @test Hop * v ≈ H(x) * v z = ones(nlp.meta.nvar) diff --git a/test/nlp/dummy-model.jl b/test/nlp/dummy-model.jl index aa0983dd..79b57fb3 100644 --- a/test/nlp/dummy-model.jl +++ b/test/nlp/dummy-model.jl @@ -26,7 +26,7 @@ end @test_throws(MethodError, jth_hess_coord!(model, [0.0], 1)) @test_throws(MethodError, jth_hprod!(model, [0.0], [1.0], 2, [3.0])) @test_throws(MethodError, ghjvprod!(model, [0.0], [1.0], [2.0], [3.0])) - @assert isa(hess_op(model, [0.0]), LinearOperator) + @assert isa(hess_op(model, [0.0]), ModelOperator) @assert isa(jac_op(model, [0.0]), LinearOperator) @assert isa(jac_lin_op(model, [0.0]), LinearOperator) @assert isa(jac_nln_op(model, [0.0]), LinearOperator)