diff --git a/NEWS.md b/NEWS.md index 16a93a529b1a9..914f70037624c 100644 --- a/NEWS.md +++ b/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 ----------------------------- diff --git a/base/abstractarray.jl b/base/abstractarray.jl index 0a4e706bd6ba0..5d595e1f1cf0c 100644 --- a/base/abstractarray.jl +++ b/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 diff --git a/base/reinterpretarray.jl b/base/reinterpretarray.jl index 126cbac9eee6b..e28d76b740662 100644 --- a/base/reinterpretarray.jl +++ b/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)) diff --git a/base/reshapedarray.jl b/base/reshapedarray.jl index bf1603e1693a9..30213fe19c70d 100644 --- a/base/reshapedarray.jl +++ b/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)) diff --git a/base/subarray.jl b/base/subarray.jl index 206abb46e1561..a469501fc1d5e 100644 --- a/base/subarray.jl +++ b/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) diff --git a/doc/src/manual/interfaces.md b/doc/src/manual/interfaces.md index 9441ae1e9303b..efbc700ede1a8 100644 --- a/doc/src/manual/interfaces.md +++ b/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. | diff --git a/stdlib/IterativeEigensolvers/test/runtests.jl b/stdlib/IterativeEigensolvers/test/runtests.jl index 6cbad717dec5d..9bad970c0ea93 100644 --- a/stdlib/IterativeEigensolvers/test/runtests.jl +++ b/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 diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index eeef2b912be53..283cfcda76a88 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/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") diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index 04c5c12448a94..896ca29747729 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -1,7 +1,7 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license using Base: @propagate_inbounds, _return_type, _default_type, @_inline_meta -import Base: length, size, axes, IndexStyle, getindex, setindex!, parent, vec, convert, similar +import Base: length, size, axes, IndexStyle, getindex, setindex!, parent, vec, convert, similar, conj ### basic definitions (types, aliases, constructors, abstractarray interface, sundry similar) @@ -118,6 +118,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)) @@ -126,6 +128,37 @@ axes(v::AdjOrTransAbsVec) = (Base.OneTo(1), axes(v.parent)...) axes(A::AdjOrTransAbsMat) = reverse(axes(A.parent)) IndexStyle(::Type{<:AdjOrTransAbsVec}) = IndexLinear() IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian() + + +# MemoryLayout of transposed and adjoint matrices +struct ConjLayout{ML<:MemoryLayout} <: MemoryLayout + layout::ML +end + +conjlayout(_1, _2) = UnknownLayout() +conjlayout(::Type{<:Complex}, M::ConjLayout) = M.layout +conjlayout(::Type{<:Complex}, M::AbstractStridedLayout) = ConjLayout(M) +conjlayout(::Type{<:Real}, M::MemoryLayout) = M + + +Base.subarraylayout(M::ConjLayout, t::Tuple) = ConjLayout(Base.subarraylayout(M.layout, t)) + +MemoryLayout(A::Transpose) = transposelayout(MemoryLayout(parent(A))) +MemoryLayout(A::Adjoint) = adjointlayout(eltype(A), MemoryLayout(parent(A))) +transposelayout(_) = UnknownLayout() +transposelayout(::StridedLayout) = StridedLayout() +transposelayout(::ColumnMajor) = RowMajor() +transposelayout(::RowMajor) = ColumnMajor() +transposelayout(::DenseColumnMajor) = DenseRowMajor() +transposelayout(::DenseRowMajor) = DenseColumnMajor() +transposelayout(M::ConjLayout) = ConjLayout(transposelayout(M.layout)) +adjointlayout(::Type{T}, M::MemoryLayout) where T = transposelayout(conjlayout(T, M)) + + +# Adjoints and transposes conform to the strided array interface if their parent does +Base.unsafe_convert(::Type{Ptr{T}}, A::AdjOrTrans{T,S}) where {T,S} = Base.unsafe_convert(Ptr{T}, parent(A)) +strides(A::AdjOrTrans) = (stride(parent(A),2), stride(parent(A),1)) + @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) @@ -137,6 +170,15 @@ IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian() # conversion of underlying storage convert(::Type{Adjoint{T,S}}, A::Adjoint) where {T,S} = Adjoint{T,S}(convert(S, A.parent)) convert(::Type{Transpose{T,S}}, A::Transpose) where {T,S} = Transpose{T,S}(convert(S, A.parent)) +convert(::Type{AbstractArray{T}}, A::Transpose{T}) where {T} = A +convert(::Type{AbstractArray{T}}, A::Adjoint{T}) where {T} = A +convert(::Type{AbstractMatrix{T}}, A::Transpose{T}) where {T} = A +convert(::Type{AbstractMatrix{T}}, A::Adjoint{T}) where {T} = A +convert(::Type{AbstractArray{T}}, A::Adjoint) where {T} = Adjoint(convert(AbstractArray{T}, A.parent)) +convert(::Type{AbstractArray{T}}, A::Transpose) where {T} = Transpose(convert(AbstractArray{T}, A.parent)) +convert(::Type{AbstractMatrix{T}}, A::Adjoint) where {T} = Adjoint(convert(AbstractArray{T}, A.parent)) +convert(::Type{AbstractMatrix{T}}, A::Transpose) where {T} = Transpose(convert(AbstractArray{T}, A.parent)) + # for vectors, the semantics of the wrapped and unwrapped types differ # so attempt to maintain both the parent and wrapper type insofar as possible @@ -193,10 +235,8 @@ broadcast(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs... # Adjoint/Transpose-vector * vector *(u::AdjointAbsVec, v::AbstractVector) = dot(u.parent, v) *(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = dot(u.parent, v) -function *(u::TransposeAbsVec, v::AbstractVector) - @boundscheck length(u) == length(v) || throw(DimensionMismatch()) - return sum(@inbounds(u[k]*v[k]) for k in 1:length(u)) -end +*(u::TransposeAbsVec, v::AbstractVector) = dotu(u.parent, v) + # vector * Adjoint/Transpose-vector *(u::AbstractVector, v::AdjOrTransAbsVec) = broadcast(*, u, v) # Adjoint/Transpose-vector * Adjoint/Transpose-vector @@ -226,16 +266,3 @@ pinv(v::TransposeAbsVec, tol::Real = 0) = pinv(conj(v.parent)).parent /(u::TransposeAbsVec, A::AbstractMatrix) = transpose(transpose(A) \ u.parent) /(u::AdjointAbsVec, A::Transpose{<:Any,<:AbstractMatrix}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent) /(u::TransposeAbsVec, A::Adjoint{<:Any,<:AbstractMatrix}) = transpose(conj(A.parent) \ u.parent) # technically should be transpose(copy(transpose(copy(A))) \ u.parent) - -# dismabiguation methods -*(A::AdjointAbsVec, B::Transpose{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::TransposeAbsVec, B::Adjoint{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractMatrix}) = copy(A) * B -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:AbstractMatrix}) = A * copy(B) -# Adj/Trans-vector * Trans/Adj-vector, shouldn't exist, here for ambiguity resolution? TODO: test removal -*(A::Adjoint{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:AbstractVector}) = throw(MethodError(*, (A, B))) -*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractVector}) = throw(MethodError(*, (A, B))) -# Adj/Trans-matrix * Trans/Adj-vector, shouldn't exist, here for ambiguity resolution? TODO: test removal -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractVector}) = throw(MethodError(*, (A, B))) -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:AbstractVector}) = throw(MethodError(*, (A, B))) -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractVector}) = throw(MethodError(*, (A, B))) diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl index 9d2d617b78023..efdff7462019c 100644 --- a/stdlib/LinearAlgebra/src/bidiag.jl +++ b/stdlib/LinearAlgebra/src/bidiag.jl @@ -341,8 +341,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))) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 101b6c7230e63..aa69d24f596d0 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -35,7 +35,7 @@ function rmul!(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 @@ -126,7 +126,11 @@ function stride(a::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray end return s end -strides(a::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}) = size_to_strides(1, size(a)...) +strides(a::DenseArray) = size_to_strides(1, size(a)...) +strides(a::ReshapedArray) = _dense_strides(size(a), MemoryLayout(parent(a))) +strides(a::ReinterpretArray) = _dense_strides(size(a), MemoryLayout(parent(a))) +_dense_strides(sz, ::DenseColumnMajor) = size_to_strides(1, sz...) +_dense_strides(sz, ::DenseRowMajor) = reverse(size_to_strides(1, sz...)) function norm(x::StridedVector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasFloat,TI<:Integer} if minimum(rx) < 1 || maximum(rx) > length(x) @@ -496,12 +500,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}, ::AbstractColumnMajor) where T<:BlasFloat n = checksquare(A) if ishermitian(A) return copytri!(parent(exp(Hermitian(A))), 'U', true) @@ -585,7 +596,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 +606,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 +636,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 +695,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 +708,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 +1158,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 @@ -1271,7 +1282,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](https://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 @@ -1300,7 +1311,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 @@ -1408,7 +1419,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) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 2ecd4220095b4..29e20d62c802b 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -599,7 +599,10 @@ norm(v::AdjointAbsVec, q::Real) = q == Inf ? norm(conj(v.parent), 1) : norm(conj norm(v::AdjointAbsVec) = norm(conj(v.parent)) norm(v::TransposeAbsVec) = norm(v.parent) -function vecdot(x::AbstractArray, y::AbstractArray) +vecdot(x::AbstractArray, y::AbstractArray) = vecdot(vec(x), vec(y)) +vecdot(x::AbstractVector, y::AbstractVector) = _vecdot(x, y, MemoryLayout(x), MemoryLayout(y)) + +function _vecdot(x::AbstractVector, y::AbstractVector, _1, _2) lx = _length(x) if lx != _length(y) throw(DimensionMismatch("first array has length $(lx) which does not match the length of the second, $(_length(y)).")) @@ -707,6 +710,12 @@ end # Call optimized BLAS methods for vectors of numbers dot(x::AbstractVector{<:Number}, y::AbstractVector{<:Number}) = vecdot(x, y) +dotu(x::AbstractVector, y::AbstractVector) = _dotu(x, y, MemoryLayout(x), MemoryLayout(y)) +_dotu(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, _1, _2) = dot(x, y) +function _dotu(x::AbstractVector, y::AbstractVector, _1, _2) + @boundscheck length(x) == length(y) || throw(DimensionMismatch()) + return sum(@inbounds(return x[k]*y[k]) for k in 1:length(x)) +end ########################################################################################### diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 1c5bcbf027322..adefea58b1e64 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -6,8 +6,16 @@ matprod(x, y) = x*y + x*y # Dot products -vecdot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y) -vecdot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotc(x, y) +for (typ, bdot) in ((:BlasReal, :(BLAS.dot)), (:BlasComplex, :(BLAS.dotc))) + @eval begin + _vecdot(x::AbstractVector{T}, y::AbstractVector{T}, + ::AbstractStridedLayout, ::AbstractStridedLayout) where T<:$typ = $bdot(x, y) + end + end + +_dotu(x::AbstractVector{T}, y::AbstractVector{T}, + ::AbstractStridedLayout, ::AbstractStridedLayout) where {T<:BlasComplex} = BLAS.dotu(x, y) + function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasReal,TI<:Integer} if length(rx) != length(ry) @@ -35,30 +43,45 @@ function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector GC.@preserve x y BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry)) end -*(transx::Transpose{<:Any,<:StridedVector{T}}, y::StridedVector{T}) where {T<:BlasComplex} = - (x = transx.parent; BLAS.dotu(x, y)) - # Matrix-vector multiplication -function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S} +(*)(A::AbstractVecOrMat, B::AbstractVecOrMat) = _mul(A, B, MemoryLayout(A), MemoryLayout(B)) + +function _mul(A::AbstractMatrix{T}, x::AbstractVector{S}, ::StridedLayout, ::StridedLayout) where {T<:BlasFloat,S} TS = promote_op(matprod, T, S) mul!(similar(x, TS, size(A,1)), A, convert(AbstractVector{TS}, x)) end -function (*)(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S} +function _mul(A::AbstractMatrix{T}, x::AbstractVector{S}, _1, _2) where {T,S} TS = promote_op(matprod, T, S) - mul!(similar(x,TS,size(A,1)),A,x) + mul!(similar(x,TS,size(A,1)), A, x) end -# these will throw a DimensionMismatch unless B has 1 row (or 1 col for transposed case): -*(a::AbstractVector, transB::Transpose{<:Any,<:AbstractMatrix}) = - (B = transB.parent; *(reshape(a,length(a),1), transpose(B))) -*(a::AbstractVector, adjB::Adjoint{<:Any,<:AbstractMatrix}) = - (B = adjB.parent; *(reshape(a,length(a),1), adjoint(B))) -(*)(a::AbstractVector, B::AbstractMatrix) = reshape(a,length(a),1)*B +""" + mul!(Y, A, B) -> Y + +Calculates the matrix-matrix or matrix-vector product ``AB`` and stores the result in `Y`, +overwriting the existing value of `Y`. Note that `Y` must not be aliased with either `A` or +`B`. + +# Examples +```jldoctest +julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); mul!(Y, A, B); + +julia> Y +2×2 Array{Float64,2}: + 3.0 3.0 + 7.0 7.0 +``` +""" +mul!(C::AbstractVecOrMat, A::AbstractVecOrMat, B::AbstractVecOrMat) = _mul!(C, A, B, MemoryLayout(C), MemoryLayout(A), MemoryLayout(B)) + +_mul!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector, _1, _2, _3) = generic_matvecmul!(y, 'N', A, x) +_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector, ::AbstractStridedLayout, ::AbstractColumnMajor, _) where {T<:BlasFloat} = + mul!(y, A, convert(Vector{T}, x)) +_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, ::AbstractStridedLayout, ::AbstractColumnMajor, ::AbstractStridedLayout) where {T<:BlasFloat} = gemv!(y, 'N', A, x) -mul!(y::StridedVector{T}, A::StridedVecOrMat{T}, x::StridedVector{T}) where {T<:BlasFloat} = gemv!(y, 'N', A, x) for elty in (Float32,Float64) @eval begin - function mul!(y::StridedVector{Complex{$elty}}, A::StridedVecOrMat{Complex{$elty}}, x::StridedVector{$elty}) + function _mul!(y::AbstractVector{Complex{$elty}}, A::AbstractMatrix{Complex{$elty}}, x::AbstractVector{$elty}, ::AbstractStridedLayout, ::AbstractColumnMajor, ::AbstractStridedLayout) Afl = reinterpret($elty,A) yfl = reinterpret($elty,y) gemv!(yfl,'N',Afl,x) @@ -66,43 +89,26 @@ for elty in (Float32,Float64) end end end -mul!(y::AbstractVector, A::AbstractVecOrMat, x::AbstractVector) = generic_matvecmul!(y, 'N', A, x) -function *(transA::Transpose{<:Any,<:StridedMatrix{T}}, x::StridedVector{S}) where {T<:BlasFloat,S} - A = transA.parent - TS = promote_op(matprod, T, S) - mul!(similar(x,TS,size(A,2)), transpose(A), convert(AbstractVector{TS}, x)) -end -function *(transA::Transpose{<:Any,<:AbstractMatrix{T}}, x::AbstractVector{S}) where {T,S} - A = transA.parent - TS = promote_op(matprod, T, S) - mul!(similar(x,TS,size(A,2)), transpose(A), x) -end -mul!(y::StridedVector{T}, transA::Transpose{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T}) where {T<:BlasFloat} = - (A = transA.parent; gemv!(y, 'T', A, x)) -mul!(y::AbstractVector, transA::Transpose{<:Any,<:AbstractVecOrMat}, x::AbstractVector) = - (A = transA.parent; generic_matvecmul!(y, 'T', A, x)) +_mul!(y::AbstractVector, transA::AbstractMatrix, x::AbstractVector, ::AbstractStridedLayout, ::AbstractRowMajor, ::AbstractStridedLayout) = + (A = transpose(transA); generic_matvecmul!(y, 'T', A, x)) +_mul!(y::AbstractVector, adjA::AbstractMatrix, x::AbstractVector, ::AbstractStridedLayout, ::ConjLayout{<:AbstractRowMajor}, ::AbstractStridedLayout) = + (A = adjoint(adjA); generic_matvecmul!(y, 'C', A, x)) +_mul!(y::AbstractVector{T}, adjA::AbstractMatrix{T}, x::AbstractVector{T}, ::AbstractStridedLayout, ::AbstractRowMajor, ::AbstractStridedLayout) where {T<:BlasFloat} = + gemv!(y, 'T', transpose(adjA), x) +_mul!(y::AbstractVector{T}, adjA::AbstractMatrix{T}, x::AbstractVector{T}, ::AbstractStridedLayout, ::ConjLayout{<:AbstractRowMajor}, ::AbstractStridedLayout) where {T<:BlasComplex} = + gemv!(y, 'C', adjoint(adjA), x) -function *(adjA::Adjoint{<:Any,<:StridedMatrix{T}}, x::StridedVector{S}) where {T<:BlasFloat,S} - A = adjA.parent - TS = promote_op(matprod, T, S) - mul!(similar(x,TS,size(A,2)), adjoint(A) ,convert(AbstractVector{TS},x)) -end -function *(adjA::Adjoint{<:Any,<:AbstractMatrix{T}}, x::AbstractVector{S}) where {T,S} - A = adjA.parent - TS = promote_op(matprod, T, S) - mul!(similar(x,TS,size(A,2)), adjoint(A), x) -end +# Vector-matrix multiplication -mul!(y::StridedVector{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T}) where {T<:BlasReal} = - (A = adjA.parent; mul!(y, transpose(A), x)) -mul!(y::StridedVector{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T}) where {T<:BlasComplex} = - (A = adjA.parent; gemv!(y, 'C', A, x)) -mul!(y::AbstractVector, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, x::AbstractVector) = - (A = adjA.parent; generic_matvecmul!(y, 'C', A, x)) +_mul(a::AbstractVector, B::AbstractMatrix, _1, _2) = reshape(a,length(a),1)*B -# Matrix-matrix multiplication +# these enable treating vectors as n x 1 matrices for important use-cases +mul!(y::AbstractMatrix, a::AbstractVector, B::AbstractMatrix) = mul!(y, reshape(a,length(a),1), B) +mul!(y::AbstractVector, a::AbstractVector, B::AbstractMatrix) = mul!(reshape(y,length(y),1), reshape(a,length(a),1), B) +mul!(y::AbstractVector, A::AbstractMatrix, B::AbstractMatrix) = mul!(reshape(y,length(b),1), A, B) +# Matrix-matrix multiplication """ *(A::AbstractMatrix, B::AbstractMatrix) @@ -116,76 +122,56 @@ julia> [1 1; 0 1] * [1 0; 1 1] 1 1 ``` """ -function (*)(A::AbstractMatrix, B::AbstractMatrix) +*(::AbstractMatrix, ::AbstractMatrix) + +function _mul(A::AbstractMatrix, B::AbstractMatrix, _1, _2) TS = promote_op(matprod, eltype(A), eltype(B)) mul!(similar(B, TS, (size(A,1), size(B,2))), A, B) end -mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = gemm_wrapper!(C, 'N', 'N', A, B) -for elty in (Float32,Float64) - @eval begin - function mul!(C::StridedMatrix{Complex{$elty}}, A::StridedVecOrMat{Complex{$elty}}, B::StridedVecOrMat{$elty}) - Afl = reinterpret($elty, A) - Cfl = reinterpret($elty, C) - gemm_wrapper!(Cfl, 'N', 'N', Afl, B) - return C - end - end -end - -""" - mul!(Y, A, B) -> Y -Calculates the matrix-matrix or matrix-vector product ``AB`` and stores the result in `Y`, -overwriting the existing value of `Y`. Note that `Y` must not be aliased with either `A` or -`B`. - -# Examples -```jldoctest -julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); mul!(Y, A, B); - -julia> Y -2×2 Array{Float64,2}: - 3.0 3.0 - 7.0 7.0 -``` -""" -mul!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat) = generic_matmatmul!(C, 'N', 'N', A, B) """ rmul!(A, B) Calculate the matrix-matrix product ``AB``, overwriting `A`, and return the result. """ -rmul!(A, B) +rmul!(A, B) = _rmul!(A, B, MemoryLayout(A), MemoryLayout(B)) +_rmul!(A, B, _1, _2) = copyto!(A, mul!(similar(A), A, B)) """ lmul!(A, B) Calculate the matrix-matrix product ``AB``, overwriting `B`, and return the result. """ -lmul!(A, B) +lmul!(A, B) = _lmul!(A, B, MemoryLayout(A), MemoryLayout(B)) +_lmul!(A, B, _1, _2) = copyto!(B, mul!(similar(B), A, B)) -function *(transA::Transpose{<:Any,<:AbstractMatrix}, B::AbstractMatrix) - A = transA.parent - TS = promote_op(matprod, eltype(A), eltype(B)) - mul!(similar(B, TS, (size(A,2), size(B,2))), transpose(A), B) -end -mul!(C::StridedMatrix{T}, transA::Transpose{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = - (A = transA.parent; A===B ? syrk_wrapper!(C, 'T', A) : gemm_wrapper!(C, 'T', 'N', A, B)) -mul!(C::AbstractMatrix, transA::Transpose{<:Any,<:AbstractVecOrMat}, B::AbstractVecOrMat) = - (A = transA.parent; generic_matmatmul!(C, 'T', 'N', A, B)) +_mul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, _1, _2, _3) = generic_matmatmul!(C, 'N', 'N', A, B) -function *(A::AbstractMatrix, transB::Transpose{<:Any,<:AbstractMatrix}) - B = transB.parent - TS = promote_op(matprod, eltype(A), eltype(B)) - mul!(similar(B, TS, (size(A,1), size(B,1))), A, transpose(B)) +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, ::AbstractStridedLayout, ::AbstractStridedLayout, ::AbstractStridedLayout) where {T<:BlasFloat} = + gemm_wrapper!(C, 'N', 'N', A, B) +for elty in (Float32,Float64) + @eval begin + function _mul!(C::AbstractMatrix{Complex{$elty}}, A::AbstractMatrix{Complex{$elty}}, B::AbstractMatrix{$elty}, ::AbstractColumnMajor, ::AbstractColumnMajor, ::AbstractColumnMajor) + Afl = reinterpret($elty, A) + Cfl = reinterpret($elty, C) + gemm_wrapper!(Cfl, 'N', 'N', Afl, B) + return C + end + end end -mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, transB::Transpose{<:Any,<:StridedVecOrMat{T}}) where {T<:BlasFloat} = - (B = transB.parent; A===B ? syrk_wrapper!(C, 'N', A) : gemm_wrapper!(C, 'N', 'T', A, B)) + +_mul!(C::AbstractMatrix{T}, transA::AbstractMatrix{T}, B::AbstractMatrix{T}, ::AbstractColumnMajor, ::AbstractRowMajor, ::AbstractColumnMajor) where {T<:BlasFloat} = + (A = transpose(transA); A === B ? syrk_wrapper!(C, 'T', A) : gemm_wrapper!(C, 'T', 'N', A, B)) +_mul!(C::AbstractMatrix, transA::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::AbstractRowMajor, ::AbstractStridedLayout) = + (A = transpose(transA); generic_matmatmul!(C, 'T', 'N', A, B)) +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, transB::AbstractMatrix{T}, ::AbstractColumnMajor, ::AbstractColumnMajor, ::AbstractRowMajor) where {T<:BlasFloat} = + (B = transpose(transB); A === B ? syrk_wrapper!(C, 'N', A) : gemm_wrapper!(C, 'N', 'T', A, B)) + for elty in (Float32,Float64) @eval begin - function mul!(C::StridedMatrix{Complex{$elty}}, A::StridedVecOrMat{Complex{$elty}}, transB::Transpose{<:Any,<:StridedVecOrMat{$elty}}) - B = transB.parent + function _mul!(C::AbstractMatrix{Complex{$elty}}, A::AbstractMatrix{Complex{$elty}}, transB::AbstractMatrix{$elty}, ::AbstractColumnMajor, ::AbstractColumnMajor, ::AbstractRowMajor) + B = transpose(transB) Afl = reinterpret($elty, A) Cfl = reinterpret($elty, C) gemm_wrapper!(Cfl, 'N', 'T', Afl, B) @@ -193,66 +179,29 @@ for elty in (Float32,Float64) end end end -# collapsing the following two defs with C::AbstractVecOrMat yields ambiguities -mul!(C::AbstractVector, A::AbstractVecOrMat, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - _disambigmul!(C, A, transB) -mul!(C::AbstractMatrix, A::AbstractVecOrMat, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - _disambigmul!(C, A, transB) -_disambigmul!(C::AbstractVecOrMat, A::AbstractVecOrMat, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (B = transB.parent; generic_matmatmul!(C, 'N', 'T', A, B)) - -# collapsing the following two defs with transB::Transpose{<:Any,<:AbstractVecOrMat{S}} yields ambiguities -*(transA::Transpose{<:Any,<:AbstractMatrix}, transB::Transpose{<:Any,<:AbstractMatrix}) = - _disambigmul(transA, transB) -*(transA::Transpose{<:Any,<:AbstractMatrix}, transB::Transpose{<:Any,<:AbstractVector}) = - _disambigmul(transA, transB) -function _disambigmul(transA::Transpose{<:Any,<:AbstractMatrix{T}}, transB::Transpose{<:Any,<:AbstractVecOrMat{S}}) where {T,S} - A, B = transA.parent, transB.parent - TS = promote_op(matprod, T, S) - mul!(similar(B, TS, (size(A,2), size(B,1))), transpose(A), transpose(B)) -end -mul!(C::StridedMatrix{T}, transA::Transpose{<:Any,<:StridedVecOrMat{T}}, transB::Transpose{<:Any,<:StridedVecOrMat{T}}) where {T<:BlasFloat} = - (A = transA.parent; B = transB.parent; gemm_wrapper!(C, 'T', 'T', A, B)) -mul!(C::AbstractMatrix, transA::Transpose{<:Any,<:AbstractVecOrMat}, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (A = transA.parent; B = transB.parent; generic_matmatmul!(C, 'T', 'T', A, B)) -mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractVecOrMat}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) - -*(adjA::Adjoint{<:Any,<:StridedMatrix{T}}, B::StridedMatrix{T}) where {T<:BlasReal} = - (A = adjA.parent; *(transpose(A), B)) -mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = - (A = adjA.parent; mul!(C, transpose(A), B)) -function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, B::AbstractMatrix) - A = adjA.parent - TS = promote_op(matprod, eltype(A), eltype(B)) - mul!(similar(B, TS, (size(A,2), size(B,2))), adjoint(A), B) -end -mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = - (A = adjA.parent; A===B ? herk_wrapper!(C,'C',A) : gemm_wrapper!(C,'C', 'N', A, B)) -mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, B::AbstractVecOrMat) = - (A = adjA.parent; generic_matmatmul!(C, 'C', 'N', A, B)) - -*(A::StridedMatrix{<:BlasFloat}, adjB::Adjoint{<:Any,<:StridedMatrix{<:BlasReal}}) = - (B = adjB.parent; *(A, transpose(B))) -mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:StridedVecOrMat{<:BlasReal}}) where {T<:BlasFloat} = - (B = adjB.parent; mul!(C, A, transpose(B))) -function *(A::AbstractMatrix, adjB::Adjoint{<:Any,<:AbstractMatrix}) - B = adjB.parent - TS = promote_op(matprod, eltype(A), eltype(B)) - mul!(similar(B,TS,(size(A,1),size(B,1))), A, adjoint(B)) -end -mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:StridedVecOrMat{T}}) where {T<:BlasComplex} = - (B = adjB.parent; A===B ? herk_wrapper!(C, 'N', A) : gemm_wrapper!(C, 'N', 'C', A, B)) -mul!(C::AbstractMatrix, A::AbstractVecOrMat, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) = - (B = adjB.parent; generic_matmatmul!(C, 'N', 'C', A, B)) - -*(adjA::Adjoint{<:Any,<:AbstractMatrix}, adjB::Adjoint{<:Any,<:AbstractMatrix}) = - (A = adjA.parent; B = adjB.parent; mul!(similar(B, promote_op(matprod, eltype(A), eltype(B)), (size(A,2), size(B,1))), adjoint(A), adjoint(B))) -mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, adjB::Adjoint{<:Any,<:StridedVecOrMat{T}}) where {T<:BlasFloat} = - (A = adjA.parent; B = adjB.parent; gemm_wrapper!(C, 'C', 'C', A, B)) -mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) = - (A = adjA.parent; B = adjB.parent; generic_matmatmul!(C, 'C', 'C', A, B)) -mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (A = adjA.parent; B = transB.parent; generic_matmatmul!(C, 'C', 'T', A, B)) + +_mul!(C::AbstractMatrix{T}, transA::AbstractMatrix{T}, transB::AbstractMatrix{T}, ::AbstractColumnMajor, ::AbstractRowMajor, ::AbstractRowMajor) where {T<:BlasFloat} = + (A = transpose(transA); B = transpose(transB); gemm_wrapper!(C, 'T', 'T', A, B)) +_mul!(C::AbstractMatrix, transA::AbstractMatrix, transB::AbstractMatrix, ::AbstractStridedLayout, ::AbstractRowMajor, ::AbstractRowMajor) = + (A = transpose(transA); B = transpose(transB); generic_matmatmul!(C, 'T', 'T', A, B)) + +_mul!(C::AbstractMatrix{T}, adjA::AbstractMatrix{T}, B::AbstractMatrix{T}, ::AbstractColumnMajor, ::ConjLayout{<:AbstractRowMajor}, ::AbstractColumnMajor) where {T<:BlasComplex} = + (A = adjoint(adjA); A===B ? herk_wrapper!(C,'C',A) : gemm_wrapper!(C,'C', 'N', A, B)) +_mul!(C::AbstractMatrix, adjA::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::ConjLayout{<:AbstractRowMajor}, ::AbstractStridedLayout) = + (A = adjoint(adjA); generic_matmatmul!(C, 'C', 'N', A, B)) + +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, adjB::AbstractMatrix{T}, ::AbstractColumnMajor, ::AbstractColumnMajor, ::ConjLayout{<:AbstractRowMajor}) where {T<:BlasComplex} = + (B = adjoint(adjB); A===B ? herk_wrapper!(C, 'N', A) : gemm_wrapper!(C, 'N', 'C', A, B)) +_mul!(C::AbstractMatrix, A::AbstractMatrix, adjB::AbstractMatrix, ::AbstractStridedLayout, ::AbstractStridedLayout, ::ConjLayout{<:AbstractRowMajor}) = + (B = adjoint(adjB); generic_matmatmul!(C, 'N', 'C', A, B)) + +_mul!(C::AbstractMatrix{T}, adjA::AbstractMatrix{T}, adjB::AbstractMatrix{T}, ::AbstractColumnMajor, ::ConjLayout{<:AbstractRowMajor}, ::ConjLayout{<:AbstractRowMajor}) where {T<:BlasFloat} = + (A = adjoint(adjA); B = adjoint(adjB); gemm_wrapper!(C, 'C', 'C', A, B)) +_mul!(C::AbstractMatrix, adjA::AbstractMatrix, adjB::AbstractMatrix, ::AbstractStridedLayout, ::ConjLayout{<:AbstractRowMajor}, ::ConjLayout{<:AbstractRowMajor}) = + (A = adjoint(adjA); B = adjoint(adjB); generic_matmatmul!(C, 'C', 'C', A, B)) +_mul!(C::AbstractMatrix, adjA::AbstractMatrix, transB::AbstractMatrix, ::AbstractStridedLayout, ::ConjLayout{<:AbstractRowMajor}, ::AbstractRowMajor) = + (A = adjoint(adjA); B = transpose(transB); generic_matmatmul!(C, 'C', 'T', A, B)) + # Supporting functions for matrix multiplication function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false) @@ -271,7 +220,7 @@ function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false) A end -function gemv!(y::StridedVector{T}, tA::AbstractChar, A::StridedVecOrMat{T}, x::StridedVector{T}) where T<:BlasFloat +function gemv!(y::AbstractVector{T}, tA::Char, A::AbstractVecOrMat{T}, x::AbstractVector{T}) where T<:BlasFloat mA, nA = lapack_size(tA, A) if nA != length(x) throw(DimensionMismatch("second dimension of A, $nA, does not match length of x, $(length(x))")) @@ -348,18 +297,19 @@ function herk_wrapper!(C::Union{StridedMatrix{T}, StridedMatrix{Complex{T}}}, tA return generic_matmatmul!(C,tA, tAt, A, A) end -function gemm_wrapper(tA::AbstractChar, tB::AbstractChar, - A::StridedVecOrMat{T}, - B::StridedVecOrMat{T}) where T<:BlasFloat +function gemm_wrapper(tA::Char, tB::Char, + A::AbstractVecOrMat{T}, + B::AbstractVecOrMat{T}) where T<:BlasFloat mA, nA = lapack_size(tA, A) mB, nB = lapack_size(tB, B) C = similar(B, T, mA, nB) gemm_wrapper!(C, tA, tB, A, B) end -function gemm_wrapper!(C::StridedVecOrMat{T}, tA::AbstractChar, tB::AbstractChar, - A::StridedVecOrMat{T}, - B::StridedVecOrMat{T}) where T<:BlasFloat +#gemm_wrapper! requires A, B, and C to conform to the strided array interface +function gemm_wrapper!(C::AbstractVecOrMat{T}, tA::Char, tB::Char, + A::AbstractVecOrMat{T}, + B::AbstractVecOrMat{T}) where T<:BlasFloat mA, nA = lapack_size(tA, A) mB, nB = lapack_size(tB, B) @@ -479,7 +429,7 @@ function generic_matvecmul!(C::AbstractVector{R}, tA, A::AbstractVecOrMat, B::Ab C end -function generic_matmatmul(tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S}) where {T,S} +function (tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S}) where {T,S} mA, nA = lapack_size(tA, A) mB, nB = lapack_size(tB, B) C = similar(B, promote_op(matprod, T, S), mA, nB) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 4c2f5e44ff507..77188527c7c4a 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -240,6 +240,9 @@ end Array(A::Union{Symmetric,Hermitian}) = convert(Matrix, A) parent(A::HermOrSym) = A.data +symmetricdata(A::Symmetric) = A.data +symmetricdata(A::Hermitian{<:Real}) = A.data +hermitiandata(A::Hermitian) = A.data Symmetric{T,S}(A::Symmetric{T,S}) where {T,S<:AbstractMatrix} = A Symmetric{T,S}(A::Symmetric) where {T,S<:AbstractMatrix} = Symmetric{T,S}(convert(S,A.data),A.uplo) AbstractMatrix{T}(A::Symmetric) where {T} = Symmetric(convert(AbstractMatrix{T}, A.data), Symbol(A.uplo)) @@ -247,6 +250,67 @@ Hermitian{T,S}(A::Hermitian{T,S}) where {T,S<:AbstractMatrix} = A Hermitian{T,S}(A::Hermitian) where {T,S<:AbstractMatrix} = Hermitian{T,S}(convert(S,A.data),A.uplo) AbstractMatrix{T}(A::Hermitian) where {T} = Hermitian(convert(AbstractMatrix{T}, A.data), Symbol(A.uplo)) +# MemoryLayout of Symmetric/Hermitian +""" + SymmetricLayout(layout, uplo) + + +is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory +as a symmetrized version of `layout`, where the entries used are dictated by the +`uplo`, which can be `'U'` or `L'`. + +A matrix that has memory layout `SymmetricLayout(layout, uplo)` must overrided +`symmetricdata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and +`A[k,j] == B[k,j]` for `j ≥ k` if `uplo == 'U'` (`j ≤ k` if `uplo == 'L'`) and +`A[k,j] == B[j,k]` for `j < k` if `uplo == 'U'` (`j > k` if `uplo == 'L'`). +""" +struct SymmetricLayout{ML<:MemoryLayout} <: MemoryLayout + layout::ML + uplo::Char +end +SymmetricLayout(layout::ML, uplo) where ML<:MemoryLayout = SymmetricLayout{ML}(layout, uplo) + +""" + HermitianLayout(layout, uplo) + + +is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory +as a hermitianized version of `layout`, where the entries used are dictated by the +`uplo`, which can be `'U'` or `L'`. + +A matrix that has memory layout `HermitianLayout(layout, uplo)` must overrided +`hermitiandata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and +`A[k,j] == B[k,j]` for `j ≥ k` if `uplo == 'U'` (`j ≤ k` if `uplo == 'L'`) and +`A[k,j] == conj(B[j,k])` for `j < k` if `uplo == 'U'` (`j > k` if `uplo == 'L'`). +""" +struct HermitianLayout{ML<:MemoryLayout} <: MemoryLayout + layout::ML + uplo::Char +end +HermitianLayout(layout::ML, uplo) where ML<:MemoryLayout = HermitianLayout{ML}(layout, uplo) + +MemoryLayout(A::Hermitian) = hermitianlayout(eltype(A), MemoryLayout(parent(A)), A.uplo) +MemoryLayout(A::Symmetric) = symmetriclayout(MemoryLayout(parent(A)), A.uplo) +hermitianlayout(_1, _2, _3) = UnknownLayout() +hermitianlayout(::Type{<:Complex}, layout::AbstractColumnMajor, uplo) = HermitianLayout(layout,uplo) +hermitianlayout(::Type{<:Real}, layout::AbstractColumnMajor, uplo) = SymmetricLayout(layout,uplo) +hermitianlayout(::Type{<:Complex}, layout::AbstractRowMajor, uplo) = HermitianLayout(layout,uplo) +hermitianlayout(::Type{<:Real}, layout::AbstractRowMajor, uplo) = SymmetricLayout(layout,uplo) +symmetriclayout(_1, _2) = UnknownLayout() +symmetriclayout(layout::AbstractColumnMajor, uplo) = SymmetricLayout(layout,uplo) +symmetriclayout(layout::AbstractRowMajor, uplo) = SymmetricLayout(layout,uplo) +transposelayout(S::SymmetricLayout) = S +adjointlayout(::Type{T}, S::SymmetricLayout) where T<:Real = S +adjointlayout(::Type{T}, S::HermitianLayout) where T = S +Base.subarraylayout(S::SymmetricLayout, ::Tuple{<:Slice,<:Slice}) = S +Base.subarraylayout(S::HermitianLayout, ::Tuple{<:Slice,<:Slice}) = S +symmetricdata(V::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = symmetricdata(parent(V)) +symmetricdata(V::Adjoint{<:Real}) = symmetricdata(parent(V)) +symmetricdata(V::Transpose) = symmetricdata(parent(V)) +hermitiandata(V::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = hermitiandata(parent(V)) +hermitiandata(V::Adjoint) = hermitiandata(parent(V)) +hermitiandata(V::Transpose{<:Real}) = hermitiandata(parent(V)) + copy(A::Symmetric{T,S}) where {T,S} = (B = copy(A.data); Symmetric{T,typeof(B)}(B,A.uplo)) copy(A::Hermitian{T,S}) where {T,S} = (B = copy(A.data); Hermitian{T,typeof(B)}(B,A.uplo)) @@ -383,52 +447,38 @@ end (-)(A::Symmetric{Tv,S}) where {Tv,S} = Symmetric{Tv,S}(-A.data, A.uplo) (-)(A::Hermitian{Tv,S}) where {Tv,S} = Hermitian{Tv,S}(-A.data, A.uplo) -## Matvec -mul!(y::StridedVector{T}, A::Symmetric{T,<:StridedMatrix}, x::StridedVector{T}) where {T<:BlasFloat} = - BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y) -mul!(y::StridedVector{T}, A::Hermitian{T,<:StridedMatrix}, x::StridedVector{T}) where {T<:BlasReal} = - BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y) -mul!(y::StridedVector{T}, A::Hermitian{T,<:StridedMatrix}, x::StridedVector{T}) where {T<:BlasComplex} = - BLAS.hemv!(A.uplo, one(T), A.data, x, zero(T), y) -## Matmat -mul!(C::StridedMatrix{T}, A::Symmetric{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} = - BLAS.symm!('L', A.uplo, one(T), A.data, B, zero(T), C) -mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Symmetric{T,<:StridedMatrix}) where {T<:BlasFloat} = - BLAS.symm!('R', B.uplo, one(T), B.data, A, zero(T), C) -mul!(C::StridedMatrix{T}, A::Hermitian{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasReal} = - BLAS.symm!('L', A.uplo, one(T), A.data, B, zero(T), C) -mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Hermitian{T,<:StridedMatrix}) where {T<:BlasReal} = - BLAS.symm!('R', B.uplo, one(T), B.data, A, zero(T), C) -mul!(C::StridedMatrix{T}, A::Hermitian{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasComplex} = - BLAS.hemm!('L', A.uplo, one(T), A.data, B, zero(T), C) -mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Hermitian{T,<:StridedMatrix}) where {T<:BlasComplex} = - BLAS.hemm!('R', B.uplo, one(T), B.data, A, zero(T), C) - -*(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B) - -# Fallbacks to avoid generic_matvecmul!/generic_matmatmul! -## Symmetric{<:Number} and Hermitian{<:Real} are invariant to transpose; peel off the t -*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractVector) = transA.parent * B -*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractMatrix) = transA.parent * B -*(A::AbstractMatrix, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = A * transB.parent -## Hermitian{<:Number} and Symmetric{<:Real} are invariant to adjoint; peel off the c -*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::AbstractVector) = adjA.parent * B -*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::AbstractMatrix) = adjA.parent * B -*(A::AbstractMatrix, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * adjB.parent - -# ambiguities with transposed AbstractMatrix methods in linalg/matmul.jl -*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = transA.parent * transB.parent -*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, transB::Transpose{<:Any,<:RealHermSymComplexHerm}) = transA.parent * transB -*(transA::Transpose{<:Any,<:RealHermSymComplexHerm}, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = transA * transB.parent -*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = adjA.parent * adjB.parent -*(adjA::Adjoint{<:Any,<:RealHermSymComplexSym}, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = adjA * adjB.parent -*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, adjB::Adjoint{<:Any,<:RealHermSymComplexSym}) = adjA.parent * adjB - -# ambiguities with AbstractTriangular -*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractTriangular) = transA.parent * B -*(A::AbstractTriangular, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = A * transB.parent -*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::AbstractTriangular) = adjA.parent * B -*(A::AbstractTriangular, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * adjB.parent +## Matrix-vector product +_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, + ::AbstractStridedLayout, AL::SymmetricLayout{<:AbstractColumnMajor}, ::AbstractStridedLayout) where {T<:BlasFloat} = + BLAS.symv!(AL.uplo, one(T), symmetricdata(A), x, zero(T), y) +_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, + ::AbstractStridedLayout, AL::SymmetricLayout{<:AbstractRowMajor}, ::AbstractStridedLayout) where {T<:BlasFloat} = + BLAS.symv!(ifelse(AL.uplo == 'L', 'U', 'L'), one(T), transpose(symmetricdata(A)), x, zero(T), y) +_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, + ::AbstractStridedLayout, AL::HermitianLayout{<:AbstractColumnMajor}, ::AbstractStridedLayout) where {T<:BlasComplex} = + BLAS.hemv!(AL.uplo, one(T), hermitiandata(A), x, zero(T), y) +_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, + ::AbstractStridedLayout, AL::HermitianLayout{<:AbstractRowMajor}, ::AbstractStridedLayout) where {T<:BlasComplex} = + BLAS.hemv!(ifelse(AL.uplo == 'L', 'U', 'L'), one(T), transpose(hermitiandata(A)), x, zero(T), y) + +## Matrix-matrix product +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, + ::AbstractStridedLayout, AL::SymmetricLayout{<:AbstractColumnMajor}, ::AbstractStridedLayout) where {T<:BlasFloat} = + BLAS.symm!('L', AL.uplo, one(T), symmetricdata(A), B, zero(T), C) +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, + ::AbstractStridedLayout, ::AbstractStridedLayout, BL::SymmetricLayout{<:AbstractColumnMajor}) where {T<:BlasFloat} = + BLAS.symm!('R', BL.uplo, one(T), symmetricdata(B), A, zero(T), C) +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, + ::AbstractStridedLayout, AL::HermitianLayout{<:AbstractColumnMajor}, ::AbstractStridedLayout) where {T<:BlasComplex} = + BLAS.hemm!('L', AL.uplo, one(T), hermitiandata(A), B, zero(T), C) +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, + ::AbstractStridedLayout, ::AbstractStridedLayout, BL::HermitianLayout{<:AbstractColumnMajor}) where {T<:BlasComplex} = + BLAS.hemm!('R', BL.uplo, one(T), hermitiandata(B), A, zero(T), C) + +_mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, + ::AbstractStridedLayout, ::Union{HermitianLayout, SymmetricLayout}, ::Union{HermitianLayout, SymmetricLayout}) where {T<:BlasFloat} = + mul!(C, A, Matrix{T}(B)) + for T in (:Symmetric, :Hermitian), op in (:*, :/) # Deal with an ambiguous case @@ -817,35 +867,3 @@ for func in (:log, :sqrt) end end end - -# dismabiguation methods: *(Adj of RealHermSymComplexHerm, Trans of RealHermSymComplexSym) and symmetric partner -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A.parent * B.parent -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A.parent * B.parent -# dismabiguation methods: *(Adj/Trans of AbsVec/AbsMat, Adj/Trans of RealHermSymComplex{Herm|Sym}) -*(A::Adjoint{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent -*(A::Adjoint{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent -*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent -*(A::Transpose{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent -# dismabiguation methods: *(Adj/Trans of RealHermSymComplex{Herm|Sym}, Adj/Trans of AbsVec/AbsMat) -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractVector}) = A.parent * B -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractMatrix}) = A.parent * B -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractVector}) = A.parent * B -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractMatrix}) = A.parent * B -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractVector}) = A.parent * B -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractMatrix}) = A.parent * B -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractVector}) = A.parent * B -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractMatrix}) = A.parent * B - -# dismabiguation methods: *(Adj/Trans of AbsTri or RealHermSymComplex{Herm|Sym}, Adj/Trans of other) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * B.parent -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = A * B.parent -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractTriangular}) = A.parent * B -*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractTriangular}) = A.parent * B -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractTriangular}) = A.parent * B -*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractTriangular}) = A.parent * B diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index c6ef222cacc7b..32295840e48ef 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -97,6 +97,7 @@ imag(A::UnitUpperTriangular) = UpperTriangular(triu!(imag(A.data),1)) Array(A::AbstractTriangular) = Matrix(A) parent(A::AbstractTriangular) = A.data +triangulardata(A::AbstractTriangular) = A.data # then handle all methods that requires specific handling of upper/lower and unit diagonal @@ -360,6 +361,108 @@ diag(A::UnitLowerTriangular) = fill(one(eltype(A)), size(A,1)) diag(A::UpperTriangular) = diag(A.data) diag(A::UnitUpperTriangular) = fill(one(eltype(A)), size(A,1)) +# MemoryLayout of triangular matrices +abstract type AbstractTriangularLayout{ML} <: MemoryLayout end + +for memlay in (:LowerTriangularLayout, :UnitLowerTriangularLayout, + :UpperTriangularLayout, :UnitUpperTriangularLayout) + @eval begin + struct $memlay{ML<:MemoryLayout} <: AbstractTriangularLayout{ML} + layout::ML + end + end +end + +""" + LowerTriangularLayout(layout) + +is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory +equivalent to a `LowerTriangular(B)` where `B` satisfies `MemoryLayout(B) == layout`. + +A matrix that has memory layout `LowerTriangularLayout(layout)` must overrided +`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and +`A[k,j] ≡ zero(eltype(A))` for `j > k` and +`A[k,j] ≡ B[k,j]` for `j ≤ k`. + +Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory +layout `UpperTriangularLayout`. +""" +LowerTriangularLayout + +""" + UnitLowerTriangularLayout(ML::MemoryLayout) + +is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory +equivalent to a `UnitLowerTriangular(B)` where `B` satisfies `MemoryLayout(B) == layout`. + +A matrix that has memory layout `UnitLowerTriangularLayout(layout)` must overrided +`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and +`A[k,j] ≡ zero(eltype(A))` for `j > k`, +`A[k,j] ≡ one(eltype(A))` for `j == k`, +`A[k,j] ≡ B[k,j]` for `j < k`. + +Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory +layout `UnitUpperTriangularLayout`. +""" +UnitLowerTriangularLayout + +""" + UpperTriangularLayout(ML::MemoryLayout) + +is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory +equivalent to a `UpperTriangularLayout(B)` where `B` satisfies `MemoryLayout(B) == ML`. + +A matrix that has memory layout `UpperTriangularLayout(layout)` must overrided +`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and +`A[k,j] ≡ B[k,j]` for `j ≥ k` and +`A[k,j] ≡ zero(eltype(A))` for `j < k`. + +Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory +layout `LowerTriangularLayout`. +""" +UpperTriangularLayout + +""" + UnitUpperTriangularLayout(ML::MemoryLayout) + +is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory +equivalent to a `UpperTriangularLayout(B)` where `B` satisfies `MemoryLayout(B) == ML`. + +A matrix that has memory layout `UnitUpperTriangularLayout(layout)` must overrided +`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and +`A[k,j] ≡ B[k,j]` for `j > k`, +`A[k,j] ≡ one(eltype(A))` for `j == k`, +`A[k,j] ≡ zero(eltype(A))` for `j < k`. + +Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory +layout `UnitLowerTriangularLayout`. +""" +UnitUpperTriangularLayout + + +MemoryLayout(A::UpperTriangular) = triangularlayout(UpperTriangularLayout, MemoryLayout(parent(A))) +MemoryLayout(A::UnitUpperTriangular) = triangularlayout(UnitUpperTriangularLayout, MemoryLayout(parent(A))) +MemoryLayout(A::LowerTriangular) = triangularlayout(LowerTriangularLayout, MemoryLayout(parent(A))) +MemoryLayout(A::UnitLowerTriangular) = triangularlayout(UnitLowerTriangularLayout, MemoryLayout(parent(A))) +triangularlayout(_, ::MemoryLayout) = UnknownLayout() +triangularlayout(::Type{Tri}, ML::AbstractColumnMajor) where {Tri} = Tri(ML) +Base.subarraylayout(layout::AbstractTriangularLayout, ::Tuple{<:Union{Slice,Base.OneTo},<:Union{Slice,Base.OneTo}}) = layout + +for (TriLayout, TriLayoutTrans) in ((UpperTriangularLayout, LowerTriangularLayout), + (UnitUpperTriangularLayout, UnitLowerTriangularLayout), + (LowerTriangularLayout, UpperTriangularLayout), + (UnitLowerTriangularLayout, UnitUpperTriangularLayout)) + @eval begin + transposelayout(ml::$TriLayout) = $TriLayoutTrans(transposelayout(ml.layout)) + conjlayout(::Type{<:Complex}, ml::$TriLayout) = $TriLayout(ConjLayout(ml.layout)) + end +end + +triangulardata(A::Adjoint) = Adjoint(triangulardata(parent(A))) +triangulardata(A::Transpose) = Transpose(triangulardata(parent(A))) +triangulardata(A::SubArray{<:Any,2,<:Any,<:Tuple{<:Union{Slice,Base.OneTo},<:Union{Slice,Base.OneTo}}}) = + view(triangulardata(parent(A)), parentindices(A)...) + # Unary operations -(A::LowerTriangular) = LowerTriangular(-A.data) -(A::UpperTriangular) = UpperTriangular(-A.data) @@ -512,72 +615,50 @@ lmul!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B) # is this necessary? mul!(C::AbstractMatrix, A::AbstractTriangular, B::Tridiagonal) = mul!(C, copyto!(similar(parent(A)), A), B) mul!(C::AbstractMatrix, A::Tridiagonal, B::AbstractTriangular) = mul!(C, A, copyto!(similar(parent(B)), B)) -mul!(C::AbstractVector, A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (B = transB.parent; lmul!(A, transpose!(C, B))) -mul!(C::AbstractMatrix, A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractVecOrMat}) = - (B = transB.parent; lmul!(A, transpose!(C, B))) -mul!(C::AbstractMatrix, A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) = - (B = adjB.parent; lmul!(A, adjoint!(C, B))) -mul!(C::AbstractVecOrMat, A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) = - (B = adjB.parent; lmul!(A, adjoint!(C, B))) - -# The three methods for each op are neceesary to avoid ambiguities with definitions in matmul.jl -mul!(C::AbstractVector , A::AbstractTriangular, B::AbstractVector) = lmul!(A, copyto!(C, B)) -mul!(C::AbstractMatrix , A::AbstractTriangular, B::AbstractVecOrMat) = lmul!(A, copyto!(C, B)) -mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = lmul!(A, copyto!(C, B)) -mul!(C::AbstractVector , adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVector) = - (A = adjA.parent; lmul!(adjoint(A), copyto!(C, B))) -mul!(C::AbstractMatrix , adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = adjA.parent; lmul!(adjoint(A), copyto!(C, B))) -mul!(C::AbstractVecOrMat, adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = adjA.parent; lmul!(adjoint(A), copyto!(C, B))) -mul!(C::AbstractVector , transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVector) = - (A = transA.parent; lmul!(transpose(A), copyto!(C, B))) -mul!(C::AbstractMatrix , transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = transA.parent; lmul!(transpose(A), copyto!(C, B))) -mul!(C::AbstractVecOrMat, transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) = - (A = transA.parent; lmul!(transpose(A), copyto!(C, B))) -mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) -mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) -mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) -mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B)) -mul!(C::AbstractVector, A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B))) -mul!(C::AbstractVector, A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B))) - -for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), - (:UnitLowerTriangular, 'L', 'U'), - (:UpperTriangular, 'U', 'N'), - (:UnitUpperTriangular, 'U', 'U')) + +for (t, memlay, memlaytrans, uploc, isunitc) in ((:LowerTriangular, :LowerTriangularLayout, :UpperTriangularLayout, 'L', 'N'), + (:UnitLowerTriangular, :UnitLowerTriangularLayout, :UnitUpperTriangularLayout, 'L', 'U'), + (:UpperTriangular, :UpperTriangularLayout, :LowerTriangularLayout, 'U', 'N'), + (:UnitUpperTriangular, :UnitUpperTriangularLayout, :UnitLowerTriangularLayout, 'U', 'U')) @eval begin # Vector multiplication - lmul!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} = - BLAS.trmv!($uploc, 'N', $isunitc, A.data, b) - lmul!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasFloat} = - (A = transA.parent; BLAS.trmv!($uploc, 'T', $isunitc, A.data, b)) - lmul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasReal} = - (A = adjA.parent; BLAS.trmv!($uploc, 'T', $isunitc, A.data, b)) - lmul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasComplex} = - (A = adjA.parent; BLAS.trmv!($uploc, 'C', $isunitc, A.data, b)) + _mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}, + ::AbstractStridedLayout, ::$memlay{<:AbstractColumnMajor}, ::AbstractStridedLayout) where T = + lmul!(A, copyto!(y, b)) + _mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}, + ::AbstractStridedLayout, ::$memlay{<:AbstractRowMajor}, ::AbstractStridedLayout) where T = + lmul!(A, copyto!(y, b)) + _mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}, + ::AbstractStridedLayout, ::$memlay{<:ConjLayout{<:AbstractRowMajor}}, ::AbstractStridedLayout) where T = + lmul!(A, copyto!(y, b)) + + _lmul!(A::AbstractMatrix{T}, b::AbstractVector{T}, ::$memlay{DC}, ::AbstractStridedLayout) where {T<:BlasFloat,DC<:AbstractColumnMajor} = + BLAS.trmv!($uploc, 'N', $isunitc, triangulardata(A), b) + _lmul!(transA::AbstractMatrix{T}, b::AbstractVector{T}, ::$memlaytrans{DR}, ::AbstractStridedLayout) where {T<:BlasFloat,DR<:AbstractRowMajor} = + (A = transpose(transA); BLAS.trmv!($uploc, 'T', $isunitc, triangulardata(A), b)) + _lmul!(adjA::AbstractMatrix{T}, b::AbstractVector{T}, ::$memlaytrans{ConjLayout{DR}}, ::AbstractStridedLayout) where {T<:BlasComplex,DR<:AbstractRowMajor} = + (A = adjoint(adjA); BLAS.trmv!($uploc, 'C', $isunitc, triangulardata(A), b)) # Matrix multiplication - lmul!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} = - BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B) - rmul!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} = + _mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, ::AbstractStridedLayout, ::$memlay, ::AbstractStridedLayout) where {T<:BlasFloat} = + lmul!(A, copyto!(C, B)) + _mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, ::AbstractStridedLayout, ::AbstractStridedLayout, ::$memlay) where {T<:BlasFloat} = + rmul!(copyto!(C, A), B) + + _lmul!(A::AbstractMatrix{T}, B::AbstractMatrix{T}, ::$memlay{DC}, ::AbstractColumnMajor) where {T<:BlasFloat,DC<:AbstractColumnMajor} = + BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), triangulardata(A), B) + _rmul!(A::AbstractMatrix{T}, B::AbstractMatrix{T}, ::AbstractColumnMajor, ::$memlay{DC}) where {T<:BlasFloat,DC<:AbstractColumnMajor} = BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A) - lmul!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasFloat} = - (A = transA.parent; BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B)) - lmul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasComplex} = - (A = adjA.parent; BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B)) - lmul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasReal} = - (A = adjA.parent; BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B)) + _lmul!(transA::AbstractMatrix{T}, B::AbstractMatrix{T}, ::$memlaytrans{DR}, ::AbstractColumnMajor) where {T<:BlasFloat,DR<:AbstractRowMajor} = + (A = transpose(transA); BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), triangulardata(A), B)) + _lmul!(adjA::AbstractMatrix{T}, B::AbstractMatrix{T}, ::$memlaytrans{ConjLayout{DR}}, ::AbstractColumnMajor) where {T<:BlasComplex,DR<:AbstractRowMajor} = + (A = adjoint(adjA); BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), triangulardata(A), B)) - rmul!(A::StridedMatrix{T}, transB::Transpose{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasFloat} = - (B = transB.parent; BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)) - rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasComplex} = - (B = adjB.parent; BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)) - rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasReal} = - (B = adjB.parent; BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)) + _rmul!(A::AbstractMatrix{T}, transB::AbstractMatrix{T}, ::AbstractColumnMajor, ::$memlaytrans{DR}) where {T<:BlasFloat,DR<:AbstractRowMajor} = + (B = transpose(transB); BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)) + _rmul!(A::AbstractMatrix{T}, adjB::AbstractMatrix{T}, ::AbstractColumnMajor, ::$memlaytrans{ConjLayout{DR}}) where {T<:BlasComplex,DR<:AbstractRowMajor} = + (B = adjoint(adjB); BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)) # Left division ldiv!(A::$t{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = @@ -712,16 +793,17 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular), end ## Generic triangular multiplication -function lmul!(A::UpperTriangular, B::StridedVecOrMat) +function _lmul!(A::AbstractMatrix, B::AbstractVecOrMat, ::UpperTriangularLayout{<:AbstractColumnMajor}, _) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = 1:m - Bij = A.data[i,i]*B[i,j] + Bij = Adata[i,i]*B[i,j] for k = i + 1:m - Bij += A.data[i,k]*B[k,j] + Bij += Adata[i,k]*B[k,j] end B[i,j] = Bij end @@ -729,16 +811,17 @@ function lmul!(A::UpperTriangular, B::StridedVecOrMat) B end -function lmul!(A::UnitUpperTriangular, B::StridedVecOrMat) +function _lmul!(A::AbstractMatrix, B::AbstractVecOrMat, ::UnitUpperTriangularLayout{<:AbstractColumnMajor}, _) where T m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = 1:m Bij = B[i,j] for k = i + 1:m - Bij += A.data[i,k]*B[k,j] + Bij += Adata[i,k]*B[k,j] end B[i,j] = Bij end @@ -746,32 +829,34 @@ function lmul!(A::UnitUpperTriangular, B::StridedVecOrMat) B end -function lmul!(A::LowerTriangular, B::StridedVecOrMat) +function _lmul!(A::AbstractMatrix, B::AbstractVecOrMat, ::LowerTriangularLayout{<:AbstractColumnMajor}, _) where T m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = m:-1:1 - Bij = A.data[i,i]*B[i,j] + Bij = Adata[i,i]*B[i,j] for k = 1:i - 1 - Bij += A.data[i,k]*B[k,j] + Bij += Adata[i,k]*B[k,j] end B[i,j] = Bij end end B end -function lmul!(A::UnitLowerTriangular, B::StridedVecOrMat) +function _lmul!(A::AbstractMatrix, B::AbstractVecOrMat, ::UnitLowerTriangularLayout{<:AbstractColumnMajor}, _) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = m:-1:1 Bij = B[i,j] for k = 1:i - 1 - Bij += A.data[i,k]*B[k,j] + Bij += Adata[i,k]*B[k,j] end B[i,j] = Bij end @@ -779,17 +864,18 @@ function lmul!(A::UnitLowerTriangular, B::StridedVecOrMat) B end -function lmul!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat) - A = adjA.parent +function _lmul!(adjA::AbstractMatrix, B::AbstractVecOrMat, ::LowerTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}, _) + A = adjoint(adjA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = m:-1:1 - Bij = A.data[i,i]'B[i,j] + Bij = Adata[i,i]'B[i,j] for k = 1:i - 1 - Bij += A.data[k,i]'B[k,j] + Bij += Adata[k,i]'B[k,j] end B[i,j] = Bij end @@ -797,17 +883,18 @@ function lmul!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat) B end -function lmul!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) - A = adjA.parent +function _lmul!(adjA::AbstractMatrix, B::AbstractVecOrMat, ::UnitLowerTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}, _) + A = adjoint(adjA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = m:-1:1 Bij = B[i,j] for k = 1:i - 1 - Bij += A.data[k,i]'B[k,j] + Bij += Adata[k,i]'B[k,j] end B[i,j] = Bij end @@ -815,34 +902,36 @@ function lmul!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) B end -function lmul!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat) - A = adjA.parent +function _lmul!(adjA::AbstractMatrix, B::AbstractVecOrMat, ::UpperTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}, _) + A = adjoint(adjA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = 1:m - Bij = A.data[i,i]'B[i,j] + Bij = Adata[i,i]'B[i,j] for k = i + 1:m - Bij += A.data[k,i]'B[k,j] + Bij += Adata[k,i]'B[k,j] end B[i,j] = Bij end end B end -function lmul!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) - A = adjA.parent +function _lmul!(adjA::AbstractMatrix, B::AbstractVecOrMat, ::UnitUpperTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}, _) + A = adjoint(adjA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = 1:m Bij = B[i,j] for k = i + 1:m - Bij += A.data[k,i]'B[k,j] + Bij += Adata[k,i]'B[k,j] end B[i,j] = Bij end @@ -850,34 +939,36 @@ function lmul!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) B end -function lmul!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat) - A = transA.parent +function _lmul!(transA::AbstractMatrix, B::AbstractVecOrMat, ::LowerTriangularLayout{<:AbstractRowMajor}, _) + A = transpose(transA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = m:-1:1 - Bij = transpose(A.data[i,i]) * B[i,j] + Bij = transpose(Adata[i,i]) * B[i,j] for k = 1:i - 1 - Bij += transpose(A.data[k,i]) * B[k,j] + Bij += transpose(Adata[k,i]) * B[k,j] end B[i,j] = Bij end end B end -function lmul!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat) - A = transA.parent +function _lmul!(transA::AbstractMatrix, B::AbstractVecOrMat, ::UnitLowerTriangularLayout{<:AbstractRowMajor}, _) + A = transpose(transA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = m:-1:1 Bij = B[i,j] for k = 1:i - 1 - Bij += transpose(A.data[k,i]) * B[k,j] + Bij += transpose(Adata[k,i]) * B[k,j] end B[i,j] = Bij end @@ -885,34 +976,36 @@ function lmul!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMa B end -function lmul!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat) - A = transA.parent +function _lmul!(transA::AbstractMatrix, B::AbstractVecOrMat, ::UpperTriangularLayout{<:AbstractRowMajor}, _) + A = transpose(transA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = 1:m - Bij = transpose(A.data[i,i]) * B[i,j] + Bij = transpose(Adata[i,i]) * B[i,j] for k = i + 1:m - Bij += transpose(A.data[k,i]) * B[k,j] + Bij += transpose(Adata[k,i]) * B[k,j] end B[i,j] = Bij end end B end -function lmul!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat) - A = transA.parent +function _lmul!(transA::AbstractMatrix, B::AbstractVecOrMat, ::UnitUpperTriangularLayout{<:AbstractRowMajor}, _) + A = transpose(transA) m, n = size(B, 1), size(B, 2) if m != size(A, 1) throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) end + Adata = triangulardata(A) for j = 1:n for i = 1:m Bij = B[i,j] for k = i + 1:m - Bij += transpose(A.data[k,i]) * B[k,j] + Bij += transpose(Adata[k,i]) * B[k,j] end B[i,j] = Bij end @@ -920,32 +1013,34 @@ function lmul!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMa B end -function rmul!(A::StridedMatrix, B::UpperTriangular) +function _rmul!(A::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::UpperTriangularLayout) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = n:-1:1 Aij = A[i,j]*B[j,j] for k = 1:j - 1 - Aij += A[i,k]*B.data[k,j] + Aij += A[i,k]*Bdata[k,j] end A[i,j] = Aij end end A end -function rmul!(A::StridedMatrix, B::UnitUpperTriangular) +function _rmul!(A::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::UnitUpperTriangularLayout) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = n:-1:1 Aij = A[i,j] for k = 1:j - 1 - Aij += A[i,k]*B.data[k,j] + Aij += A[i,k]*Bdata[k,j] end A[i,j] = Aij end @@ -953,32 +1048,34 @@ function rmul!(A::StridedMatrix, B::UnitUpperTriangular) A end -function rmul!(A::StridedMatrix, B::LowerTriangular) +function _rmul!(A::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::LowerTriangularLayout) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = 1:n Aij = A[i,j]*B[j,j] for k = j + 1:n - Aij += A[i,k]*B.data[k,j] + Aij += A[i,k]*Bdata[k,j] end A[i,j] = Aij end end A end -function rmul!(A::StridedMatrix, B::UnitLowerTriangular) +function _rmul!(A::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::UnitLowerTriangularLayout) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = 1:n Aij = A[i,j] for k = j + 1:n - Aij += A[i,k]*B.data[k,j] + Aij += A[i,k]*Bdata[k,j] end A[i,j] = Aij end @@ -986,17 +1083,18 @@ function rmul!(A::StridedMatrix, B::UnitLowerTriangular) A end -function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular}) - B = adjB.parent +function _rmul!(A::AbstractMatrix, adjB::AbstractMatrix, ::AbstractStridedLayout, ::LowerTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}) + B = adjoint(adjB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = 1:n - Aij = A[i,j]*B.data[j,j]' + Aij = A[i,j]*Bdata[j,j]' for k = j + 1:n - Aij += A[i,k]*B.data[j,k]' + Aij += A[i,k]*Bdata[j,k]' end A[i,j] = Aij end @@ -1004,17 +1102,18 @@ function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular}) A end -function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular}) - B = adjB.parent +function _rmul!(A::AbstractMatrix, adjB::AbstractMatrix, ::AbstractStridedLayout, ::UnitLowerTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}) + B = adjoint(adjB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = 1:n Aij = A[i,j] for k = j + 1:n - Aij += A[i,k]*B.data[j,k]' + Aij += A[i,k]*Bdata[j,k]' end A[i,j] = Aij end @@ -1022,17 +1121,18 @@ function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular}) A end -function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular}) - B = adjB.parent +function _rmul!(A::AbstractMatrix, adjB::AbstractMatrix, ::AbstractStridedLayout, ::UpperTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}) + B = adjoint(adjB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = n:-1:1 - Aij = A[i,j]*B.data[j,j]' + Aij = A[i,j]*Bdata[j,j]' for k = 1:j - 1 - Aij += A[i,k]*B.data[j,k]' + Aij += A[i,k]*Bdata[j,k]' end A[i,j] = Aij end @@ -1040,17 +1140,18 @@ function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular}) A end -function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular}) - B = adjB.parent +function _rmul!(A::AbstractMatrix, adjB::AbstractMatrix, ::AbstractStridedLayout, ::UnitUpperTriangularLayout{<:ConjLayout{<:AbstractRowMajor}}) + B = adjoint(adjB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = n:-1:1 Aij = A[i,j] for k = 1:j - 1 - Aij += A[i,k]*B.data[j,k]' + Aij += A[i,k]*Bdata[j,k]' end A[i,j] = Aij end @@ -1058,34 +1159,36 @@ function rmul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular}) A end -function rmul!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular}) - B = transB.parent +function _rmul!(A::AbstractMatrix, transB::AbstractMatrix, ::AbstractStridedLayout, ::LowerTriangularLayout{<:AbstractRowMajor}) + B = transpose(transB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = 1:n - Aij = A[i,j] * transpose(B.data[j,j]) + Aij = A[i,j] * transpose(Bdata[j,j]) for k = j + 1:n - Aij += A[i,k] * transpose(B.data[j,k]) + Aij += A[i,k] * transpose(Bdata[j,k]) end A[i,j] = Aij end end A end -function rmul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular}) - B = transB.parent +function _rmul!(A::AbstractMatrix, transB::AbstractMatrix, ::AbstractStridedLayout, ::UnitLowerTriangularLayout{<:AbstractRowMajor}) + B = transpose(transB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = 1:n Aij = A[i,j] for k = j + 1:n - Aij += A[i,k] * transpose(B.data[j,k]) + Aij += A[i,k] * transpose(Bdata[j,k]) end A[i,j] = Aij end @@ -1093,17 +1196,18 @@ function rmul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular}) A end -function rmul!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular}) - B = transB.parent +function _rmul!(A::AbstractMatrix, transB::AbstractMatrix, ::AbstractStridedLayout, ::UpperTriangularLayout{<:AbstractRowMajor}) + B = transpose(transB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = n:-1:1 - Aij = A[i,j] * transpose(B.data[j,j]) + Aij = A[i,j] * transpose(Bdata[j,j]) for k = 1:j - 1 - Aij += A[i,k] * transpose(B.data[j,k]) + Aij += A[i,k] * transpose(Bdata[j,k]) end A[i,j] = Aij end @@ -1111,17 +1215,18 @@ function rmul!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular}) A end -function rmul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular}) - B = transB.parent +function _rmul!(A::AbstractMatrix, transB::AbstractMatrix, ::AbstractStridedLayout, ::UnitUpperTriangularLayout{<:AbstractRowMajor}) + B = transpose(transB) m, n = size(A) if size(B, 1) != n throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))")) end + Bdata = triangulardata(B) for i = 1:m for j = n:-1:1 Aij = A[i,j] for k = 1:j - 1 - Aij += A[i,k] * transpose(B.data[j,k]) + Aij += A[i,k] * transpose(Bdata[j,k]) end A[i,j] = Aij end @@ -1589,7 +1694,27 @@ rdiv!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUp ## for these cases, but I'm not sure it is worth it. (*)(A::Union{Tridiagonal,SymTridiagonal}, B::AbstractTriangular) = rmul!(Matrix(A), B) -for (f, f2!) in ((:*, :lmul!), (:\, :ldiv!)) + +function _mul(A::AbstractMatrix, B::AbstractMatrix, ::Union{LowerTriangularLayout,UnitLowerTriangularLayout}, ::LowerTriangularLayout) + TAB = typeof(*(zero(eltype(A)), zero(eltype(B))) + + *(zero(eltype(A)), zero(eltype(B)))) + BB = similar(B, TAB, size(B)) + copyto!(BB, B) + return LowerTriangular(lmul!(convert(AbstractMatrix{TAB}, A), BB)) +end + + +function _mul(A::AbstractMatrix, B::AbstractMatrix, ::Union{UpperTriangularLayout,UnitUpperTriangularLayout}, ::UpperTriangularLayout) + TAB = typeof(*(zero(eltype(A)), zero(eltype(B))) + + *(zero(eltype(A)), zero(eltype(B)))) + BB = similar(B, TAB, size(B)) + copyto!(BB, B) + return UpperTriangular(lmul!(convert(AbstractMatrix{TAB}, A), BB)) +end + + + +for (f, f2!) in ((:\, :ldiv!),) @eval begin function ($f)(A::LowerTriangular, B::LowerTriangular) TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) + @@ -1626,8 +1751,6 @@ for (f, f2!) in ((:*, :lmul!), (:\, :ldiv!)) end for (ipop, op, xformtype, xformop) in ( - (:lmul!, :*, :Adjoint, :adjoint), - (:lmul!, :*, :Transpose, :transpose), (:ldiv!, :\, :Adjoint, :adjoint), (:ldiv!, :\, :Transpose, :transpose)) @eval begin @@ -1699,8 +1822,6 @@ function (/)(A::UpperTriangular, B::UnitUpperTriangular) end for (ipop, op, xformtype, xformop) in ( - (:rmul!, :*, :Adjoint, :adjoint), - (:rmul!, :*, :Transpose, :transpose), (:rdiv!, :/, :Adjoint, :adjoint), (:rdiv!, :/, :Transpose, :transpose)) @eval begin @@ -1743,66 +1864,22 @@ for (ipop, op, xformtype, xformop) in ( end ## The general promotion methods -function *(A::AbstractTriangular, B::AbstractTriangular) +function _mul(A::AbstractMatrix, B::AbstractMatrix, ::AbstractTriangularLayout, ::AbstractTriangularLayout) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) lmul!(convert(AbstractArray{TAB}, A), BB) end -function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractTriangular) - A = adjA.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - BB = similar(B, TAB, size(B)) - copyto!(BB, B) - lmul!(adjoint(convert(AbstractArray{TAB}, A)), BB) -end -function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractTriangular) - A = transA.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - BB = similar(B, TAB, size(B)) - copyto!(BB, B) - lmul!(transpose(convert(AbstractArray{TAB}, A)), BB) -end - -function *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractTriangular}) - B = adjB.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - AA = similar(A, TAB, size(A)) - copyto!(AA, A) - rmul!(AA, adjoint(convert(AbstractArray{TAB}, B))) -end -function *(A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractTriangular}) - B = transB.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - AA = similar(A, TAB, size(A)) - copyto!(AA, A) - rmul!(AA, transpose(convert(AbstractArray{TAB}, B))) -end for mat in (:AbstractVector, :AbstractMatrix) ### Multiplication with triangle to the left and hence rhs cannot be transposed. - @eval begin - function *(A::AbstractTriangular, B::$mat) - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - BB = similar(B, TAB, size(B)) - copyto!(BB, B) - lmul!(convert(AbstractArray{TAB}, A), BB) - end - function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::$mat) - A = adjA.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - BB = similar(B, TAB, size(B)) - copyto!(BB, B) - lmul!(adjoint(convert(AbstractArray{TAB}, A)), BB) - end - function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::$mat) - A = transA.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - BB = similar(B, TAB, size(B)) - copyto!(BB, B) - lmul!(transpose(convert(AbstractArray{TAB}, A)), BB) - end + @eval function _mul(A::AbstractMatrix, B::$mat, ::AbstractTriangularLayout, ::AbstractStridedLayout) + TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) + BB = similar(B, TAB, size(B)) + copyto!(BB, B) + lmul!(convert(AbstractArray{TAB}, A), BB) end + ### Left division with triangle to the left hence rhs cannot be transposed. No quotients. @eval begin function \(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::$mat) @@ -1898,48 +1975,12 @@ for mat in (:AbstractVector, :AbstractMatrix) end ### Multiplication with triangle to the right and hence lhs cannot be transposed. # Only for AbstractMatrix, hence outside the above loop. -function *(A::AbstractMatrix, B::AbstractTriangular) +function _mul(A::AbstractMatrix, B::AbstractMatrix, ::AbstractStridedLayout, ::AbstractTriangularLayout) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) rmul!(AA, convert(AbstractArray{TAB}, B)) end -function *(A::AbstractMatrix, adjB::Adjoint{<:Any,<:AbstractTriangular}) - B = adjB.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - AA = similar(A, TAB, size(A)) - copyto!(AA, A) - rmul!(AA, adjoint(convert(AbstractArray{TAB}, B))) -end -function *(A::AbstractMatrix, transB::Transpose{<:Any,<:AbstractTriangular}) - B = transB.parent - TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) - AA = similar(A, TAB, size(A)) - copyto!(AA, A) - rmul!(AA, transpose(convert(AbstractArray{TAB}, B))) -end -# ambiguity resolution with definitions in linalg/rowvector.jl -*(v::AdjointAbsVec, A::AbstractTriangular) = adjoint(adjoint(A) * v.parent) -*(v::TransposeAbsVec, A::AbstractTriangular) = transpose(transpose(A) * v.parent) -*(v::AdjointAbsVec, A::Adjoint{<:Any,<:AbstractTriangular}) = adjoint(A.parent * v.parent) -*(v::TransposeAbsVec, A::Transpose{<:Any,<:AbstractTriangular}) = transpose(A.parent * v.parent) - - -# If these are not defined, they will fallback to the versions in matmul.jl -# and dispatch to generic_matmatmul! which is very costly to compile. The methods -# below might compute an unnecessary copy. Eliminating the copy requires adding -# all the promotion logic here once again. Since these methods are probably relatively -# rare, we chose not to bother for now. -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::AbstractTriangular) = copy(A) * B -*(A::Transpose{<:Any,<:AbstractMatrix}, B::AbstractTriangular) = copy(A) * B -*(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::AbstractTriangular, B::Transpose{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractTriangular}) = A * copy(B) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractTriangular}) = copy(A) * B -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractTriangular}) = A * copy(B) -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:AbstractTriangular}) = copy(A) * B # Complex matrix power for upper triangular factor, see: # Higham and Lin, "A Schur-Padé algorithm for fractional powers of a Matrix", @@ -2443,25 +2484,6 @@ end factorize(A::AbstractTriangular) = A -# dismabiguation methods: *(AbstractTriangular, Adj/Trans of AbstractVector) -*(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractVector}) = adjoint(adjoint(B) * adjoint(A)) -*(A::AbstractTriangular, B::Transpose{<:Any,<:AbstractVector}) = transpose(transpose(B) * transpose(A)) -# dismabiguation methods: *(Adj/Trans of AbstractTriangular, Trans/Ajd of AbstractTriangular) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractTriangular}) = copy(A) * B -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractTriangular}) = copy(A) * B -# dismabiguation methods: *(Adj/Trans of AbstractTriangular, Adj/Trans of AbsVec or AbsMat) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVector}) = adjoint(adjoint(B) * adjoint(A)) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractMatrix}) = A * copy(B) -*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVector}) = transpose(transpose(B) * transpose(A)) -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVector}) = transpose(transpose(B) * transpose(A)) -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVector}) = adjoint(adjoint(B) * adjoint(A)) -*(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractMatrix}) = A * copy(B) -# dismabiguation methods: *(Adj/Trans of AbsVec or AbsMat, Adj/Trans of AbstractTriangular) -*(A::Adjoint{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:AbstractTriangular}) = adjoint(adjoint(B) * adjoint(A)) -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:AbstractTriangular}) = copy(A) * B -*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractTriangular}) = transpose(transpose(B) * transpose(A)) -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractTriangular}) = copy(A) * B - # disambiguation methods: /(Adjoint of AbsVec, <:AbstractTriangular) /(u::AdjointAbsVec, A::Union{LowerTriangular,UpperTriangular}) = adjoint(adjoint(A) \ u.parent) /(u::AdjointAbsVec, A::Union{UnitLowerTriangular,UnitUpperTriangular}) = adjoint(adjoint(A) \ u.parent) diff --git a/stdlib/LinearAlgebra/test/adjtrans.jl b/stdlib/LinearAlgebra/test/adjtrans.jl index 9dd5b068a5dcd..7e7cb60f312ee 100644 --- a/stdlib/LinearAlgebra/test/adjtrans.jl +++ b/stdlib/LinearAlgebra/test/adjtrans.jl @@ -447,6 +447,42 @@ end @test adjoint!(b, a) === b end +@testset "adjoint and transpose MemoryLayout" begin + A = [1.0 2; 3 4] + @test Base.MemoryLayout(A') == Base.DenseRowMajor() + @test Base.MemoryLayout(transpose(A)) == Base.DenseRowMajor() + B = [1.0+im 2; 3 4] + @test Base.MemoryLayout(B') == LinearAlgebra.ConjLayout(Base.DenseRowMajor()) + @test Base.MemoryLayout(transpose(B)) == Base.DenseRowMajor() + VA = view(A, 1:1, 1:1) + @test Base.MemoryLayout(VA') == Base.RowMajor() + @test Base.MemoryLayout(transpose(VA)) == Base.RowMajor() + VB = view(B, 1:1, 1:1) + @test Base.MemoryLayout(VB') == LinearAlgebra.ConjLayout(Base.RowMajor()) + @test Base.MemoryLayout(transpose(VB)) == Base.RowMajor() + VA = view(A, 1:2:2, 1:2:2) + @test Base.MemoryLayout(VA') == Base.StridedLayout() + @test Base.MemoryLayout(transpose(VA)) == Base.StridedLayout() + VB = view(B, 1:2:2, 1:2:2) + @test Base.MemoryLayout(VB') == LinearAlgebra.ConjLayout(Base.StridedLayout()) + @test Base.MemoryLayout(transpose(VB)) == Base.StridedLayout() + VA2 = view(A, [1,2], :) + @test Base.MemoryLayout(VA2') == Base.UnknownLayout() + @test Base.MemoryLayout(transpose(VA2)) == Base.UnknownLayout() + VB2 = view(B, [1,2], :) + @test Base.MemoryLayout(VB2') == Base.UnknownLayout() + @test Base.MemoryLayout(transpose(VB2)) == Base.UnknownLayout() + VAc = view(A', 1:1, 1:1) + @test Base.MemoryLayout(VAc) == Base.RowMajor() + VAt = view(transpose(A), 1:1, 1:1) + @test Base.MemoryLayout(VAt) == Base.RowMajor() + VBc = view(B', 1:1, 1:1) + @test Base.MemoryLayout(VBc) == LinearAlgebra.ConjLayout(Base.RowMajor()) + VBt = view(transpose(B), 1:1, 1:1) + @test Base.MemoryLayout(VBt) == Base.RowMajor() +end + + @testset "aliasing with adjoint and transpose" begin A = collect(reshape(1:25, 5, 5)) .+ rand.().*im B = copy(A) diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl index 9bbd914047462..178e4ebdf40da 100644 --- a/stdlib/LinearAlgebra/test/dense.jl +++ b/stdlib/LinearAlgebra/test/dense.jl @@ -847,6 +847,25 @@ end end end + +@testset "Dispatch to BLAS routines" begin + A = rand(100,100) + x = rand(100) + @test all(A*x .=== view(A,:,:)*x .=== view(A',:,:)'*x .=== + BLAS.gemv!('N', 1.0, A, x, 0.0, similar(x))) + @test all(A'*x .=== view(A',:,:)*x .=== view(A,:,:)'*x .=== + BLAS.gemv!('T', 1.0, A, x, 0.0, similar(x))) + + B = rand(ComplexF64,100,100) + y = rand(ComplexF64,100) + @test all(B*y .=== view(B,:,:)*y .=== view(B',:,:)'*y .=== transpose(view(transpose(B),:,:))*y .=== + BLAS.gemv!('N', one(ComplexF64), B, y, zero(ComplexF64), similar(y))) + @test all(B'*y .=== view(B',:,:)*y .=== view(B,:,:)'*y .=== + BLAS.gemv!('C', one(ComplexF64), B, y, zero(ComplexF64), similar(y))) + @test all(transpose(B)*y .=== view(transpose(B),:,:)*y .=== transpose(view(B,:,:))*y .=== + BLAS.gemv!('T', one(ComplexF64), B, y, zero(ComplexF64), similar(y))) +end + @testset "inverse of Adjoint" begin A = randn(n, n) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 5f3d09561bcbc..9955ae5e0fe16 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -250,6 +250,11 @@ vecdot_(x,y) = invoke(vecdot, Tuple{Any,Any}, x,y) end end +@testset "array vecdot" begin + a = rand(3,3,3) + @test vecdot(a,a) === vecdot(vec(a), a) === vecdot(a, vec(a)) === dot(vec(a), vec(a)) === BLAS.dot(a,a) +end + @testset "Issue 11978" begin A = Matrix{Matrix{Float64}}(undef, 2, 2) A[1,1] = Matrix(1.0I, 3, 3) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 2a6a279a397a0..2c05ccb13db04 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -487,6 +487,81 @@ end end end +@testset "Symmetric/Hermitian MemoryLayout" begin + A = [1.0 2; 3 4] + @test Base.MemoryLayout(Symmetric(A)) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Hermitian(A)) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Transpose(Symmetric(A))) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Transpose(Hermitian(A))) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Adjoint(Symmetric(A))) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Adjoint(Hermitian(A))) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(view(Symmetric(A),:,:)) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(view(Hermitian(A),:,:)) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Symmetric(A')) == LinearAlgebra.SymmetricLayout(Base.DenseRowMajor(),'U') + @test Base.MemoryLayout(Hermitian(A')) == LinearAlgebra.SymmetricLayout(Base.DenseRowMajor(),'U') + @test Base.MemoryLayout(Symmetric(transpose(A))) == LinearAlgebra.SymmetricLayout(Base.DenseRowMajor(),'U') + @test Base.MemoryLayout(Hermitian(transpose(A))) == LinearAlgebra.SymmetricLayout(Base.DenseRowMajor(),'U') + + @test LinearAlgebra.symmetricdata(Symmetric(A)) ≡ A + @test LinearAlgebra.symmetricdata(Hermitian(A)) ≡ A + @test LinearAlgebra.symmetricdata(Transpose(Symmetric(A))) ≡ A + @test LinearAlgebra.symmetricdata(Transpose(Hermitian(A))) ≡ A + @test LinearAlgebra.symmetricdata(Adjoint(Symmetric(A))) ≡ A + @test LinearAlgebra.symmetricdata(Adjoint(Hermitian(A))) ≡ A + @test LinearAlgebra.symmetricdata(view(Symmetric(A),:,:)) ≡ A + @test LinearAlgebra.symmetricdata(view(Hermitian(A),:,:)) ≡ A + @test LinearAlgebra.symmetricdata(Symmetric(A')) ≡ A' + @test LinearAlgebra.symmetricdata(Hermitian(A')) ≡ A' + @test LinearAlgebra.symmetricdata(Symmetric(transpose(A))) ≡ transpose(A) + @test LinearAlgebra.symmetricdata(Hermitian(transpose(A))) ≡ transpose(A) + + B = [1.0+im 2; 3 4] + @test Base.MemoryLayout(Symmetric(B)) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Hermitian(B)) == LinearAlgebra.HermitianLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Transpose(Symmetric(B))) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Transpose(Hermitian(B))) == Base.UnknownLayout() + @test Base.MemoryLayout(Adjoint(Symmetric(B))) == Base.UnknownLayout() + @test Base.MemoryLayout(Adjoint(Hermitian(B))) == LinearAlgebra.HermitianLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(view(Symmetric(B),:,:)) == LinearAlgebra.SymmetricLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(view(Hermitian(B),:,:)) == LinearAlgebra.HermitianLayout(Base.DenseColumnMajor(),'U') + @test Base.MemoryLayout(Symmetric(B')) == Base.UnknownLayout() + @test Base.MemoryLayout(Hermitian(B')) == Base.UnknownLayout() + @test Base.MemoryLayout(Symmetric(transpose(B))) == LinearAlgebra.SymmetricLayout(Base.DenseRowMajor(),'U') + @test Base.MemoryLayout(Hermitian(transpose(B))) == LinearAlgebra.HermitianLayout(Base.DenseRowMajor(),'U') + + @test LinearAlgebra.symmetricdata(Symmetric(B)) ≡ B + @test LinearAlgebra.hermitiandata(Hermitian(B)) ≡ B + @test LinearAlgebra.symmetricdata(Transpose(Symmetric(B))) ≡ B + @test LinearAlgebra.hermitiandata(Adjoint(Hermitian(B))) ≡ B + @test LinearAlgebra.symmetricdata(view(Symmetric(B),:,:)) ≡ B + @test LinearAlgebra.hermitiandata(view(Hermitian(B),:,:)) ≡ B + @test LinearAlgebra.symmetricdata(Symmetric(B')) ≡ B' + @test LinearAlgebra.hermitiandata(Hermitian(B')) ≡ B' + @test LinearAlgebra.symmetricdata(Symmetric(transpose(B))) ≡ transpose(B) + @test LinearAlgebra.hermitiandata(Hermitian(transpose(B))) ≡ transpose(B) + + A = randn(100,100) + x = randn(100) + @test all(Symmetric(A)*x .=== Hermitian(A)*x .=== Symmetric(A)'*x .=== + Symmetric(view(A,:,:)',:L)*x .=== view(Symmetric(A),:,:)*x .=== + BLAS.symv!('U', 1.0, A, x, 0.0, similar(x))) + @test all(Symmetric(A,:L)*x .=== Hermitian(A,:L)*x .=== Symmetric(A,:L)'*x .=== + Symmetric(view(A,:,:)',:U)*x .=== view(Hermitian(A,:L),:,:)*x .=== + BLAS.symv!('L', 1.0, A, x, 0.0, similar(x))) + A = randn(ComplexF64,100,100) + x = randn(ComplexF64,100) + @test all(Symmetric(A)*x .=== transpose(Symmetric(A))*x .=== + Symmetric(transpose(view(A,:,:)),:L)*x .=== view(Symmetric(A),:,:)*x .=== + BLAS.symv!('U', one(ComplexF64), A, x, zero(ComplexF64), similar(x))) + @test all(Symmetric(A,:L)*x .=== transpose(Symmetric(A,:L))*x .=== + Symmetric(transpose(view(A,:,:)),:U)*x .=== view(Symmetric(A,:L),:,:)*x .=== + BLAS.symv!('L', one(ComplexF64), A, x, zero(ComplexF64), similar(x))) + @test all(Hermitian(A)*x .=== Hermitian(A)'*x .=== view(Hermitian(A),:,:)*x .=== + BLAS.hemv!('U', one(ComplexF64), A, x, zero(ComplexF64), similar(x))) + @test all(Hermitian(A,:L)*x .=== Hermitian(A,:L)'*x .=== view(Hermitian(A,:L),:,:)*x .=== + BLAS.hemv!('L', one(ComplexF64), A, x, zero(ComplexF64), similar(x))) +end + @testset "#25625 recursive transposition" begin A = Matrix{Matrix{Int}}(undef, 2, 2) A[1,1] = [1 2; 2 3] @@ -521,4 +596,4 @@ end @test Hermitian(A, :U)[1,1] == Hermitian(A, :L)[1,1] == real(A[1,1]) end -end # module TestSymmetric +end # module TestSymmetricx diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 50b94cc6f7764..bffb68c5b9116 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -549,6 +549,93 @@ end sprint(show, MIME"text/plain"(), UnitUpperTriangular(2ones(Int64,3,3)))) end +@testset "triangular MemoryLayout" begin + A = [1.0 2; 3 4] + B = [1.0+im 2; 3 4] + for (TriType, TriLayout, TriLayoutTrans) in ((UpperTriangular, LinearAlgebra.UpperTriangularLayout, LinearAlgebra.LowerTriangularLayout), + (UnitUpperTriangular, LinearAlgebra.UnitUpperTriangularLayout, LinearAlgebra.UnitLowerTriangularLayout), + (LowerTriangular, LinearAlgebra.LowerTriangularLayout, LinearAlgebra.UpperTriangularLayout), + (UnitLowerTriangular, LinearAlgebra.UnitLowerTriangularLayout, LinearAlgebra.UnitUpperTriangularLayout)) + @test Base.MemoryLayout(TriType(A)) == TriLayout(Base.DenseColumnMajor()) + @test Base.MemoryLayout(TriType(transpose(A))) == Base.UnknownLayout() + @test Base.MemoryLayout(TriType(A')) == Base.UnknownLayout() + @test Base.MemoryLayout(transpose(TriType(A))) == TriLayoutTrans(Base.DenseRowMajor()) + @test Base.MemoryLayout(TriType(A)') == TriLayoutTrans(Base.DenseRowMajor()) + + @test Base.MemoryLayout(TriType(B)) == TriLayout(Base.DenseColumnMajor()) + @test Base.MemoryLayout(TriType(transpose(B))) == Base.UnknownLayout() + @test Base.MemoryLayout(TriType(B')) == Base.UnknownLayout() + @test Base.MemoryLayout(transpose(TriType(B))) == TriLayoutTrans(Base.DenseRowMajor()) + @test Base.MemoryLayout(TriType(B)') == TriLayoutTrans(LinearAlgebra.ConjLayout(Base.DenseRowMajor())) + end + + A = randn(Float64, 100, 100) + x = randn(Float64, 100) + @test all(UpperTriangular(A)*x .=== UpperTriangular(view(A',:,:)')*x .=== + view(UpperTriangular(A), Base.OneTo(100), :)*x .=== + BLAS.trmv!('U', 'N', 'N', A, copy(x))) + @test all(UnitUpperTriangular(A)*x .=== UnitUpperTriangular(view(A',:,:)')*x .=== + view(UnitUpperTriangular(A), Base.OneTo(100), :)*x .=== + BLAS.trmv!('U', 'N', 'U', A, copy(x))) + @test all(LowerTriangular(A)*x .=== LowerTriangular(view(A',:,:)')*x .=== + view(LowerTriangular(A), Base.OneTo(100), :)*x .=== + BLAS.trmv!('L', 'N', 'N', A, copy(x))) + @test all(UnitLowerTriangular(A)*x .=== UnitLowerTriangular(view(A',:,:)')*x .=== + view(UnitLowerTriangular(A), Base.OneTo(100), :)*x .=== + BLAS.trmv!('L', 'N', 'U', A, copy(x))) + @test all(UpperTriangular(A)'*x .=== UpperTriangular(view(A',:,:)')'*x .=== + view(UpperTriangular(A)', Base.OneTo(100), :)*x .=== + BLAS.trmv!('U', 'T', 'N', A, copy(x))) + @test all(UnitUpperTriangular(A)'*x .=== UnitUpperTriangular(view(A',:,:)')'*x .=== + view(UnitUpperTriangular(A)', Base.OneTo(100), :)*x .=== + BLAS.trmv!('U', 'T', 'U', A, copy(x))) + @test all(LowerTriangular(A)'*x .=== LowerTriangular(view(A',:,:)')'*x .=== + view(LowerTriangular(A)', Base.OneTo(100), :)*x .=== + BLAS.trmv!('L', 'T', 'N', A, copy(x))) + @test all(UnitLowerTriangular(A)'*x .=== UnitLowerTriangular(view(A',:,:)')'*x .=== + view(UnitLowerTriangular(A)', Base.OneTo(100), :)*x .=== + BLAS.trmv!('L', 'T', 'U', A, copy(x))) + + A = randn(ComplexF64, 100, 100) + x = randn(ComplexF64, 100) + @test all(UpperTriangular(A)*x .=== UpperTriangular(view(A',:,:)')*x .=== + view(UpperTriangular(A), :, Base.OneTo(100))*x .=== + BLAS.trmv!('U', 'N', 'N', A, copy(x))) + @test all(UnitUpperTriangular(A)*x .=== UnitUpperTriangular(view(A',:,:)')*x .=== + view(UnitUpperTriangular(A), :, Base.OneTo(100))*x .=== + BLAS.trmv!('U', 'N', 'U', A, copy(x))) + @test all(LowerTriangular(A)*x .=== LowerTriangular(view(A',:,:)')*x .=== + view(LowerTriangular(A), :, Base.OneTo(100))*x .=== + BLAS.trmv!('L', 'N', 'N', A, copy(x))) + @test all(UnitLowerTriangular(A)*x .=== UnitLowerTriangular(view(A',:,:)')*x .=== + view(UnitLowerTriangular(A), :, Base.OneTo(100))*x .=== + BLAS.trmv!('L', 'N', 'U', A, copy(x))) + @test all(transpose(UpperTriangular(A))*x .=== transpose(UpperTriangular(view(A',:,:)'))*x .=== + view(transpose(UpperTriangular(A)), :, Base.OneTo(100))*x .=== + BLAS.trmv!('U', 'T', 'N', A, copy(x))) + @test all(transpose(UnitUpperTriangular(A))*x .=== transpose(UnitUpperTriangular(view(A',:,:)'))*x .=== + view(transpose(UnitUpperTriangular(A)), :, Base.OneTo(100))*x .=== + BLAS.trmv!('U', 'T', 'U', A, copy(x))) + @test all(transpose(LowerTriangular(A))*x .=== transpose(LowerTriangular(view(A',:,:)'))*x .=== + view(transpose(LowerTriangular(A)), :, Base.OneTo(100))*x .=== + BLAS.trmv!('L', 'T', 'N', A, copy(x))) + @test all(transpose(UnitLowerTriangular(A))*x .=== transpose(UnitLowerTriangular(view(A',:,:)'))*x .=== + view(transpose(UnitLowerTriangular(A)), :, Base.OneTo(100))*x .=== + BLAS.trmv!('L', 'T', 'U', A, copy(x))) + @test all(UpperTriangular(A)'*x .=== UpperTriangular(view(A',:,:)')'*x .=== + view(UpperTriangular(A)', :, Base.OneTo(100))*x .=== + BLAS.trmv!('U', 'C', 'N', A, copy(x))) + @test all(UnitUpperTriangular(A)'*x .=== UnitUpperTriangular(view(A',:,:)')'*x .=== + view(UnitUpperTriangular(A)', :, Base.OneTo(100))*x .=== + BLAS.trmv!('U', 'C', 'U', A, copy(x))) + @test all(LowerTriangular(A)'*x .=== LowerTriangular(view(A',:,:)')'*x .=== + view(LowerTriangular(A)', :, Base.OneTo(100))*x .=== + BLAS.trmv!('L', 'C', 'N', A, copy(x))) + @test all(UnitLowerTriangular(A)'*x .=== UnitLowerTriangular(view(A',:,:)')'*x .=== + view(UnitLowerTriangular(A)', :, Base.OneTo(100))*x .=== + BLAS.trmv!('L', 'C', 'U', A, copy(x))) +end + @testset "adjoint/transpose triangular/vector multiplication" begin for elty in (Float64, ComplexF64), trity in (UpperTriangular, LowerTriangular) A1 = trity(rand(elty, 1, 1)) diff --git a/stdlib/SharedArrays/src/SharedArrays.jl b/stdlib/SharedArrays/src/SharedArrays.jl index 6874c4dc27235..dc42aa00db2ba 100644 --- a/stdlib/SharedArrays/src/SharedArrays.jl +++ b/stdlib/SharedArrays/src/SharedArrays.jl @@ -10,7 +10,8 @@ module SharedArrays using Mmap, Distributed, Random import Base: length, size, ndims, IndexStyle, reshape, convert, deepcopy_internal, - show, getindex, setindex!, fill!, similar, reduce, map!, copyto!, unsafe_convert + show, getindex, setindex!, fill!, similar, reduce, map!, copyto!, unsafe_convert, + strides, stride import Random using Serialization using Serialization: serialize_cycle_header, serialize_type, writetag, UNDEFREF_TAG, serialize, deserialize @@ -342,6 +343,9 @@ for each worker process. """ localindices(S::SharedArray) = S.pidx > 0 ? range_1dim(S, S.pidx) : 1:0 +strides(S::SharedArray) = strides(sdata(S)) +stride(S::SharedArray, i::Int) = stride(sdata(S), i) + unsafe_convert(::Type{Ptr{T}}, S::SharedArray{T}) where {T} = unsafe_convert(Ptr{T}, sdata(S)) unsafe_convert(::Type{Ptr{T}}, S::SharedArray ) where {T} = unsafe_convert(Ptr{T}, sdata(S)) diff --git a/test/abstractarray.jl b/test/abstractarray.jl index 2183e4013a77a..442562e43fca9 100644 --- a/test/abstractarray.jl +++ b/test/abstractarray.jl @@ -894,3 +894,106 @@ end @test hcat(1:2, fill(1, (2,1))) == hcat([1:2;], fill(1, (2,1))) == reshape([1,2,1,1],2,2) @test [(1:3) (4:6); fill(1, (3,2))] == reshape([1,2,3,1,1,1,4,5,6,1,1,1], 6,2) end + +struct MyDenseArray{T,N} <: DenseArray{T,N} end +struct RowMajorArray{T,N} <: AbstractArray{T,N} + A::Array{T,N} +end + +Base.getindex(A::RowMajorArray{T,N}, inds::Vararg{Int,N}) where {T,N} = A.A[reverse(inds)...] +Base.setindex!(A::RowMajorArray, v, inds::Vararg{Int,N}) where {T,N} = (A.A[reverse(inds)...] = v) +Base.MemoryLayout(::RowMajorArray) = Base.DenseRowMajor() +Base.size(A::RowMajorArray) = reverse(size(A.A)) + +@testset "MemoryLayout for Array, SubArray, and ReinterpretArray" begin + let A = randn(2,2) + @test Base.MemoryLayout(A) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:,:)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:,1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:,1:1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1:1,1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,2)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,1:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1,:)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,:,1:-2:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:1,1:-2:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:1,1:2)) == Base.ColumnMajor() + @test Base.MemoryLayout(view(A,1:1,:)) == Base.ColumnMajor() + @test Base.MemoryLayout(view(A,1:2:1,1:2:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:2:1,:)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:2:1,1:2)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,[1,2],:)) == Base.UnknownLayout() + @test Base.MemoryLayout(view(A,:,[1,2])) == Base.UnknownLayout() + + @test Base.MemoryLayout(Base.ReshapedArray(A,(4,),())) == Base.DenseColumnMajor() + @test Base.MemoryLayout(Base.ReshapedArray(view(A,:,:),(4,),())) == Base.DenseColumnMajor() + @test Base.MemoryLayout(Base.ReshapedArray(view(A,:,1),(1,2),())) == Base.DenseColumnMajor() + + @test Base.MemoryLayout(reinterpret(ComplexF64,A)) == Base.DenseColumnMajor() + end + + let A = randn(2,2,2) + @test Base.MemoryLayout(A) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:,:,:)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:,:,1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,:,:,1:1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1:1,1,1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,2,2)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,1:1,1:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1,:,:)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:1,1:2,1:2)) == Base.ColumnMajor() + @test Base.MemoryLayout(view(A,1:1,:,:)) == Base.ColumnMajor() + @test Base.MemoryLayout(view(A,1:2:1,1:2:1,1:2:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:2:1,:,:)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,[1,2],:,:)) == Base.UnknownLayout() + + @test Base.MemoryLayout(Base.ReshapedArray(A,(3^3,),())) == Base.DenseColumnMajor() + @test Base.MemoryLayout(Base.ReshapedArray(view(A,:,:,:),(3^3,),())) == Base.DenseColumnMajor() + @test Base.MemoryLayout(Base.ReshapedArray(view(A,:,1,1),(1,3),())) == Base.DenseColumnMajor() + + @test Base.MemoryLayout(reinterpret(ComplexF64,A)) == Base.DenseColumnMajor() + end + + let A = RowMajorArray(randn(2,2)) + @test Base.MemoryLayout(A) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,:,:)) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,:)) == Base.UnknownLayout() + @test Base.MemoryLayout(view(A,:,1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,:,1:1)) == Base.RowMajor() + @test Base.MemoryLayout(view(A,1:1,1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1,2)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,1:1)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,:)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1:1,1:2)) == Base.RowMajor() + @test Base.MemoryLayout(view(A,1:1,:)) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,1:2:1,1:2:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:-2:1,:)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:-2:1,1:2)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,[1,2],:)) == Base.UnknownLayout() + @test Base.MemoryLayout(view(A,:,[1,2])) == Base.UnknownLayout() + end + + let A = RowMajorArray(randn(2,2,2)) + @test Base.MemoryLayout(A) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,:,:,:)) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,:)) == Base.UnknownLayout() + @test Base.MemoryLayout(view(A,:,:,1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,:,:,1:1)) == Base.RowMajor() + @test Base.MemoryLayout(view(A,1:1,1,1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1,2,2)) == Base.DenseColumnMajor() + @test Base.MemoryLayout(view(A,1,1:1,1:1)) == Base.RowMajor() + @test Base.MemoryLayout(view(A,1,:,:)) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,1:1,1:2,1:2)) == Base.RowMajor() + @test Base.MemoryLayout(view(A,1:1,:,:)) == Base.DenseRowMajor() + @test Base.MemoryLayout(view(A,1:2:1,1:2:1,1:2:1)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,1:2:1,:,:)) == Base.StridedLayout() + @test Base.MemoryLayout(view(A,[1,2],:,:)) == Base.UnknownLayout() + end + + @test Base.MemoryLayout(MyDenseArray{Float64,1}()) == Base.DenseColumnMajor() + @test Base.MemoryLayout(MyDenseArray{Float64,3}()) == Base.DenseColumnMajor() + + @test Base.MemoryLayout(BitArray([true,true,false])) == Base.UnknownLayout() +end