Remove A.tprod

[amontoison]

close #135

I remove A.tprod but I don't understand something.
In the REPL with the following code :

using LinearOperators
M = rand(100,100)
x = rand(100)
A = PreallocatedLinearOperator(M)
A2 = A'
@allocated A2 * x
A3 = A'
@allocated A3 * x

Nothing is allocated when I compute A3 * x (cost of compilation of * method) whereas Aᵀ * x seems to be recompiled inside methods even if the method was compiled before.
And Aᵀ = A' is not allocating if the function has been "warm up"...

[dpo]

whereas Aᵀ * x seems to be recompiled inside methods even if the method was compiled before.
And Aᵀ = A' is not allocating if the function has been "warm up"...

I don't understand. Could you give an example?

[amontoison]

This example shows that we don't allocate with the new transpose linear operator.

using LinearOperators

n  = 10^4
A  = rand(n, n)
x  = rand(n)
y  = rand(n)
op = PreallocatedLinearOperator(A) 

function test1(op, x, y, p=10)
	opᵀ = op'
	for i = 1:p
		p = op  * x
		q = opᵀ * y
	end
end

function test2(op, x, y, p=10)
	for i = 1:p
	  p = op * x
	  q = op.tprod(y)
	end
end

test1(op, x, y)  # warm up 
@allocated test1(op, x, y)

test2(op, x, y)  # warm up
@allocated test2(op, x, y)

m = 100
A2 = rand(m, n)
op2 = PreallocatedLinearOperator(A2)
x2 = rand(n)
y2 = rand(m)

test1(op2, x2, y2)  # warm up 
@allocated test1(op2, x2, y2)

test2(op2, x2, y2)  # warm up
@allocated test2(op2, x2, y2)

[amontoison]

And this example shows that something else is allocated whereas I only changed A.tprod between crmr1 and crmr2...

using LinearOperators, Printf, LinearAlgebra, SparseArrays

include("aux.jl")
include("crmr1.jl")
include("crmr2.jl")

L   = get_div_grad(32, 32, 32)
n   = size(L, 1)
m   = div(n, 2)
A   = PreallocatedLinearOperator(L)        # Dimension n x n
Au  = PreallocatedLinearOperator(L[1:m,:]) # Dimension m x n
Ao  = PreallocatedLinearOperator(L[:,1:m]) # Dimension n x m
b   = ones(n)                              # Dimension n
c   = ones(m)                              # Dimension m

crmr1(Au, c)  # warmup
@allocated crmr1(Au, c)

crmr2(Au, c)  # warmup
@allocated crmr2(Au, c)

files.zip

[dpo]

According to the new TimedLinearOperators, there is no extra allocation:

julia> op = TimedLinearOperator(A);
julia> crmr1(op, b);
julia> op.timer
──────────────────────────────────────────────────────────────────
                           Time                   Allocations
                   ──────────────────────   ───────────────────────
 Tot / % measured:      7.02s / 3.72%            974KiB / 0.71%

 Section   ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────
 prod         441    139ms  53.1%   315μs     0.00B  0.00%    0.00B
 tprod        442    123ms  46.9%   277μs   6.91KiB  100%     16.0B
 ──────────────────────────────────────────────────────────────────
julia> op = TimedLinearOperator(A);
julia> crmr2(op, b);
julia> op.timer
──────────────────────────────────────────────────────────────────
                           Time                   Allocations
                   ──────────────────────   ───────────────────────
 Tot / % measured:      7.22s / 3.40%           16.5MiB / 0.04%

 Section   ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────
 prod         441    130ms  53.1%   295μs     0.00B  0.00%    0.00B
 ctprod       442    115ms  46.9%   260μs   6.91KiB  100%     16.0B
 ──────────────────────────────────────────────────────────────────

If I change Aᵀ = A' to Aᵀ = transpose(A) in crmr2.jl, I get the same:

julia> op = TimedLinearOperator(A);
julia> crmr2(op, b);
julia> op.timer
──────────────────────────────────────────────────────────────────
                           Time                   Allocations
                   ──────────────────────   ───────────────────────
 Tot / % measured:      88.8s / 0.29%           18.3MiB / 0.04%

 Section   ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────
 prod         441    135ms  52.8%   307μs     0.00B  0.00%    0.00B
 tprod        442    121ms  47.2%   273μs   6.91KiB  100%     16.0B
 ──────────────────────────────────────────────────────────────────

I do observe the same extra allocations as you do with @allocated and also with @btime. I'm not sure where they come from.

[amontoison]

I rebased my pull request. We should be able to use Krylov.jl with the last release of LinearOperators.jl (v0.7.1) with this commit.

[abelsiqueira]

Use LinearOperators = "0.7.1, 1", so it has an upper limit also.

[dpo]

🤞

[coveralls]

Coverage increased (+0.1%) to 96.371% when pulling ace45b8 on amontoison:transpose_product into e5aca02 on JuliaSmoothOptimizers:master.

[codecov]

Codecov Report

Merging #151 into master will increase coverage by 0.01%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master     #151      +/-   ##
==========================================
+ Coverage   97.15%   97.16%   +0.01%     
==========================================
  Files          26       26              
  Lines        2146     2155       +9     
==========================================
+ Hits         2085     2094       +9     
  Misses         61       61

Impacted Files	Coverage Δ
src/minres_qlp.jl	100% <ø> (ø)	⬆️
src/cg.jl	100% <ø> (ø)	⬆️
src/dqgmres.jl	100% <ø> (ø)	⬆️
src/symmlq.jl	99.31% <ø> (ø)	⬆️
src/minres.jl	100% <ø> (ø)	⬆️
src/diom.jl	100% <ø> (ø)	⬆️
src/cr.jl	59.54% <ø> (-1.2%)	⬇️
src/lslq.jl	100% <100%> (ø)	⬆️
src/cgs.jl	100% <100%> (ø)	⬆️
src/usymqr.jl	100% <100%> (ø)	⬆️
... and 13 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 e5aca02...ace45b8. Read the comment docs.

[amontoison]

Wow, Travis / Appveyor tests time has been reduced by a factor 4 with our recent updates of LinearOperators and this commit !

[dpo]

The type annotations must be doing something right.

[dpo]

🤘