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

base/abstractarray.jl

@@ -333,6 +333,7 @@ IndexStyle(A::AbstractArray, B::AbstractArray...) = IndexStyle(IndexStyle(A), In
 IndexStyle(::IndexLinear, ::IndexLinear) = IndexLinear()
 IndexStyle(::IndexStyle, ::IndexStyle) = IndexCartesian()
 
+
 ## Bounds checking ##
 
 # The overall hierarchy is

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(eigvals(A, B))[2:4]
     end
 end
 

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -18,7 +18,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     setindex!, show, similar, sin, sincos, sinh, size, size_to_strides, sqrt, StridedReinterpretArray,
     StridedReshapedArray, strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec
 using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof,
-    @propagate_inbounds, @pure, reduce, typed_vcat
+    @propagate_inbounds, @pure, reduce, typed_vcat, AbstractCartesianIndex, RangeIndex
 # 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))`.
 
@@ -160,6 +160,54 @@ else
     const BlasInt = Int32
 end
 
+## Traits for memory layouts ##
+
+abstract type MemoryLayout{T} end
+struct UnknownLayout{T} <: MemoryLayout{T} end
+struct StridedLayout{T} <: MemoryLayout{T} end
+struct TransposeLayout{T} <: MemoryLayout{T} end
+struct CTransposeLayout{T} <: MemoryLayout{T} end
+struct LowerTriangularLayout{trans,T} <: MemoryLayout{T} end
+struct UnitLowerTriangularLayout{trans,T} <: MemoryLayout{T} end
+struct UpperTriangularLayout{trans,T} <: MemoryLayout{T} end
+struct UnitUpperTriangularLayout{trans,T} <: MemoryLayout{T} end
+struct SymmetricLayout{T} <: MemoryLayout{T}
+    uplo::Char
+end
+struct HermitianLayout{T} <: MemoryLayout{T}
+    uplo::Char
+end
+
+"""
+    MemoryLayout(A)
+    MemoryLayout(typeof(A))
+
+`MemoryLayout` specifies the layout in memory for an array `A`. When
+you define a new `AbstractArray` type, you can choose to implement
+memory layout to indicate that an array is strided in memory. If you decide to
+implement memory layout, then you must set this trait for your array
+type:
+
+    Base.MemoryLayout(::Type{M}) where M <: MyArray{T,N} where {T,N} = Base.StridedLayout{T,N}()
+
+The default is `Base.UnknownLayout{T,N}()` 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 BLAS 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::MyArray)` and `unknown_convert(::Type{Ptr{T}}, A::MyArray)`
+"""
+MemoryLayout(A::AbstractArray{T,N}) where {T,N} = UnknownLayout{T}()
+MemoryLayout(A::Vector{T}) where {T} = StridedLayout{T}()
+MemoryLayout(A::Matrix{T}) where {T} = StridedLayout{T}()
+MemoryLayout(A::SubArray) = submemorylayout(MemoryLayout(parent(A)), parentindices(A))
+
+submemorylayout(::MemoryLayout{T}, _) where T = UnknownLayout{T}()
+submemorylayout(S::StridedLayout{T}, ::Tuple{I}) where {T,I<:Union{RangeIndex, AbstractCartesianIndex}} = S
+submemorylayout(S::StridedLayout{T}, ::NTuple{2,I}) where {T,I<:Union{RangeIndex, AbstractCartesianIndex}} = S
+
 # Check that stride of matrix/vector is 1
 # Writing like this to avoid splatting penalty when called with multiple arguments,
 # see PR 16416

stdlib/LinearAlgebra/src/adjtrans.jl

@@ -94,6 +94,8 @@ const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix}
 wrapperop(A::Adjoint) = adjoint
 wrapperop(A::Transpose) = transpose
 
+
+
 # AbstractArray interface, basic definitions
 length(A::AdjOrTrans) = length(A.parent)
 size(v::AdjOrTransAbsVec) = (1, length(v.parent))
@@ -102,6 +104,17 @@ axes(v::AdjOrTransAbsVec) = (Base.OneTo(1), axes(v.parent)...)
 axes(A::AdjOrTransAbsMat) = reverse(axes(A.parent))
 IndexStyle(::Type{<:AdjOrTransAbsVec}) = IndexLinear()
 IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian()
+
+transpose(::MemoryLayout{T}) where {T} = UnknownLayout{T}()
+transpose(::StridedLayout{T}) where {T} = TransposeLayout{T}()
+transpose(::TransposeLayout{T}) where {T} = StridedLayout{T}()
+adjoint(::MemoryLayout{T}) where {T} = UnknownLayout{T}()
+adjoint(M::MemoryLayout{T}) where {T<:Real} = transpose(M)
+adjoint(::StridedLayout{T}) where {T<:Complex} = CTransposeLayout{T}()
+adjoint(::CTransposeLayout{T}) where {T<:Complex} = StridedLayout{T}()
+MemoryLayout(::Type{A}) where A<:Adjoint{T,S} where {T,S} = adjoint(MemoryLayout(S))
+MemoryLayout(::Type{A}) where A<:Transpose{T,S} where {T,S} = transpose(MemoryLayout(S))
+
 @propagate_inbounds getindex(v::AdjOrTransAbsVec, i::Int) = wrapperop(v)(v.parent[i])
 @propagate_inbounds getindex(A::AdjOrTransAbsMat, i::Int, j::Int) = wrapperop(A)(A.parent[j, i])
 @propagate_inbounds setindex!(v::AdjOrTransAbsVec, x, i::Int) = (setindex!(v.parent, wrapperop(v)(x), i); v)

stdlib/LinearAlgebra/src/bidiag.jl

@@ -350,8 +350,7 @@ mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::BiTriSym) =
 mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractVecOrMat}, B::BiTriSym) = A_mul_B_td!(C, A, B)
 mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractVecOrMat}, B::BiTriSym) = A_mul_B_td!(C, A, B)
 mul!(C::AbstractVector,   A::BiTri,              B::AbstractVector) = A_mul_B_td!(C, A, B)
-mul!(C::AbstractMatrix,   A::BiTri,              B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
-mul!(C::AbstractVecOrMat, A::BiTri,              B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
+mul!(C::AbstractMatrix,   A::BiTri,              B::AbstractMatrix) = A_mul_B_td!(C, A, B)
 mul!(C::AbstractMatrix, A::BiTri, B::Transpose{<:Any,<:AbstractVecOrMat}) = A_mul_B_td!(C, A, B) # around bidiag line 330
 mul!(C::AbstractMatrix, A::BiTri, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A_mul_B_td!(C, A, B)
 mul!(C::AbstractVector, A::BiTri, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B)))

stdlib/LinearAlgebra/src/dense.jl

@@ -35,7 +35,7 @@ function scale!(X::Array{T}, s::Real) where T<:BlasComplex
 end
 
 
-function isone(A::StridedMatrix)
+function isone(A::AbstractMatrix)
     m, n = size(A)
     m != n && return false # only square matrices can satisfy x == one(x)
     if sizeof(A) < ISONE_CUTOFF
@@ -45,7 +45,7 @@ function isone(A::StridedMatrix)
     end
 end
 
-@inline function _isone_triacheck(A::StridedMatrix, m::Int)
+@inline function _isone_triacheck(A::AbstractMatrix, m::Int)
     @inbounds for i in 1:m, j in i:m
         if i == j
             isone(A[i,i]) || return false
@@ -57,7 +57,7 @@ end
 end
 
 # Inner loop over rows to be friendly to the CPU cache
-@inline function _isone_cachefriendly(A::StridedMatrix, m::Int)
+@inline function _isone_cachefriendly(A::AbstractMatrix, m::Int)
     @inbounds for i in 1:m, j in 1:m
         if i == j
             isone(A[i,i]) || return false
@@ -496,12 +496,19 @@ julia> exp(A)
  0.0      2.71828
 ```
 """
-exp(A::StridedMatrix{<:BlasFloat}) = exp!(copy(A))
-exp(A::StridedMatrix{<:Union{Integer,Complex{<:Integer}}}) = exp!(float.(A))
+exp(A::AbstractMatrix{<:BlasFloat}) = exp!(copy(A))
+exp(A::AbstractMatrix{<:Union{Integer,Complex{<:Integer}}}) = exp!(float.(A))
 
 ## Destructive matrix exponential using algorithm from Higham, 2008,
 ## "Functions of Matrices: Theory and Computation", SIAM
-function exp!(A::StridedMatrix{T}) where T<:BlasFloat
+exp!(A::AbstractMatrix) = _exp!(A, MemoryLayout(A))
+exp!(A::AbstractMatrix, _) =
+    throw(ArgumentError("exp! not implemented for non-BlasFloat types."))
+function exp!(A::AbstractMatrix{T}, _) where T<:BlasFloat
+    A .= exp!(Matrix{T}(A))
+    A
+end
+function _exp!(A::AbstractMatrix{T}, ::StridedLayout{T}) where T<:BlasFloat
     n = checksquare(A)
     if ishermitian(A)
         return copytri!(parent(exp(Hermitian(A))), 'U', true)
@@ -585,7 +592,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat
 end
 
 ## Swap rows i and j and columns i and j in X
-function rcswap!(i::Integer, j::Integer, X::StridedMatrix{<:Number})
+function rcswap!(i::Integer, j::Integer, X::AbstractMatrix{<:Number})
     for k = 1:size(X,1)
         X[k,i], X[k,j] = X[k,j], X[k,i]
     end
@@ -595,7 +602,7 @@ function rcswap!(i::Integer, j::Integer, X::StridedMatrix{<:Number})
 end
 
 """
-    log(A{T}::StridedMatrix{T})
+    log(A{T}::AbstractMatrix{T})
 
 If `A` has no negative real eigenvalue, compute the principal matrix logarithm of `A`, i.e.
 the unique matrix ``X`` such that ``e^X = A`` and ``-\\pi < Im(\\lambda) < \\pi`` for all
@@ -625,7 +632,7 @@ julia> log(A)
  0.0  1.0
 ```
 """
-function log(A::StridedMatrix)
+function log(A::AbstractMatrix)
     # If possible, use diagonalization
     if ishermitian(A)
         logHermA = log(Hermitian(A))
@@ -684,7 +691,7 @@ julia> sqrt(A)
  0.0  2.0
 ```
 """
-function sqrt(A::StridedMatrix{<:Real})
+function sqrt(A::AbstractMatrix{<:Real})
     if issymmetric(A)
         return copytri!(parent(sqrt(Symmetric(A))), 'U')
     end
@@ -697,7 +704,7 @@ function sqrt(A::StridedMatrix{<:Real})
         return SchurF.vectors * R * SchurF.vectors'
     end
 end
-function sqrt(A::StridedMatrix{<:Complex})
+function sqrt(A::AbstractMatrix{<:Complex})
     if ishermitian(A)
         sqrtHermA = sqrt(Hermitian(A))
         return isa(sqrtHermA, Hermitian) ? copytri!(parent(sqrtHermA), 'U', true) : parent(sqrtHermA)
@@ -1147,7 +1154,7 @@ julia> factorize(A) # factorize will check to see that A is already factorized
 This returns a `5×5 Bidiagonal{Float64}`, which can now be passed to other linear algebra functions
 (e.g. eigensolvers) which will use specialized methods for `Bidiagonal` types.
 """
-function factorize(A::StridedMatrix{T}) where T
+function factorize(A::AbstractMatrix{T}) where T
     m, n = size(A)
     if m == n
         if m == 1 return A[1] end
@@ -1269,7 +1276,7 @@ julia> M * N
 
 [^KY88]: Konstantinos Konstantinides and Kung Yao, "Statistical analysis of effective singular values in matrix rank determination", IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. [doi:10.1109/29.1585](http://dx.doi.org/10.1109/29.1585)
 """
-function pinv(A::StridedMatrix{T}, tol::Real) where T
+function pinv(A::AbstractMatrix{T}, tol::Real) where T
     m, n = size(A)
     Tout = typeof(zero(T)/sqrt(one(T) + one(T)))
     if m == 0 || n == 0
@@ -1298,7 +1305,7 @@ function pinv(A::StridedMatrix{T}, tol::Real) where T
     Sinv[findall(.!isfinite.(Sinv))] = zero(Stype)
     return SVD.Vt' * (Diagonal(Sinv) * SVD.U')
 end
-function pinv(A::StridedMatrix{T}) where T
+function pinv(A::AbstractMatrix{T}) where T
     tol = eps(real(float(one(T))))*min(size(A)...)
     return pinv(A, tol)
 end
@@ -1329,7 +1336,7 @@ julia> nullspace(M)
  1.0
 ```
 """
-function nullspace(A::StridedMatrix{T}) where T
+function nullspace(A::AbstractMatrix{T}) where T
     m, n = size(A)
     (m == 0 || n == 0) && return Matrix{T}(I, n, n)
     SVD = svdfact(A, full = true)
@@ -1396,7 +1403,7 @@ julia> A*X + X*B + C
  -3.77476e-15  4.44089e-16
 ```
 """
-function sylvester(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T}) where T<:BlasFloat
+function sylvester(A::StridedMatrix{T}, B::StridedMatrix{T}, C::StridedMatrix{T}) where T<:BlasFloat
     RA, QA = schur(A)
     RB, QB = schur(B)