# LinearAlgebra.BLAS does not accept non-dense Julia Matrix type?

**URL:** <https://discourse.julialang.org/t/linearalgebra-blas-does-not-accept-non-dense-julia-matrix-type/59921>\
**Category:** Performance\
**Tags:** blas, linearalgebra\
**Created:** [April 24, 2021, 8:37am UTC](https://discourse.julialang.org/t/linearalgebra-blas-does-not-accept-non-dense-julia-matrix-type/59921 "2021-04-24T08:37:36Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![JianghuiDu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghuidu/32/8028_2.png) [@JianghuiDu](https://discourse.julialang.org/u/JianghuiDu)\
**Post date:** [April 24, 2021, 8:37am UTC](https://discourse.julialang.org/t/linearalgebra-blas-does-not-accept-non-dense-julia-matrix-type/59921/1 "2021-04-24T08:37:37Z")

</div>

```julia
LinearAlgebra.BLAS.gbmv!('N', 100, 1, 1, 1.0, A_POC, ones(100), 0.0, ones(100))

```

throws error

```julia
julia> LinearAlgebra.BLAS.gbmv!('N', 100, 1, 1, 1.0, (A_POC), ones(100), 0.0, on
es(100))
ERROR: LoadError: MethodError: no method matching strides(::Tridiagonal{Float64,
 Vector{Float64}})
Closest candidates are:
  strides(::LoopVectorization.LowDimArray) at C:\Users\Jianghui\.julia\packages\
LoopVectorization\O8WW6\src\broadcast.jl:10
  strides(::SubArray) at subarray.jl:350
  strides(::ArrayInterface.AbstractArray2) at C:\Users\Jianghui\.julia\packages\
ArrayInterface\VEFPY\src\ArrayInterface.jl:810
  ...
Stacktrace:
 [1] stride(A::Tridiagonal{Float64, Vector{Float64}}, k::Int64)
   @ Base .\abstractarray.jl:501
 [2] stride1(x::Tridiagonal{Float64, Vector{Float64}})
   @ LinearAlgebra ~\AppData\Local\Programs\Julia-1.6.0\share\julia\stdlib\v1.6\
LinearAlgebra\src\LinearAlgebra.jl:197
 [3] _chkstride1
   @ ~\AppData\Local\Programs\Julia-1.6.0\share\julia\stdlib\v1.6\LinearAlgebra\
src\LinearAlgebra.jl:203 [inlined]
 [4] chkstride1
   @ ~\AppData\Local\Programs\Julia-1.6.0\share\julia\stdlib\v1.6\LinearAlgebra\
src\LinearAlgebra.jl:201 [inlined]
 [5] gbmv!(trans::Char, m::Int64, kl::Int64, ku::Int64, alpha::Float64, A::Tridi
agonal{Float64, Vector{Float64}}, x::Vector{Float64}, beta::Float64, y::Vector{F
loat64})
   @ LinearAlgebra.BLAS ~\AppData\Local\Programs\Julia-1.6.0\share\julia\stdlib\
v1.6\LinearAlgebra\src\blas.jl:783
 [6] top-level scope
   @ ~\Documents\Manuscripts\SBB\v12\main_fvcf.jl:29
in expression starting at c:\Users\Jianghui\Documents\Manuscripts\SBB\v12\main_f
vcf.jl:29

```

when A is a tridiagonal matrix

```julia
julia> A_POC
100×100 Tridiagonal{Float64, Vector{Float64}}:
 -4736.81 2509.53 ⋅ … ⋅ ⋅ ⋅
  4286.01 -6265.57 1979.56 ⋅ ⋅ ⋅
      ⋅ 3614.68 -5169.42 ⋅ ⋅ ⋅
      ⋅ ⋅ 3056.74 ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
     ⋮ ⋱
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
      ⋅ ⋅ ⋅ 0.0 ⋅ ⋅
      ⋅ ⋅ ⋅ -0.0712289 0.0 ⋅
      ⋅ ⋅ ⋅ 0.061917 -0.061917 0.
      ⋅ ⋅ ⋅ ⋅ 0.0537804 -0.0537804

```

but it works after converting A\_POC to a dense Matrix

```julia
LinearAlgebra.BLAS.gbmv!('N', 100, 1, 1, 1.0, Array(A_POC), ones(100), 0.0, ones(100))

```

Does this mean the `BLAS` functions only accept dense matrices?

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [April 24, 2021, 9:01am UTC](https://discourse.julialang.org/t/linearalgebra-blas-does-not-accept-non-dense-julia-matrix-type/59921/2 "2021-04-24T09:01:08Z")

</div>

If you directly call the LAPACK functions, the answer is likely yes.

Is there a good reason why you’re doing that though? There is a five-argument `mul!` which is what you probably want to use instead. It automatically should dispatch to the best available method for the given input types (including `Tridiagonal`):

```julia
mul!(C, A, B, α, β) -> C

  Combined inplace matrix-matrix or matrix-vector multiply-add A B α + C β. The result is stored
  in C by overwriting it. Note that C must not be aliased with either A or B.

```

Demonstration:

```julia
julia> A = rand(100,100); b = rand(100); c = rand(100); α = rand(); β = rand();

julia> T = Tridiagonal(A);

julia> @which mul!(c, T, b, α, β) # dispatches to `mul!` in bidiag.jl
mul!(C::AbstractVector{T} where T, A::Union{Bidiagonal, SymTridiagonal, Tridiagonal}, B::AbstractVector{T} where T, alpha::Number, beta::Number) in LinearAlgebra at /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/stdlib/v1.6/LinearAlgebra/src/bidiag.jl:398

julia> @which mul!(c, A, b, α, β) # dispatches to `mul!` in generic matmul.jl
mul!(y::StridedVector{T}, A::StridedVecOrMat{T}, x::StridedVector{T}, alpha::Number, beta::Number) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} in LinearAlgebra at /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:66

julia> @btime mul!($c, $A, $b, $α, $β);
  2.322 μs (0 allocations: 0 bytes)

julia> @btime mul!($c, $T, $b, $α, $β);
  155.034 ns (1 allocation: 32 bytes)

```

I have to note, though, that for the dense matrix, an explicit call into `LinearAlgebra.BLAS.gbmv!` seems faster (don’t know why that is):

```julia
julia> @btime LinearAlgebra.BLAS.gbmv!('N', 100, 1, 1, $α, $A, $b, $β, $c);
  1.044 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [April 24, 2021, 9:09am UTC](https://discourse.julialang.org/t/linearalgebra-blas-does-not-accept-non-dense-julia-matrix-type/59921/3 "2021-04-24T09:09:14Z")

</div>

Inspecting the `mul!` method in the generic case, it doesn’t dispatch to `gbmv` but to `gemv`

```julia
@inline mul!(y::StridedVector{T}, A::StridedVecOrMat{T}, x::StridedVector{T},
             alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemv!(y, 'N', A, x, alpha, beta)

```

But some LAPACK expert can probably comment on whether this is good (for some reason) or a “bug”.

UPDATE: Ah, of course, we generally can not assume that `kl == ku == 1`. So the `gbmv` method has more information than `gemv` to operate on. Importantly, the Julia implementation for the `Tridiagonal` case is most efficient for this matrix structure. So, as long as one indicates the structure through the type, one “doesn’t have to care” and automatically gets the best speed (as it should be 😄 ).

---

<div class="post-metadata">

**Author:** ![JianghuiDu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghuidu/32/8028_2.png) [@JianghuiDu](https://discourse.julialang.org/u/JianghuiDu)\
**Post date:** [April 24, 2021, 6:55pm UTC](https://discourse.julialang.org/t/linearalgebra-blas-does-not-accept-non-dense-julia-matrix-type/59921/4 "2021-04-24T18:55:37Z")

</div>

I was just curious if the BLAS version is faster than `mul!` I suppose `gbmv` must exist for a reason of better matrix vector multiplication than `gemv`?