3-args operators

[geoffroyleconte]

@dpo @amontoison what do you think about this to define operators that only use the 3-args mul!? I can change the name of the constructor Should I also add the @closure macro to the functions inside the new constructor?

[codecov[bot]]

Codecov Report

Merging #190 (191b6ba) into main (85c85d3) will increase coverage by 0.29%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main     #190      +/-   ##
==========================================
+ Coverage   97.45%   97.75%   +0.29%     
==========================================
  Files          13       13              
  Lines         787      890     +103     
==========================================
+ Hits          767      870     +103     
  Misses         20       20              

Impacted Files	Coverage Δ
src/TimedOperators.jl	93.33% <ø> (+0.47%)	:arrow_up:
src/abstract.jl	100.00% <100.00%> (ø)
src/adjtrans.jl	99.07% <100.00%> (+0.14%)	:arrow_up:
src/cat.jl	100.00% <100.00%> (ø)
src/constructors.jl	93.33% <100.00%> (ø)
src/lbfgs.jl	97.14% <100.00%> (+0.30%)	:arrow_up:
src/lsr1.jl	97.95% <100.00%> (+0.28%)	:arrow_up:
src/operations.jl	100.00% <100.00%> (ø)
src/special-operators.jl	99.18% <100.00%> (+0.07%)	:arrow_up:
... and 8 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 85c85d3...191b6ba. Read the comment docs.

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCI.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[amontoison]

@geoffroyleconte, I discovered the applicable function today. It should allow us to use the same constructor LinearOperator(....) for both operators. Instead of using a args5 boolean we could check that 3-args / 5-args mul! are working or not with applicable.

[geoffroyleconte]

@geoffroyleconte, I discovered the applicable function today. It should allow us to use the same constructor LinearOperator(....) for both operators. Instead of using a args5 boolean we could check that 3-args / 5-args mul! are working or not with applicable.

It is not very fast: for an operation of size 10x10:

julia> A = rand(10,10); b = rand(10); res = rand(10);
julia> @benchmark applicable(mul!, res, A, b)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     80.209 ns (0.00% GC)
  median time:      89.110 ns (0.00% GC)
  mean time:        91.245 ns (0.00% GC)
  maximum time:     194.764 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     955

julia> @benchmark mul!(res, A, b)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     114.521 ns (0.00% GC)
  median time:      122.992 ns (0.00% GC)
  mean time:        125.009 ns (0.00% GC)
  maximum time:     320.792 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     909

[amontoison]

I don't understand @geoffroyleconte, with applicable it's faster.

[geoffroyleconte]

I don't understand @geoffroyleconte, with applicable it's faster.

applicable only returns a boolean saying that it's possible to call mul!, but you still have to call mul! after.

[amontoison]

I don't understand @geoffroyleconte, with applicable it's faster.

applicable only returns a boolean saying that it's possible to call mul!, but you still have to call mul! after.

Oh ok. But can't we use applicable inside a method has_arg5(op) and only use it to define the args5 attribute?

[geoffroyleconte]

I don't understand @geoffroyleconte, with applicable it's faster.

applicable only returns a boolean saying that it's possible to call mul!, but you still have to call mul! after.

Oh ok. But can't we use applicable inside a method has_arg5(op) and only use it to define the args5 attribute?

To do so I think you'll need vectors, but we do not want to pre-allocate vectors. And the use of applicable would still be slow for small problems. In my opinion, the use of a separate constructor is a better option, and it is more explicit

[amontoison]

I don't understand @geoffroyleconte, with applicable it's faster.

applicable only returns a boolean saying that it's possible to call mul!, but you still have to call mul! after.

We can call it at the first mul! with the function proposed by Dominique (has_arg5(op)) to determine what kind of operator we have and save it as an attribute of the operator.

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCI.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[amontoison]

@geoffroyleconte bump

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[amontoison]

Is it possible to have LinearOperator instead of LinearOperator_args3 for the name of the constructor ?

[geoffroyleconte]

@amontoison @dpo I rebased, if you are OK with the constructor I'll add some documentation about it

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[geoffroyleconte]

@amontoison I don't think your changes would work with functions such as the ones in cat.jl and many others. It would require duplicating a lot of code for the 3-args mul!. I think that the easiest way is still to "convert" the 3-args mul! functions to 5-args.

[amontoison]

@amontoison I don't think your changes would work with functions such as the ones in cat.jl and many others. It would require duplicating a lot of code for the 3-args mul!. I think that the easiest way is still to "convert" the 3-args mul! functions to 5-args.

Ok, I understand the problem. I would like to avoid the compilation of the two functions tprod! and tprod5!. We can add a @inline macro or keep things like this.

[amontoison]

@geoffroyleconte Can you add the option args5 for BlockLinearOperator ? It is relevant for block-Jacobi preconditioners.

[geoffroyleconte]

@amontoison I don't think it's needed, you just have to specify it to the operators you're using to construct your BlockDiagonalOperator

[amontoison]

You will an issue here https://github.com/JuliaSmoothOptimizers/LinearOperators.jl/blob/main/src/special-operators.jl#L275, no ?

[amontoison]

I give you an example :

function block_jacobi(A, b, N)
  ops = LinearOperator[]
  graph = Metis.graph(A, check_hermitian=false)
  partitions = Metis.partition(graph, N)
  size_partition = [count(partitions .== i) for i=1:N]

  p = sortperm(partitions)
  Abis = A[p,p]
  bbis = b[p]

  debut = 1
  for i = 1 : N
    np = size_partition[i]
    fin = debut + np - 1
    D = Abis[debut:fin, debut:fin]
    F = lu(D)
    push!(ops, LinearOperator(Float64, np, np, false, false, (y, v) -> ldiv!(y, F, v)), args5=false))
    debut = fin + 1
  end
  return Abis, bbis, BlockDiagonalOperator(ops...)
end

A = sparse(rand(10,10))
b = rand(10)
N = 3
Abis, bbis, P = block_jacobi(Abis, bbis, N)
Matrix(P * Abis)

[geoffroyleconte]

using LinearOperators, LinearAlgebra
n = 10;
A, A2 = rand(n, n), rand(n, n);
opA = LinearOperator(Float64, n, n, false, false, (res, v) -> mul!(res, A, v), args5 = false);
opA2 = LinearOperator(Float64, n, n, false, false, (res, v) -> mul!(res, A2, v), args5 = false);
opAbloc = BlockDiagonalOperator(opA, opA2)
vb, vb2 = rand(2*n), rand(2*n)

julia> mul!(vb2, opAbloc, vb, 1.0, 0.0)

julia> mul!(vb2, opAbloc, vb, 1.0, 1.0)
ERROR: LinearOperatorException("5-args product not defined")
Stacktrace:
 [1] mul!
   @ c:\Users\Geoffroy Leconte\.julia\dev\LinearOperators\src\operations.jl:4 [inlined]
 [2] (::LinearOperators.var"#prod!#135"{Tuple{LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#5#6"}, Nothing, Nothing}, LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#7#8"}, Nothing, Nothing}}})(y::Vector{Float64}, x::Vector{Float64}, α::Float64, β::Float64)
   @ LinearOperators c:\Users\Geoffroy Leconte\.julia\dev\LinearOperators\src\special-operators.jl:245
 [3] mul!(res::Vector{Float64}, op::LinearOperator{Float64, Int64, LinearOperators.var"#prod!#135"{Tuple{LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#5#6"}, Nothing, Nothing}, LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#7#8"}, Nothing, Nothing}}}, LinearOperators.var"#tprod!#136"{Tuple{LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#5#6"}, Nothing, Nothing}, LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#7#8"}, Nothing, Nothing}}}, LinearOperators.var"#ctprod!#137"{Tuple{LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#5#6"}, Nothing, Nothing}, LinearOperator{Float64, Int64, LinearOperators.var"#10#13"{var"#7#8"}, Nothing, Nothing}}}}, v::Vector{Float64}, α::Float64, β::Float64)
   @ LinearOperators c:\Users\Geoffroy Leconte\.julia\dev\LinearOperators\src\operations.jl:7
 [4] top-level scope
   @ REPL[31]:1

Is this the behaviour you want for BlockDiagonalOperators? The first mul! works because it is actually a 3-args mul! (alpha=1, beta=0)

[amontoison]

Yes, If I only use 3args linear operators, I want a 3args block diagonal operator and check it with has_5args() function. has_5args(opAbloc) should return false and not true here.

[geoffroyleconte]

I can add this but there are a fiew LinearOperators for which this will not work by construction (sum of 2 operators, - operator). I will specify this in the documentation for the 3-args mul!

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[amontoison]

LGTM, if @dpo is ok we can merge it.

[dpo]

I've no major complaints here but I really find it clunky that we have to check whether an operator defines mul5 or not. It seems to me that it should be possible for operators to always have mul5, whether the user instantiates the operator with a 3-arg mul or not.

[geoffroyleconte]

It's not always possible, using functions to define an operator such as ldiv! for example you would have to allocate another vector

[dpo]

Sorry I don't understand.

[geoffroyleconte]

In some cases you cannot define directly a 5-args LinearOperator without allocating an additional vector

[dpo]

That's ok; you can define a closure.

[amontoison]

I've no major complaints here but I really find it clunky that we have to check whether an operator defines mul5 or not. It seems to me that it should be possible for operators to always have mul5, whether the user instantiates the operator with a 3-arg mul or not.

I agree with you Dominique, but this PR gives the possibility to define a 3-arg mul!, which a missing feature. I propose to define mul5 for all operators with an other PR. It will be easier to review too. :smiley:

[geoffroyleconte]

@dpo to do what you want I will need to add 3 more vectors in the non symmetric case to the LinearOperator struct. Since we have

julia> isconcretetype(Union{Nothing, Vector{Float64}})
false

, and we need concrete types to maximize the performance (I think), I will have to define something like:

mutable struct LinearOperator{T, I <: Integer, F, Ft, Fct, S} <: AbstractLinearOperator{T}
  nrow::I
  ncol::I
  symmetric::Bool
  hermitian::Bool
  prod!::F
  tprod!::Ft
  ctprod!::Fct
  nprod::I
  ntprod::I
  nctprod::I
  args5::Bool
  Mv5::S
  Mtu5::S
  Maw5::S
end

and in the case of a 5-args mul! already defined by the user, initialize

Mv5, Mtu5, Maw5 = T[], T[], T[]

Do you agree with this solution?

Another issue is that we loose 3 vectors in the case of a 3-args mul! defined by the user, but maybe I can initialize some empty vectors as well, and change them at the first use of a 5-args mul! function

[dpo]

Yes. Let's not be worried about 3 vectors.

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[geoffroyleconte]

@dpo @amontoison I made the changes, and I added a bool allocated5 to see if some vectors have already been allocated when using the 3-args mul!. The bool use_prod5! is usefull when creating LinearOperators from operations. I created a new constructor CompositeLinearOperator so that it is clearer, but I did not export it. For example, if opA = opA1 + opA2 and opA1 and/or opA2 have args5 = false, then opA also has args5 = false, because at some point if the user wants to use 5-args mul! he will need to allocated vectors, but use_prod5 = true because the operation defined in + is a 5-args mul!.

From the user point of view, he should only be aware that if he uses an operator op, and that has_args5(op) == false, then there will be allocations if he uses a 5-args mul!

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[github-actions[bot]]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
JSOSolvers.jl
Krylov.jl
NLPModels.jl
NLPModelsModifiers.jl
PROPACK.jl
Percival.jl
QuadraticModels.jl
SolverTools.jl

[geoffroyleconte]

@dpo I added get_nargs following your comments

[dpo]

OUF !!!