Skip to content

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

@goerz

Description

@goerz

Maybe this is a case of "you're using it wrong", but the following MWE runs fine on Julia 1.8

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 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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions