Use hasmethod for checking mul! implementation

src/DiagonalHessianApproximation.jl

@@ -20,26 +20,23 @@ positive definite.
 """
 mutable struct DiagonalPSB{T <: Real, I <: Integer, V <: AbstractVector{T}, F} <:
                AbstractDiagonalQuasiNewtonOperator{T}
-  d::V # Diagonal of the operator
-  nrow::I
-  ncol::I
-  symmetric::Bool
-  hermitian::Bool
-  prod!::F
-  tprod!::F
-  ctprod!::F
+  const d::V # Diagonal of the operator
+  const nrow::I
+  const ncol::I
+  const symmetric::Bool
+  const hermitian::Bool
+  const prod!::F
+  const tprod!::F
+  const ctprod!::F
   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!
-  allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
 end
 
 @doc (@doc DiagonalPSB) function DiagonalPSB(d::AbstractVector{T}) where {T <: Real}
   prod = (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
   n = length(d)
-  DiagonalPSB(d, n, n, true, true, prod, prod, prod, 0, 0, 0, true, true, true)
+  DiagonalPSB(d, n, n, true, true, prod, prod, prod, 0, 0, 0)
 end
 
 # update function
@@ -98,26 +95,23 @@ positive definite.
 """
 mutable struct DiagonalAndrei{T <: Real, I <: Integer, V <: AbstractVector{T}, F} <:
                AbstractDiagonalQuasiNewtonOperator{T}
-  d::V # Diagonal of the operator
-  nrow::I
-  ncol::I
-  symmetric::Bool
-  hermitian::Bool
-  prod!::F
-  tprod!::F
-  ctprod!::F
+  const d::V # Diagonal of the operator
+  const nrow::I
+  const ncol::I
+  const symmetric::Bool
+  const hermitian::Bool
+  const prod!::F
+  const tprod!::F
+  const ctprod!::F
   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!
-  allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
 end
 
 @doc (@doc DiagonalAndrei) function DiagonalAndrei(d::AbstractVector{T}) where {T <: Real}
   prod = (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
   n = length(d)
-  DiagonalAndrei(d, n, n, true, true, prod, prod, prod, 0, 0, 0, true, true, true)
+  DiagonalAndrei(d, n, n, true, true, prod, prod, prod, 0, 0, 0)
 end
 
 # update function
@@ -155,20 +149,17 @@ https://doi.org/10.18637/jss.v060.i03
 """
 mutable struct SpectralGradient{T <: Real, I <: Integer, F} <:
                AbstractDiagonalQuasiNewtonOperator{T}
-  d::Vector{T} # Diagonal coefficient of the operator (multiple of the identity)
-  nrow::I
-  ncol::I
-  symmetric::Bool
-  hermitian::Bool
-  prod!::F
-  tprod!::F
-  ctprod!::F
+  const d::Vector{T} # Diagonal coefficient of the operator (multiple of the identity)
+  const nrow::I
+  const ncol::I
+  const symmetric::Bool
+  const hermitian::Bool
+  const prod!::F
+  const tprod!::F
+  const ctprod!::F
   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!
-  allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
 end
 
 """
@@ -186,7 +177,7 @@ function SpectralGradient(σ::T, n::I) where {T <: Real, I <: Integer}
   @assert σ > 0
   d = [σ]
   prod = (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
-  SpectralGradient(d, n, n, true, true, prod, prod, prod, 0, 0, 0, true, true, true)
+  SpectralGradient(d, n, n, true, true, prod, prod, prod, 0, 0, 0)
 end
 
 # update function
@@ -207,37 +198,34 @@ end
 """
     DiagonalBFGS(d)
 
-A diagonal approximation of the BFGS update inspired by 
-Marnissi, Y., Chouzenoux, E., Benazza-Benyahia, A., & Pesquet, J. C. (2020). 
+A diagonal approximation of the BFGS update inspired by
+Marnissi, Y., Chouzenoux, E., Benazza-Benyahia, A., & Pesquet, J. C. (2020).
 Majorize–minimize adapted Metropolis–Hastings algorithm.
-https://ieeexplore.ieee.org/abstract/document/9050537. 
+https://ieeexplore.ieee.org/abstract/document/9050537.
 
 # Arguments
 
 - `d::AbstractVector`: initial diagonal approximation.
 """
 mutable struct DiagonalBFGS{T <: Real, I <: Integer, V <: AbstractVector{T}, F} <:
                AbstractDiagonalQuasiNewtonOperator{T}
-  d::V # Diagonal of the operator
-  nrow::I
-  ncol::I
-  symmetric::Bool
-  hermitian::Bool
-  prod!::F
-  tprod!::F
-  ctprod!::F
+  const d::V # Diagonal of the operator
+  const nrow::I
+  const ncol::I
+  const symmetric::Bool
+  const hermitian::Bool
+  const prod!::F
+  const tprod!::F
+  const ctprod!::F
   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!
-  allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
 end
 
 @doc (@doc DiagonalBFGS) function DiagonalBFGS(d::AbstractVector{T}) where {T <: Real}
   prod = (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
   n = length(d)
-  DiagonalBFGS(d, n, n, true, true, prod, prod, prod, 0, 0, 0, true, true, true)
+  DiagonalBFGS(d, n, n, true, true, prod, prod, prod, 0, 0, 0)
 end
 
 # update function
@@ -258,3 +246,10 @@ function push!(
   B.d .*= sum(B.d) / sT_y
   return B
 end
+
+for op in (DiagonalPSB, DiagonalAndrei, SpectralGradient, DiagonalBFGS)
+  @eval begin
+    isallocated5(::$op) = true
+    has_args5(::$op) = true
+  end
+end

src/TimedOperators.jl

@@ -5,10 +5,10 @@ export TimedLinearOperator
 mutable struct TimedLinearOperator{T, OP <: AbstractLinearOperator{T}, F, Ft, Fct} <:
                AbstractLinearOperator{T}
   timer::TimerOutput
-  op::OP
-  prod!::F
-  tprod!::Ft
-  ctprod!::Fct
+  const op::OP
+  const prod!::F
+  const tprod!::Ft
+  const ctprod!::Fct
 end
 
 TimedLinearOperator{T}(
@@ -26,9 +26,9 @@ Creates a linear operator instrumented with timers from TimerOutputs.
 """
 function TimedLinearOperator(op::AbstractLinearOperator{T}) where {T}
   timer = TimerOutput()
-  prod!(res, x, α, β) = @timeit timer "prod" op.prod!(res, x, α, β)
-  tprod!(res, x, α, β) = @timeit timer "tprod" op.tprod!(res, x, α, β)
-  ctprod!(res, x, α, β) = @timeit timer "ctprod" op.ctprod!(res, x, α, β)
+  prod!(res, x, α, β) = @timeit timer "prod" mul!(res, op, x, α, β)
+  tprod!(res, x, α, β) = @timeit timer "tprod" mul!(res, transpose(op),x, α, β)
+  ctprod!(res, x, α, β) = @timeit timer "ctprod" mul!(res, op', x, α, β)
   TimedLinearOperator{T}(timer, op, prod!, tprod!, ctprod!)
 end
 
@@ -42,7 +42,6 @@ for fn ∈ (
   :issymmetric,
   :ishermitian,
   :has_args5,
-  :use_prod5!,
   :isallocated5,
   :allocate_vectors_args3!,
   :nprod,

src/abstract.jl

@@ -7,14 +7,14 @@ export AbstractLinearOperator,
   ishermitian,
   symmetric,
   issymmetric,
-  has_args5,
+  has_args5,     # TODO: deprecate
   isallocated5,
   nprod,
   ntprod,
   nctprod,
   reset!
 
-mutable struct LinearOperatorException <: Exception
+struct LinearOperatorException <: Exception
   msg::AbstractString
 end
 
@@ -44,21 +44,18 @@ other operators, with matrices and with scalars. Operators may
 be transposed and conjugate-transposed using the usual Julia syntax.
 """
 mutable struct LinearOperator{T, S, I <: Integer, F, Ft, Fct} <: AbstractLinearOperator{T}
-  nrow::I
-  ncol::I
-  symmetric::Bool
-  hermitian::Bool
-  prod!::F
-  tprod!::Ft
-  ctprod!::Fct
+  const nrow::I
+  const ncol::I
+  const symmetric::Bool
+  const hermitian::Bool
+  const prod!::F
+  const tprod!::Ft
+  const 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
+  Mv::S
+  Mtu::S
 end
 
 function LinearOperator{T, S}(
@@ -73,12 +70,6 @@ function LinearOperator{T, S}(
   ntprod::I,
   nctprod::I,
 ) where {T, S, I <: Integer, F, Ft, Fct}
-  Mv5, Mtu5 = S(undef, 0), S(undef, 0)
-  nargs = get_nargs(prod!)
-  args5 = (nargs == 4)
-  (args5 == false) || (nargs != 2) || throw(LinearOperatorException("Invalid number of arguments"))
-  allocated5 = args5 ? true : false
-  use_prod5! = args5 ? true : false
   return LinearOperator{T, S, I, F, Ft, Fct}(
     nrow,
     ncol,
@@ -90,11 +81,8 @@ function LinearOperator{T, S}(
     nprod,
     ntprod,
     nctprod,
-    args5,
-    use_prod5!,
-    Mv5,
-    Mtu5,
-    allocated5,
+    S(undef, 0),
+    S(undef, 0),
   )
 end
 
@@ -147,56 +135,7 @@ LinearOperator{T}(
   S::Type = Vector{T},
 ) where {T} = LinearOperator{T, S}(nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!)
 
-# create operator from other operators with +, *, vcat,...
-# TODO: this is not a type, so it should not be uppercase
-function CompositeLinearOperator(
-  T::Type,
-  nrow::I,
-  ncol::I,
-  symmetric::Bool,
-  hermitian::Bool,
-  prod!::F,
-  tprod!::Ft,
-  ctprod!::Fct,
-  args5::Bool,
-  ::Type{S},
-) where {S, I <: Integer, F, Ft, Fct}
-  Mv5, Mtu5 = S(undef, 0), S(undef, 0)
-  allocated5 = true
-  use_prod5! = true
-  return LinearOperator{T, S, I, F, Ft, Fct}(
-    nrow,
-    ncol,
-    symmetric,
-    hermitian,
-    prod!,
-    tprod!,
-    ctprod!,
-    0,
-    0,
-    0,
-    args5,
-    use_prod5!,
-    Mv5,
-    Mtu5,
-    allocated5,
-  )
-end
-
-# backward compatibility (not inferrable)
-CompositeLinearOperator(
-  T::Type,
-  nrow::I,
-  ncol::I,
-  symmetric::Bool,
-  hermitian::Bool,
-  prod!::F,
-  tprod!::Ft,
-  ctprod!::Fct,
-  args5::Bool;
-  S::Type = Vector{T},
-) where {I <: Integer, F, Ft, Fct} =
-  CompositeLinearOperator(T, nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, args5, S)
+const CompositeLinearOperator = LinearOperator   # backwards compatibility
 
 nprod(op::AbstractLinearOperator) = op.nprod
 ntprod(op::AbstractLinearOperator) = op.ntprod
@@ -213,18 +152,24 @@ Determine whether the operator can work with the 5-args `mul!`.
 If `false`, storage vectors will be generated at the first call of
 the 5-args `mul!`.
 No additional vectors are generated when using the 3-args `mul!`.
+
+!!! warning
+    `has_args5` can be very slow. A better option is to use Julia's `hasmethod`
+    at points in the code where the concrete types of objects used in `mul!` are known.
+
+    `has_args5` may be removed in a future release.
 """
-has_args5(op::AbstractLinearOperator) = op.args5
-use_prod5!(op::AbstractLinearOperator) = op.use_prod5!
-isallocated5(op::AbstractLinearOperator) = op.allocated5
+has_args5(op::AbstractLinearOperator) = get_nargs(op.prod!) == 4
+
+isallocated5(op::LinearOperator) = !(isempty(op.Mv) || isempty(op.Mtu))
 
 has_args5(op::AbstractMatrix) = true  # Needed for BlockDiagonalOperator
 
 # Alert user of the need for storage_type method definition for arbitrary, user defined operators
 storage_type(op::AbstractLinearOperator) = error("please implement storage_type for $(typeof(op))")
 
-storage_type(op::LinearOperator) = typeof(op.Mv5)
-storage_type(M::AbstractMatrix{T}) where {T} = Vector{T}
+storage_type(::LinearOperator{T, S}) where {T, S} = S
+storage_type(::AbstractMatrix{T}) where {T} = Vector{T}
 
 # Lazy wrappers
 storage_type(op::Adjoint) = storage_type(parent(op))