Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modified CG to use axpy!/axpby! directly instead of broadcasting #334

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

vlc1
Copy link

@vlc1 vlc1 commented Mar 17, 2023

When using IterativeSolvers.jl on custom vectors and matrices, I find it easier to add methods to generic BLAS routines (scal!, mul!, axpy!...) than customizing the broadcasting machinery on said types.

This is a very small modification that simply requires changing lines like

it.r .-= alpha .* it.c

to

axpy!(-alpha, it.c, it.r)

and so on as illustrated in the cg.jl file in this PR.

If this is a modification that is acceptable, I'll work my way through the other solvers (bicgstabl...).

@coveralls
Copy link

Pull Request Test Coverage Report for Build 4446585905

  • 7 of 7 (100.0%) changed or added relevant lines in 1 file are covered.
  • No unchanged relevant lines lost coverage.
  • Overall coverage decreased (-0.0%) to 97.331%

Totals Coverage Status
Change from base Build 2991758705: -0.0%
Covered Lines: 1933
Relevant Lines: 1986

💛 - Coveralls

@haampie
Copy link
Member

haampie commented Mar 17, 2023

I believe we've gone back and forth on this. Long ago there was overhead to broadcasting, so the BLAS routines were called. Then Julia improved and folks wanted broadcasting because it's more generic. I think I'm happy with the current broadcasting stuff, also because for example CG would need an axpby somewhere, which I think is not standard BLAS.

@dkarrasch
Copy link
Member

These functions are now owned by LinearAlgebra (and not LinearAlgebra.BLAS) and admit generic fallbacks, likely using broadcasting but I don't remember. Another option could be to introduce _axp[b]y!, that use broadcasting generically but may be overloaded with a custom axpy!-like implementation.

@vlc1
Copy link
Author

vlc1 commented Mar 17, 2023

These functions are now owned by LinearAlgebra (and not LinearAlgebra.BLAS) and admit generic fallbacks, likely using broadcasting but I don't remember. Another option could be to introduce _axp[b]y!, that use broadcasting generically but may be overloaded with a custom axpy!-like implementation.

Fair enough, I might have to dig into how to customize broadcasting for my own datatypes...

These functions are now owned by LinearAlgebra (and not LinearAlgebra.BLAS) and admit generic fallbacks, likely using broadcasting but I don't remember. Another option could be to introduce _axp[b]y!, that use broadcasting generically but may be overloaded with a custom axpy!-like implementation.

That's pretty much what Krylov.jl does with some macros (@kaxpy...).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants