search

JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
399 stars 105 forks source link

Can `adjoint(A)` be placed outside of the loop? #240

Closed vincentcp closed 5 years ago

vincentcp commented 5 years ago

I have a linear object A that allocates some memory when creating the adjoint. When I apply lsqr, this adjoint is created in every iteration.

Can adjoint(A) be placed outside of the loop?

andreasnoack commented 5 years ago

Yes. That should be possible. The best way to make this change happen is if you could prepare the PR. It should be a simple one.

mschauer commented 5 years ago

The author's of lsqr give a helpful comment:

            # Note that the following three lines are a band aid for a GEMM: X: C := αA'B + βC.
            # This is already supported in mul! for sparse and distributed matrices, but not yet dense
            mul!(tmpn, adjoint(A), u)
            v .= -beta .* v .+ tmpn
andreasnoack commented 5 years ago

The author's of lsqr give a helpful comment:

I think I'm missing something. Why would that affect moving the call toadjoint outside the loop and save the result as e.g. adjA?

vincentcp commented 5 years ago

If adjoint(A) is an expensive operation why do I do need to do it every iteration while its result stays the same?

mschauer commented 5 years ago

One can certainly take the adjoint out of the loop. The comment just points at the underlying problem. You (@andreasnoack) wrote this 4 years ago. Maybe it is fixed by now?

andreasnoack commented 5 years ago

Maybe it is fixed by now?

Unfortunate not, see https://github.com/JuliaLang/julia/issues/23919