What's the difference between LinearAlgebra.axpy! and BLAS.axpy!, and which one should packages extend?

jishnub · December 13, 2022, 5:15am

BLAS sounds too low level, so intuitively I would think it’s LinearAlgebra.axpy! that packages should extend. However, what is the intended difference between the two? And when should one use one over the other?

julia> VERSION
v"1.9.0-alpha1"

julia> LinearAlgebra.axpy! === BLAS.axpy!
false

github.com/JuliaLang/julia

`LinearAlgebra.BLAS.axpy!` breaks for sparse matrices in Julia 1.9

opened 02:08AM - 07 Sep 22 UTC

closed 10:06AM - 12 Oct 22 UTC

goerz

linear algebra

Maybe this is a case of "you're using it wrong", but the following MWE runs fine… on Julia 1.8 ~~~julia using SparseArrays using LinearAlgebra.BLAS: axpy! # breaks # using LinearAlgebra: axpy! # works N = 5 sparsity = 0.3 H0 = sprand(ComplexF64, N, N, sparsity) H1 = sprand(ComplexF64, N, N, sparsity) axpy!(1.0, H0, H1) ~~~ On Julia 1.9 (nightly of 2022-09-06), I get ~~~ ERROR: LoadError: MethodError: no method matching strides(::SparseMatrixCSC{ComplexF64, Int64}) Closest candidates are: strides(::SubArray) @ Base subarray.jl:360 strides(::Union{Base.ReinterpretArray{T, N, S, A, IsReshaped} where {T, N, A<:Union{SubArray{T, N, A, I, true} where {T, N, A<:DenseArray, I<:Union{Tuple{Vararg{Real}}, Tuple{AbstractUnitRange, Vararg{Any}}}}, DenseArray}, IsReshaped, S}, Base.ReshapedArray{T, N, A} where {T, N, A<:Union{Base.ReinterpretArray{T, N, S, A, IsReshaped} where {T, N, A<:Union{SubArray{T, N, A, I, true} where {T, N, A<:DenseArray, I<:Union{Tuple{Vararg{Real}}, Tuple{AbstractUnitRange, Vararg{Any}}}}, DenseArray}, IsReshaped, S}, SubArray{T, N, A, I, true} where {T, N, A<:DenseArray, I<:Union{Tuple{Vararg{Real}}, Tuple{AbstractUnitRange, Vararg{Any}}}}, DenseArray}}, DenseArray}) @ Base reinterpretarray.jl:151 strides(::Base.ReshapedArray) @ Base reshapedarray.jl:305 ... Stacktrace: [1] checkedstride(x::SparseMatrixCSC{ComplexF64, Int64}) @ LinearAlgebra.BLAS /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:178 [2] vec_pointer_stride @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:171 [inlined] [3] vec_pointer_stride @ /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:170 [inlined] [4] axpy!(alpha::Float64, x::SparseMatrixCSC{ComplexF64, Int64}, y::SparseMatrixCSC{ComplexF64, Int64}) @ LinearAlgebra.BLAS /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:515 [5] top-level scope @ ~/Desktop/DEBUG/debug2.jl:10 ~~~ The problem is that in Julia 1.9, different methods get called depending on whether I import `axpy!` from `LinearAlgebra` or from `LinearAlgebra.BLAS`. Maybe you say I should stick to the [manual](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.axpy!) and only ever use `LinearAlgebra.axpy!`, but the example right there in the manual actually also shows `BLAS.axpy!`. My assumption was that the two were interchangeable. This seems to have been the case in Julia 1.8, but not in Julia 1.9. I'd be okay with `LinearAlgebra.BLAS` only being for dense matrices and `LinearAlgebra` being for `AbstractArrays`, as long as that is very clearly documented. Also, the error I got here is not very illuminating. It took me quite a while to find out what was happening here.

makes it sound like BLAS.axpy! is meant for strided arrays. Is this the distinction?

Edit: After reading some more, it seems BLAS.axpy! will necessarily call BLAS, whereas LinearAlgebra.axpy! is more liberal in what it does. It may choose to call BLAS.axpy! if a suitable BLAS call exists.