In BLAS overrides, remove restriction to StridedMatrix

[dlfivefifty]

Definitions such as gemv! (https://github.com/JuliaLang/julia/blob/master/base/linalg/blas.jl#L553) restrict to StridedVecOrMat. This prevents user-defined types from calling these functions.

I think this restriction is unnecessary as a user AbstractMatrix which overrides the (as-of-yet not explicitly defined) strided matrix "interface" of unsafe_convert(::Type{Ptr{T}},A), and stride should be able to call gemv!.

To overcome this restriction, I end up having to copy and paste Base's code to make my own gemv!:

https://github.com/JuliaMatrices/BlockBandedMatrices.jl/blob/master/src/lapack.jl#L19

This is not too burdensome but seems unnecessary.

[andreasnoack]

Why is your custom type not a DenseArray?

[dlfivefifty]

Because it’s a SubArray (representing a block of a block matrix)

[timholy]

That method you defined is going to give you wrong answers for quite a few matrices. For example, try it with view(rand(5,5), [1,2,5], :).

[dlfivefifty]

One could a check that stride(A,1) == 1 .

EDIT: I think I mean chkstride1(A)

[timholy]

Unfortunately there's some confusion about what stride means, and it has fallbacks that it probably shouldn't have:

A = view(rand(5,5), [1,2,5], :)
b = vec(A)

julia> stride(b, 1)
1

[timholy]

#17812

[timholy]

@dlfivefifty, want to tackle that for 1.0? 😄

[dlfivefifty]

Sure, I'm happy to put together a pull request. When's the deadline?

[StefanKarpinski]

Last week, but we're still wrapping some things up. This is a fairly non-disruptive change as I understand it though, so it could probably sneak into the alpha period.

[dlfivefifty]

OK. While I'm at it I might also include some of the missing routines (sbgst) for banded matrix routines, which are currently in

https://github.com/JuliaMatrices/BandedMatrices.jl/blob/master/src/lapack.jl

[dlfivefifty]

One issue here is with axpy! and axpby!, which exist for both generic arrays and BLAS arrays. I'd almost say deprecate axpy! for generic arrays in favour of y .= a.*x.

[andreasnoack]

Made a comment about this in JuliaLinearAlgebra/IterativeSolvers.jl#184. For new array types, it is much more demanding to support general broadcasting so I think we should keep axpy! and friends.

[dlfivefifty]

Yes I did a quick timing of axpy! with sparse vectors and it's a factor of 10 slower to use broadcasting.

[Sacha0]

Yes I did a quick timing of axpy! with sparse vectors and it's a factor of 10 slower to use broadcasting.

Modulo benchmark details, broadcast should not be so much slower; perhaps something has regressed?