Rewrite mul! to dispatch based on memory layout, not matrix type

NEWS.md

@@ -650,6 +650,11 @@ Library improvements
     keyword arguments. ([#26156], [#26670])
 
 
+    * `Base.MemoryLayout(A::AbstractArray)` introduced to specify the layout in
+      memory for an array `A`, by returning `DenseColumnMajor()`, `ColumnMajor()`,
+      `DenseRowMajor()`, `RowMajor()` or other specialized memory layouts. The return
+      value is used in `mul!` to dispatch to the correct BLAS or LAPACK routines. ([#25558])
+
 Compiler/Runtime improvements
 -----------------------------
 

base/abstractarray.jl

@@ -345,6 +345,138 @@ size_to_strides(s, d) = (s,)
 size_to_strides(s) = ()
 
 
+abstract type MemoryLayout end
+struct UnknownLayout <: MemoryLayout end
+abstract type AbstractStridedLayout <: MemoryLayout end
+abstract type AbstractIncreasingStrides <: AbstractStridedLayout end
+abstract type AbstractColumnMajor <: AbstractIncreasingStrides end
+struct DenseColumnMajor <: AbstractColumnMajor end
+struct ColumnMajor <: AbstractColumnMajor end
+struct IncreasingStrides <: AbstractIncreasingStrides end
+abstract type AbstractDecreasingStrides <: AbstractStridedLayout end
+abstract type AbstractRowMajor <: AbstractDecreasingStrides end
+struct DenseRowMajor <: AbstractRowMajor end
+struct RowMajor <: AbstractRowMajor end
+struct DecreasingStrides <: AbstractIncreasingStrides end
+struct StridedLayout <: AbstractStridedLayout end
+
+"""
+    UnknownLayout()
+
+is returned by `MemoryLayout(A)` if it is unknown how the entries of an array `A`
+are stored in memory.
+"""
+UnknownLayout
+
+"""
+    AbstractStridedLayout
+
+is an abstract type whose subtypes are returned by `MemoryLayout(A)`
+if an array `A` has storage laid out at regular offsets in memory,
+and which can therefore be passed to external C and Fortran functions expecting
+this memory layout.
+
+Julia's internal linear algebra machinery will automatically (and invisibly)
+dispatch to BLAS and LAPACK routines if the memory layout is BLAS compatible and
+the element type is a `Float32`, `Float64`, `ComplexF32`, or `ComplexF64`.
+In this case, one must implement the strided array interface, which requires
+overrides of `strides(A::MyMatrix)` and `unknown_convert(::Type{Ptr{T}}, A::MyMatrix)`.
+"""
+AbstractStridedLayout
+
+"""
+    DenseColumnMajor()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage in memory
+equivalent to an `Array`, so that `stride(A,1) == 1` and
+`stride(A,i) ≡ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`. In particular,
+if `A` is a matrix then `strides(A) == `(1, size(A,1))`.
+
+Arrays with `DenseColumnMajor` memory layout must conform to the `DenseArray` interface.
+"""
+DenseColumnMajor
+
+"""
+    ColumnMajor()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage in memory
+as a column major array, so that `stride(A,1) == 1` and
+`stride(A,i) ≥ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`.
+
+Arrays with `ColumnMajor` memory layout must conform to the `DenseArray` interface.
+"""
+ColumnMajor
+
+"""
+    IncreasingStrides()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage in memory
+as a strided array with  increasing strides, so that `stride(A,1) ≥ 1` and
+`stride(A,i) ≥ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`.
+"""
+IncreasingStrides
+
+"""
+    DenseRowMajor()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage in memory
+as a row major array with dense entries, so that `stride(A,ndims(A)) == 1` and
+`stride(A,i) ≡ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`. In particular,
+if `A` is a matrix then `strides(A) == `(size(A,2), 1)`.
+"""
+DenseRowMajor
+
+"""
+    RowMajor()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage in memory
+as a row major array, so that `stride(A,ndims(A)) == 1` and
+stride(A,i) ≥ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`.
+
+If `A` is a matrix  with `RowMajor` memory layout, then
+`transpose(A)` should return a matrix whose layout is `ColumnMajor`.
+"""
+RowMajor
+
+"""
+    DecreasingStrides()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage in memory
+as a strided array with decreasing strides, so that `stride(A,ndims(A)) ≥ 1` and
+stride(A,i) ≥ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`.
+"""
+DecreasingStrides
+
+"""
+    StridedLayout()
+
+is returned by `MemoryLayout(A)` if an array `A` has storage laid out at regular
+offsets in memory. `Array`s with `StridedLayout` must conform to the `DenseArray` interface.
+"""
+StridedLayout
+
+"""
+    MemoryLayout(A)
+
+specifies the layout in memory for an array `A`. When
+you define a new `AbstractArray` type, you can choose to override
+`MemoryLayout` to indicate how an array is stored in memory.
+For example, if your matrix is column major with `stride(A,2) == size(A,1)`,
+then override as follows:
+
+    Base.MemoryLayout(::MyMatrix) = Base.DenseColumnMajor()
+
+The default is `Base.UnknownLayout()` to indicate that the layout
+in memory is unknown.
+
+Julia's internal linear algebra machinery will automatically (and invisibly)
+dispatch to BLAS and LAPACK routines if the memory layout is compatible.
+"""
+MemoryLayout(A::AbstractArray{T}) where T = UnknownLayout()
+MemoryLayout(A::DenseArray{T}) where T = DenseColumnMajor()
+
+
+
 function isassigned(a::AbstractArray, i::Int...)
     try
         a[i...]
@@ -379,6 +511,50 @@ function trailingsize(inds::Indices)
     prod(map(unsafe_length, inds))
 end
 
+## Traits for array types ##
+
+abstract type IndexStyle end
+struct IndexLinear <: IndexStyle end
+struct IndexCartesian <: IndexStyle end
+
+"""
+    IndexStyle(A)
+    IndexStyle(typeof(A))
+
+`IndexStyle` specifies the "native indexing style" for array `A`. When
+you define a new `AbstractArray` type, you can choose to implement
+either linear indexing or cartesian indexing.  If you decide to
+implement linear indexing, then you must set this trait for your array
+type:
+
+    Base.IndexStyle(::Type{<:MyArray}) = IndexLinear()
+
+The default is `IndexCartesian()`.
+
+Julia's internal indexing machinery will automatically (and invisibly)
+convert all indexing operations into the preferred style. This allows users
+to access elements of your array using any indexing style, even when explicit
+methods have not been provided.
+
+If you define both styles of indexing for your `AbstractArray`, this
+trait can be used to select the most performant indexing style. Some
+methods check this trait on their inputs, and dispatch to different
+algorithms depending on the most efficient access pattern. In
+particular, [`eachindex`](@ref) creates an iterator whose type depends
+on the setting of this trait.
+"""
+IndexStyle(A::AbstractArray) = IndexStyle(typeof(A))
+IndexStyle(::Type{Union{}}) = IndexLinear()
+IndexStyle(::Type{<:AbstractArray}) = IndexCartesian()
+IndexStyle(::Type{<:Array}) = IndexLinear()
+IndexStyle(::Type{<:AbstractRange}) = IndexLinear()
+
+IndexStyle(A::AbstractArray, B::AbstractArray) = IndexStyle(IndexStyle(A), IndexStyle(B))
+IndexStyle(A::AbstractArray, B::AbstractArray...) = IndexStyle(IndexStyle(A), IndexStyle(B...))
+IndexStyle(::IndexLinear, ::IndexLinear) = IndexLinear()
+IndexStyle(::IndexStyle, ::IndexStyle) = IndexCartesian()
+
+
 ## Bounds checking ##
 
 # The overall hierarchy is

base/reinterpretarray.jl

@@ -44,6 +44,12 @@ function size(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
     tuple(size1, tail(psize)...)
 end
 
+
+MemoryLayout(A::ReinterpretArray) = reinterpretedmemorylayout(MemoryLayout(parent(A)))
+reinterpretedmemorylayout(::MemoryLayout) = UnknownLayout()
+reinterpretedmemorylayout(::DenseColumnMajor) = DenseColumnMajor()
+
+
 elsize(::Type{<:ReinterpretArray{T}}) where {T} = sizeof(T)
 unsafe_convert(::Type{Ptr{T}}, a::ReinterpretArray{T,N,S} where N) where {T,S} = Ptr{T}(unsafe_convert(Ptr{S},a.parent))
 

base/reshapedarray.jl

@@ -251,4 +251,9 @@ setindex!(A::ReshapedRange, val, index::ReshapedIndex) = _rs_setindex!_err()
 
 @noinline _rs_setindex!_err() = error("indexed assignment fails for a reshaped range; consider calling collect")
 
+
+MemoryLayout(A::ReshapedArray) = reshapedmemorylayout(MemoryLayout(parent(A)))
+reshapedmemorylayout(::MemoryLayout) = UnknownLayout()
+reshapedmemorylayout(::DenseColumnMajor) = DenseColumnMajor()
+
 unsafe_convert(::Type{Ptr{T}}, a::ReshapedArray{T}) where {T} = unsafe_convert(Ptr{T}, parent(a))

base/subarray.jl

@@ -320,6 +320,82 @@ find_extended_inds(::ScalarIndex, I...) = (@_inline_meta; find_extended_inds(I..
 find_extended_inds(i1, I...) = (@_inline_meta; (i1, find_extended_inds(I...)...))
 find_extended_inds() = ()
 
+MemoryLayout(A::SubArray) = subarraylayout(MemoryLayout(parent(A)), parentindices(A))
+subarraylayout(_1, _2) = UnknownLayout()
+subarraylayout(_1, _2, _3)= UnknownLayout()
+subarraylayout(::DenseColumnMajor, ::Tuple{I}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
+    DenseColumnMajor()  # A[:] is DenseColumnMajor if A is DenseColumnMajor
+subarraylayout(ml::AbstractColumnMajor, inds) = _column_subarraylayout1(ml, inds)
+subarraylayout(::AbstractRowMajor, ::Tuple{I}) where I =
+    UnknownLayout()  # A[:] does not have any structure if A is AbstractRowMajor
+subarraylayout(ml::AbstractRowMajor, inds) = _row_subarraylayout1(ml, reverse(inds))
+subarraylayout(ml::AbstractStridedLayout, inds) = _strided_subarraylayout(ml, inds)
+
+_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:Union{Int,AbstractCartesianIndex} =
+    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
+_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:Slice =
+    DenseColumnMajor() # view(A,:,1,2) is a DenseColumnMajor vector
+_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:AbstractUnitRange{Int} =
+    DenseColumnMajor() # view(A,1:3,1,2) is a DenseColumnMajor vector
+_column_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:Union{Int,AbstractCartesianIndex} =
+    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
+_column_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:AbstractUnitRange{Int} =
+    DenseColumnMajor() # view(A,1:3,1,2) is a DenseColumnMajor vector
+_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Any}}) where I<:Slice =
+    _column_subarraylayout(DenseColumnMajor(), DenseColumnMajor(), tail(inds))
+_column_subarraylayout1(par::DenseColumnMajor, inds::Tuple{I,Vararg{Any}}) where I<:AbstractUnitRange{Int} =
+    _column_subarraylayout(par, ColumnMajor(), tail(inds))
+_column_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:AbstractUnitRange{Int} =
+    _column_subarraylayout(par, ColumnMajor(), tail(inds))
+_column_subarraylayout1(par::DenseColumnMajor, inds::Tuple{I,Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
+    _column_subarraylayout(par, StridedLayout(), tail(inds))
+_column_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
+    _column_subarraylayout(par, StridedLayout(), tail(inds))
+_column_subarraylayout1(par, inds) = UnknownLayout()
+_column_subarraylayout(par, ret, ::Tuple{}) = ret
+_column_subarraylayout(par, ret, ::Tuple{I}) where I = UnknownLayout()
+_column_subarraylayout(::DenseColumnMajor, ::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
+    DenseColumnMajor() # A[:,1:3,1,2] is DenseColumnMajor if A is DenseColumnMajor
+_column_subarraylayout(par::DenseColumnMajor, ::DenseColumnMajor, inds::Tuple{I, Vararg{Int}}) where I<:Slice =
+    DenseColumnMajor()
+_column_subarraylayout(par::DenseColumnMajor, ::DenseColumnMajor, inds::Tuple{I, Vararg{Any}}) where I<:Slice =
+    _column_subarraylayout(par, DenseColumnMajor(), tail(inds))
+_column_subarraylayout(par, ::AbstractColumnMajor, inds::Tuple{I, Vararg{Any}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
+    _column_subarraylayout(par, ColumnMajor(), tail(inds))
+_column_subarraylayout(par, ::AbstractStridedLayout, inds::Tuple{I, Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
+    _column_subarraylayout(par, StridedLayout(), tail(inds))
+
+_row_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:Union{Int,AbstractCartesianIndex} =
+    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
+_row_subarraylayout1(::DenseRowMajor, inds::Tuple{I,Vararg{Int}}) where I<:Slice =
+    DenseColumnMajor() # view(A,1,2,:) is a DenseColumnMajor vector
+_row_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:AbstractUnitRange{Int} =
+    DenseColumnMajor() # view(A,1,2,1:3) is a DenseColumnMajor vector
+_row_subarraylayout1(::DenseRowMajor, inds::Tuple{I,Vararg{Any}}) where I<:Slice =
+    _row_subarraylayout(DenseRowMajor(), DenseRowMajor(), tail(inds))
+_row_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:AbstractUnitRange{Int} =
+    _row_subarraylayout(par, RowMajor(), tail(inds))
+_row_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
+    _row_subarraylayout(par, StridedLayout(), tail(inds))
+_row_subarraylayout1(par, inds) = UnknownLayout()
+_row_subarraylayout(par, ret, ::Tuple{}) = ret
+_row_subarraylayout(par, ret, ::Tuple{I}) where I = UnknownLayout()
+_row_subarraylayout(::DenseRowMajor, ::DenseRowMajor, inds::Tuple{I,Vararg{Int}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
+    DenseRowMajor() # A[1,2,1:3,:] is DenseRowMajor if A is DenseRowMajor
+_row_subarraylayout(par::DenseRowMajor, ::DenseRowMajor, inds::Tuple{I, Vararg{Int}}) where I<:Slice =
+    DenseRowMajor()
+_row_subarraylayout(par::DenseRowMajor, ::DenseRowMajor, inds::Tuple{I, Vararg{Any}}) where I<:Slice =
+    _row_subarraylayout(par, DenseRowMajor(), tail(inds))
+_row_subarraylayout(par::AbstractRowMajor, ::AbstractRowMajor, inds::Tuple{I, Vararg{Any}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
+    _row_subarraylayout(par, RowMajor(), tail(inds))
+_row_subarraylayout(par::AbstractRowMajor, ::AbstractStridedLayout, inds::Tuple{I, Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
+    _row_subarraylayout(par, StridedLayout(), tail(inds))
+
+_strided_subarraylayout(par, inds) = UnknownLayout()
+_strided_subarraylayout(par, ::Tuple{}) = StridedLayout()
+_strided_subarraylayout(par, inds::Tuple{I, Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
+    _strided_subarraylayout(par, tail(inds))
+
 unsafe_convert(::Type{Ptr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{RangeIndex}}}) where {T,N,P} =
     unsafe_convert(Ptr{T}, V.parent) + (first_index(V)-1)*sizeof(T)
 

doc/src/manual/interfaces.md

@@ -227,6 +227,7 @@ ourselves, we can officially define it as a subtype of an [`AbstractArray`](@ref
 | `similar(A, ::Type{S})`                         | `similar(A, S, size(A))`               | Return a mutable array with the same shape and the specified element type             |
 | `similar(A, dims::NTuple{Int})`                 | `similar(A, eltype(A), dims)`          | Return a mutable array with the same element type and size *dims*                     |
 | `similar(A, ::Type{S}, dims::NTuple{Int})`      | `Array{S}(undef, dims)`               | Return a mutable array with the specified element type and size                       |
+| `Base.MemoryLayout(A)`                            | `Base.UnknownLayout()`                     | Return a subtype of `Base.MemoryLayout`, e.g. `Base.ColumnMajor()`      |
 | **Non-traditional indices**                     | **Default definition**                 | **Brief description**                                                                 |
 | `axes(A)`                                    | `map(OneTo, size(A))`                  | Return the `AbstractUnitRange` of valid indices                                       |
 | `Base.similar(A, ::Type{S}, inds::NTuple{Ind})` | `similar(A, S, Base.to_shape(inds))`   | Return a mutable array with the specified indices `inds` (see below)                  |
@@ -401,6 +402,7 @@ perhaps range-types `Ind` of your own design. For more information, see
 |:----------------------------------------------- |:-------------------------------------- |:------------------------------------------------------------------------------------- |
 | `strides(A)`                             |                                        | Return the distance in memory (in number of elements) between adjacent elements in each dimension as a tuple. If `A` is an `AbstractArray{T,0}`, this should return an empty tuple.    |
 | `Base.unsafe_convert(::Type{Ptr{T}}, A)`        |                                        | Return the native address of an array.                                            |
+| `Base.MemoryLayout(A)`                 |                                        | Return a subtype of `Base.AbstractStridedLayout`,  such as `Base.DenseColumnMajor()`,`Base.DenseRowMajor()`, or `Base.StridedLayout()`.
 | **Optional methods**                            | **Default definition**                 | **Brief description**                                                                 |
 | `stride(A, i::Int)`                             |     `strides(A)[i]`                                   | Return the distance in memory (in number of elements) between adjacent elements in dimension k.    |
 

stdlib/IterativeEigensolvers/test/runtests.jl

@@ -101,7 +101,7 @@ using Test, LinearAlgebra, SparseArrays, Random
         k = 3
         A = randn(n,n); A = A'A
         B = randn(n,k); B = B*B'
-        @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k])
+        @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(filter(isfinite, eigvals(A, B)))
     end
 end
 

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -16,11 +16,14 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     getproperty, imag, inv, isapprox, isone, IndexStyle, kron, length, log, map, ndims,
     oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round, sec, sech,
     setindex!, show, similar, sin, sincos, sinh, size, size_to_strides, sqrt, StridedReinterpretArray,
-    StridedReshapedArray, strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec
+    StridedReshapedArray, ReshapedArray, ReinterpretArray, strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec,
+    MemoryLayout, UnknownLayout, AbstractStridedLayout, AbstractRowMajor, AbstractColumnMajor, DenseRowMajor, DenseColumnMajor,
+    ColumnMajor, RowMajor, StridedLayout
 using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof,
     @propagate_inbounds, @pure, reduce, typed_vcat
 using Base.Broadcast: Broadcasted
 
+
 # We use `_length` because of non-1 indices; releases after julia 0.5
 # can go back to `length`. `_length(A)` is equivalent to `length(LinearIndices(A))`.
 
@@ -201,8 +204,9 @@ julia> LinearAlgebra.stride1(B)
 ```
 """
 stride1(x) = stride(x,1)
-stride1(x::Array) = 1
-stride1(x::DenseArray) = stride(x, 1)::Int
+stride1(x::AbstractArray) = _stride1(x, MemoryLayout(x))
+_stride1(x, _) = stride(x, 1)::Int
+_stride1(x, ::AbstractColumnMajor) = 1
 
 @inline chkstride1(A...) = _chkstride1(true, A...)
 @noinline _chkstride1(ok::Bool) = ok || error("matrix does not have contiguous columns")