diff --git a/Project.toml b/Project.toml index 4db785c..124a6e1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ArrayLayouts" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" authors = ["Sheehan Olver "] -version = "0.3.4" +version = "0.3.5" [deps] FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" diff --git a/src/ArrayLayouts.jl b/src/ArrayLayouts.jl index 2e33bb8..a2ec30b 100644 --- a/src/ArrayLayouts.jl +++ b/src/ArrayLayouts.jl @@ -104,15 +104,8 @@ include("factorizations.jl") @inline layout_getindex(A, I...) = sub_materialize(view(A, I...)) - - -macro layoutmatrix(Typ) +macro _layoutgetindex(Typ) esc(quote - ArrayLayouts.@layoutldiv $Typ - ArrayLayouts.@layoutmul $Typ - ArrayLayouts.@layoutlmul $Typ - ArrayLayouts.@layoutfactorizations $Typ - @inline Base.getindex(A::$Typ, kr::Colon, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr) @inline Base.getindex(A::$Typ, kr::Colon, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr) @inline Base.getindex(A::$Typ, kr::AbstractUnitRange, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr) @@ -123,6 +116,24 @@ macro layoutmatrix(Typ) end) end +macro layoutgetindex(Typ) + esc(quote + ArrayLayouts.@_layoutgetindex $Typ + ArrayLayouts.@_layoutgetindex LinearAlgebra.AbstractTriangular{<:Any,<:$Typ} + end) +end + + +macro layoutmatrix(Typ) + esc(quote + ArrayLayouts.@layoutldiv $Typ + ArrayLayouts.@layoutmul $Typ + ArrayLayouts.@layoutlmul $Typ + ArrayLayouts.@layoutfactorizations $Typ + ArrayLayouts.@layoutgetindex $Typ + end) +end + @layoutmatrix LayoutMatrix getindex(A::LayoutVector, kr::AbstractVector) = layout_getindex(A, kr) diff --git a/src/diagonal.jl b/src/diagonal.jl index 2aa7d40..336407e 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1,7 +1,4 @@ -rowsupport(::DiagonalLayout, _, k) = isempty(k) ? (1:0) : minimum(k):maximum(k) -colsupport(::DiagonalLayout, _, j) = isempty(j) ? (1:0) : minimum(j):maximum(j) - ### # Lmul #### diff --git a/src/memorylayout.jl b/src/memorylayout.jl index 689bc77..b8224b5 100644 --- a/src/memorylayout.jl +++ b/src/memorylayout.jl @@ -500,32 +500,47 @@ abstract type AbstractBandedLayout <: MemoryLayout end abstract type AbstractTridiagonalLayout <: AbstractBandedLayout end struct DiagonalLayout{ML} <: AbstractBandedLayout end +struct BidiagonalLayout{ML} <: AbstractBandedLayout end struct SymTridiagonalLayout{ML} <: AbstractTridiagonalLayout end struct TridiagonalLayout{ML} <: AbstractTridiagonalLayout end diagonallayout(_) = DiagonalLayout{UnknownLayout}() diagonallayout(::ML) where ML<:AbstractStridedLayout = DiagonalLayout{ML}() MemoryLayout(D::Type{Diagonal{T,P}}) where {T,P} = diagonallayout(MemoryLayout(P)) -diagonaldata(D::Diagonal) = parent(D) - +MemoryLayout(::Type{Bidiagonal{T,V}}) where {T,V} = BidiagonalLayout{typeof(MemoryLayout(V))}() MemoryLayout(::Type{SymTridiagonal{T,P}}) where {T,P} = SymTridiagonalLayout{typeof(MemoryLayout(P))}() +MemoryLayout(::Type{Tridiagonal{T,P}}) where {T,P} = TridiagonalLayout{typeof(MemoryLayout(P))}() + +bidiagonaluplo(A::Bidiagonal) = A.uplo +bidiagonaluplo(A::AdjOrTrans) = bidiagonaluplo(parent(A)) == 'L' ? 'U' : 'L' + +diagonaldata(D::Diagonal) = parent(D) +diagonaldata(D::Bidiagonal) = D.dv diagonaldata(D::SymTridiagonal) = D.dv +diagonaldata(D::Tridiagonal) = D.d + +supdiagonaldata(D::Bidiagonal) = D.uplo == 'U' ? D.ev : throw(ArgumentError("$D is lower-bidiagonal")) +subdiagonaldata(D::Bidiagonal) = D.uplo == 'L' ? D.ev : throw(ArgumentError("$D is upper-bidiagonal")) + supdiagonaldata(D::SymTridiagonal) = D.ev subdiagonaldata(D::SymTridiagonal) = D.ev -MemoryLayout(::Type{Tridiagonal{T,P}}) where {T,P} = TridiagonalLayout{typeof(MemoryLayout(P))}() -diagonaldata(D::Tridiagonal) = D.d subdiagonaldata(D::Tridiagonal) = D.dl supdiagonaldata(D::Tridiagonal) = D.du - transposelayout(ml::DiagonalLayout) = ml +transposelayout(ml::BidiagonalLayout) = ml transposelayout(ml::SymTridiagonalLayout) = ml transposelayout(ml::TridiagonalLayout) = ml transposelayout(ml::ConjLayout{DiagonalLayout}) = ml adjointlayout(::Type{<:Real}, ml::SymTridiagonalLayout) = ml adjointlayout(::Type{<:Real}, ml::TridiagonalLayout) = ml +adjointlayout(::Type{<:Real}, ml::BidiagonalLayout) = ml + +symmetriclayout(B::BidiagonalLayout{ML}) where ML = SymTridiagonalLayout{ML}() +hermitianlayout(::Type{<:Real}, B::BidiagonalLayout{ML}) where ML = SymTridiagonalLayout{ML}() +hermitianlayout(_, B::BidiagonalLayout) = HermitianLayout{typeof(B)}() subdiagonaldata(D::Transpose) = supdiagonaldata(parent(D)) supdiagonaldata(D::Transpose) = subdiagonaldata(parent(D)) @@ -535,6 +550,10 @@ subdiagonaldata(D::Adjoint{<:Real}) = supdiagonaldata(parent(D)) supdiagonaldata(D::Adjoint{<:Real}) = subdiagonaldata(parent(D)) diagonaldata(D::Adjoint{<:Real}) = diagonaldata(parent(D)) +diagonaldata(S::HermOrSym) = diagonaldata(parent(S)) +subdiagonaldata(S::HermOrSym) = symmetricuplo(S) == 'L' ? subdiagonaldata(parent(S)) : supdiagonaldata(parent(S)) +supdiagonaldata(S::HermOrSym) = symmetricuplo(S) == 'L' ? subdiagonaldata(parent(S)) : supdiagonaldata(parent(S)) + ### # Fill #### @@ -575,3 +594,20 @@ colsupport(A) = colsupport(A, axes(A,2)) rowsupport(::ZerosLayout, A, _) = 1:0 colsupport(::ZerosLayout, A, _) = 1:0 + +rowsupport(::DiagonalLayout, _, k) = isempty(k) ? (1:0) : minimum(k):maximum(k) +colsupport(::DiagonalLayout, _, j) = isempty(j) ? (1:0) : minimum(j):maximum(j) + +colsupport(::BidiagonalLayout, A, j) = + bidiagonaluplo(A) == 'L' ? (minimum(j):min(size(A,1),maximum(j)+1)) : (max(minimum(j)-1,1):maximum(j)) +rowsupport(::BidiagonalLayout, A, j) = + bidiagonaluplo(A) == 'U' ? (minimum(j):min(size(A,2),maximum(j)+1)) : (max(minimum(j)-1,1):maximum(j)) + +colsupport(::AbstractTridiagonalLayout, A, j) = max(minimum(j)-1,1):min(size(A,1),maximum(j)+1) +rowsupport(::AbstractTridiagonalLayout, A, j) = max(minimum(j)-1,1):min(size(A,2),maximum(j)+1) + +colsupport(::SymmetricLayout, A, j) = first(colsupport(symmetricdata(A),j)):last(rowsupport(symmetricdata(A),j)) +rowsupport(::SymmetricLayout, A, j) = colsupport(A, j) + +colsupport(::HermitianLayout, A, j) = first(colsupport(hermitiandata(A),j)):last(rowsupport(hermitiandata(A),j)) +rowsupport(::HermitianLayout, A, j) = colsupport(A, j) diff --git a/src/muladd.jl b/src/muladd.jl index 580149f..391b9c4 100644 --- a/src/muladd.jl +++ b/src/muladd.jl @@ -194,7 +194,7 @@ function default_blasmul!(α, A::AbstractMatrix, B::AbstractMatrix, β, C::Abstr C end -function default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector) +function _default_blasmul!(::IndexLinear, α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector) mA, nA = size(A) mB = length(B) nA == mB || throw(DimensionMismatch("Dimensions must match")) @@ -217,6 +217,30 @@ function default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::Abstr C end +function _default_blasmul!(::IndexCartesian, α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector) + mA, nA = size(A) + mB = length(B) + nA == mB || throw(DimensionMismatch("Dimensions must match")) + length(C) == mA || throw(DimensionMismatch("Dimensions must match")) + + lmul!(β, C) + (nA == 0 || mB == 0) && return C + + z = zero(A[1,1]*B[1] + A[1,1]*B[1]) + + @inbounds for k in colsupport(B,1) + b = B[k] + for i = colsupport(A,k) + C[i] += α * A[i,k] * b + end + end + + C +end + +default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector) = + _default_blasmul!(Base.IndexStyle(typeof(A)), α, A, B, β, C) + function materialize!(M::MatMulMatAdd) α, A, B, β, C = M.α, M.A, M.B, M.β, M.C if C ≡ B diff --git a/src/triangular.jl b/src/triangular.jl index fc86258..04a14b8 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -118,7 +118,7 @@ materialize!(M::MatRmulMat{<:AbstractStridedLayout,<:TriangularLayout}) = Linear @inline function copyto!(dest::AbstractArray, M::Ldiv{<:TriangularLayout}) A, B = M.A, M.B - dest ≡ B || (dest .= B) + dest ≡ B || copyto!(dest, B) ldiv!(A, dest) end diff --git a/test/test_layouts.jl b/test/test_layouts.jl index 9e06203..c4be72d 100644 --- a/test/test_layouts.jl +++ b/test/test_layouts.jl @@ -6,7 +6,7 @@ import ArrayLayouts: MemoryLayout, DenseRowMajor, DenseColumnMajor, StridedLayou UnitLowerTriangularLayout, ScalarLayout, UnknownLayout, hermitiandata, symmetricdata, FillLayout, ZerosLayout, DiagonalLayout, TridiagonalLayout, SymTridiagonalLayout, colsupport, rowsupport, - diagonaldata, subdiagonaldata, supdiagonaldata + diagonaldata, subdiagonaldata, supdiagonaldata, BidiagonalLayout, bidiagonaluplo struct FooBar end struct FooNumber <: Number end @@ -80,7 +80,45 @@ struct FooNumber <: Number end @test MemoryLayout(view(randn(5)',[1,3])) == UnknownLayout() end - @testset "Symmetric/Hermitian MemoryLayout" begin + @testset "Bi/Tridiagonal" begin + T = Tridiagonal(randn(5),randn(6),randn(5)) + S = SymTridiagonal(T.d, T.du) + Bl = Bidiagonal(T.d, T.dl, :L) + Bu = Bidiagonal(T.d, T.du, :U) + + @test MemoryLayout(T) isa TridiagonalLayout + @test MemoryLayout(Adjoint(T)) isa TridiagonalLayout + @test MemoryLayout(Transpose(T)) isa TridiagonalLayout + @test MemoryLayout(S) isa SymTridiagonalLayout + @test MemoryLayout(Adjoint(S)) isa SymTridiagonalLayout + @test MemoryLayout(Transpose(S)) isa SymTridiagonalLayout + @test MemoryLayout(Bl) isa BidiagonalLayout + @test MemoryLayout(Adjoint(Bl)) isa BidiagonalLayout + @test MemoryLayout(Transpose(Bl)) isa BidiagonalLayout + @test MemoryLayout(Bu) isa BidiagonalLayout + @test MemoryLayout(Adjoint(Bu)) isa BidiagonalLayout + @test MemoryLayout(Transpose(Bu)) isa BidiagonalLayout + + @test bidiagonaluplo(Bl) == bidiagonaluplo(Adjoint(Bu)) == 'L' + @test bidiagonaluplo(Bu) == bidiagonaluplo(Adjoint(Bl)) == 'U' + + @test diagonaldata(T) == diagonaldata(T') == diagonaldata(S) == diagonaldata(Bl) == diagonaldata(Bu) + @test supdiagonaldata(T) == subdiagonaldata(Adjoint(T)) == subdiagonaldata(Transpose(T)) == + supdiagonaldata(S) == subdiagonaldata(S) == + supdiagonaldata(Bu) == subdiagonaldata(Adjoint(Bu)) == subdiagonaldata(Transpose(Bu)) + @test subdiagonaldata(T) == supdiagonaldata(Adjoint(T)) == supdiagonaldata(Transpose(T)) == + subdiagonaldata(Bl) == supdiagonaldata(Adjoint(Bl)) == supdiagonaldata(Transpose(Bl)) == + T.dl + + @test colsupport(T,3) == rowsupport(T,3) == colsupport(S,3) == rowsupport(S,3) == 2:4 + @test colsupport(T,3:6) == rowsupport(T,3:6) == colsupport(S,3:6) == rowsupport(S,3:6) == 2:6 + @test colsupport(Bl,3) == rowsupport(Bu,3) == rowsupport(Adjoint(Bl),3) == 3:4 + @test rowsupport(Bl,3) == colsupport(Bu,3) == colsupport(Adjoint(Bl),3) == 2:3 + @test colsupport(Bl,3:6) == rowsupport(Bu,3:6) == 3:6 + @test colsupport(Bu,3:6) == rowsupport(Bl,3:6) == 2:6 + end + + @testset "Symmetric/Hermitian" begin A = [1.0 2; 3 4] @test MemoryLayout(Symmetric(A)) == SymmetricLayout{DenseColumnMajor}() @test MemoryLayout(Hermitian(A)) == SymmetricLayout{DenseColumnMajor}() @@ -108,6 +146,11 @@ struct FooNumber <: Number end @test symmetricdata(Symmetric(transpose(A))) ≡ transpose(A) @test symmetricdata(Hermitian(transpose(A))) ≡ transpose(A) + @test colsupport(Symmetric(A),2) ≡ colsupport(Symmetric(A),1:2) ≡ + rowsupport(Symmetric(A),2) ≡ rowsupport(Symmetric(A),1:2) ≡ 1:2 + @test colsupport(Hermitian(A),2) ≡ colsupport(Hermitian(A),1:2) ≡ + rowsupport(Hermitian(A),2) ≡ rowsupport(Hermitian(A),1:2) ≡ 1:2 + B = [1.0+im 2; 3 4] @test MemoryLayout(Symmetric(B)) == SymmetricLayout{DenseColumnMajor}() @test MemoryLayout(Hermitian(B)) == HermitianLayout{DenseColumnMajor}() @@ -132,6 +175,26 @@ struct FooNumber <: Number end @test hermitiandata(Hermitian(B')) ≡ B' @test symmetricdata(Symmetric(transpose(B))) ≡ transpose(B) @test hermitiandata(Hermitian(transpose(B))) ≡ transpose(B) + + @testset "Bidiagonal" begin + B = Bidiagonal(randn(6),randn(5),:U) + Bc = Bidiagonal(randn(6) .+ 0im,randn(5) .+ 1im,:U) + S = Symmetric(B) + H = Hermitian(B) + Sc = Symmetric(Bc) + Hc = Hermitian(Bc) + + @test MemoryLayout(S) isa SymTridiagonalLayout + @test MemoryLayout(H) isa SymTridiagonalLayout + @test MemoryLayout(Sc) isa SymTridiagonalLayout + @test MemoryLayout(Hc) isa HermitianLayout + + @test diagonaldata(S) == diagonaldata(B) + @test subdiagonaldata(S) == supdiagonaldata(S) == supdiagonaldata(B) + + @test colsupport(S,3) == colsupport(H,3) == colsupport(Sc,3) == colsupport(Hc,3) == 2:4 + @test rowsupport(S,3) == rowsupport(H,3) == rowsupport(Sc,3) == rowsupport(Hc,3) == 2:4 + end end @testset "triangular MemoryLayout" begin @@ -233,19 +296,4 @@ struct FooNumber <: Number end MemoryLayout(revD) @test 0 == @allocated MemoryLayout(revD) end - - @testset "Tridiagonal" begin - T = Tridiagonal(randn(5),randn(6),randn(5)) - S = SymTridiagonal(T.d, T.du) - @test MemoryLayout(T) isa TridiagonalLayout - @test MemoryLayout(Adjoint(T)) isa TridiagonalLayout - @test MemoryLayout(Transpose(T)) isa TridiagonalLayout - @test MemoryLayout(S) isa SymTridiagonalLayout - @test MemoryLayout(Adjoint(S)) isa SymTridiagonalLayout - @test MemoryLayout(Transpose(S)) isa SymTridiagonalLayout - - @test diagonaldata(T) == diagonaldata(T') == diagonaldata(S) - @test supdiagonaldata(T) == subdiagonaldata(Adjoint(T)) == subdiagonaldata(Transpose(T)) == supdiagonaldata(S) == subdiagonaldata(S) - @test subdiagonaldata(T) == supdiagonaldata(Adjoint(T)) == supdiagonaldata(Transpose(T)) == T.dl - end end diff --git a/test/test_muladd.jl b/test/test_muladd.jl index 0fb0a80..d0861f9 100644 --- a/test/test_muladd.jl +++ b/test/test_muladd.jl @@ -2,542 +2,551 @@ using ArrayLayouts, FillArrays, Random, Test import ArrayLayouts: DenseColumnMajor, AbstractStridedLayout, AbstractColumnMajor, DiagonalLayout, mul Random.seed!(0) +@testset "Multiplication" begin + @testset "Matrix * Vector" begin + @testset "eltype" begin + v = MulAdd(1,zeros(Int,2,2), zeros(Float64,2),0,zeros(2)) + A = MulAdd(1,zeros(Int,2,2), zeros(Float64,2,2), 0,zeros(2,2)) + @test eltype(v) == eltype(A) == Float64 + @test @inferred(axes(v)) == (@inferred(axes(v,1)),) == (Base.OneTo(2),) + @test @inferred(size(v)) == (@inferred(size(v,1)),) == (2,) + @test @inferred(axes(A)) == (@inferred(axes(A,1)),@inferred(axes(A,2))) == (Base.OneTo(2),Base.OneTo(2)) + @test @inferred(size(A)) == (@inferred(size(A,1)),size(A,2)) == (2,2) + + à = materialize(A) + @test à isa Matrix{Float64} + @test à == zeros(2,2) + + c = similar(v) + fill!(c,NaN) + @test copyto!(c,v) == zeros(2) + fill!(c,NaN) + c .= v + @test c == zeros(2) + end -@testset "Matrix * Vector" begin - @testset "eltype" begin - v = MulAdd(1,zeros(Int,2,2), zeros(Float64,2),0,zeros(2)) - A = MulAdd(1,zeros(Int,2,2), zeros(Float64,2,2), 0,zeros(2,2)) - @test eltype(v) == eltype(A) == Float64 - @test @inferred(axes(v)) == (@inferred(axes(v,1)),) == (Base.OneTo(2),) - @test @inferred(size(v)) == (@inferred(size(v,1)),) == (2,) - @test @inferred(axes(A)) == (@inferred(axes(A,1)),@inferred(axes(A,2))) == (Base.OneTo(2),Base.OneTo(2)) - @test @inferred(size(A)) == (@inferred(size(A,1)),size(A,2)) == (2,2) - - à = materialize(A) - @test à isa Matrix{Float64} - @test à == zeros(2,2) - - c = similar(v) - fill!(c,NaN) - @test copyto!(c,v) == zeros(2) - fill!(c,NaN) - c .= v - @test c == zeros(2) - end - - @testset "gemv Float64" begin - for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), - view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)), - b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)) - c = similar(b); - - muladd!(1.0,A,b,0.0,c) - @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c))) - copyto!(c, MulAdd(1.0,A,b,0.0,c)) - @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c))) - c .= MulAdd(1.0,A,b,0.0,c) - @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c))) - - M = MulAdd(2.0, A,b, 0.0, c) - @test M isa MulAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractColumnMajor} - c .= M - @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c))) - copyto!(c, M) - @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c))) - - c = copy(b) - M = MulAdd(1.0, A,b, 1.0, c) - @test M isa MulAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractColumnMajor} - copyto!(c, M) - @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b))) + @testset "gemv Float64" begin + for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), + view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)), + b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)) + c = similar(b); + + muladd!(1.0,A,b,0.0,c) + @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c))) + copyto!(c, MulAdd(1.0,A,b,0.0,c)) + @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c))) + c .= MulAdd(1.0,A,b,0.0,c) + @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c))) + + M = MulAdd(2.0, A,b, 0.0, c) + @test M isa MulAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractColumnMajor} + c .= M + @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c))) + copyto!(c, M) + @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c))) + + c = copy(b) + M = MulAdd(1.0, A,b, 1.0, c) + @test M isa MulAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractColumnMajor} + copyto!(c, M) + @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b))) + end end - end - @testset "gemv mixed array types" begin - (A, b, c) = (randn(5,5), randn(5), 1.0:5.0) - d = similar(b) - d .= MulAdd(1.0,A,b,1.0,c) - @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 1.0, Vector{Float64}(c))) - end + @testset "gemv mixed array types" begin + (A, b, c) = (randn(5,5), randn(5), 1.0:5.0) + d = similar(b) + d .= MulAdd(1.0,A,b,1.0,c) + @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 1.0, Vector{Float64}(c))) + end - @testset "gemv Complex" begin - for T in (ComplexF64,), - A in (randn(T,5,5), view(randn(T,5,5),:,:)), - b in (randn(T,5), view(randn(T,5),:)) - c = similar(b); + @testset "gemv Complex" begin + for T in (ComplexF64,), + A in (randn(T,5,5), view(randn(T,5,5),:,:)), + b in (randn(T,5), view(randn(T,5),:)) + c = similar(b); - c .= MulAdd(one(T), A, b, zero(T), c) - @test all(c .=== BLAS.gemv!('N', one(T), A, b, zero(T), similar(c))) + c .= MulAdd(one(T), A, b, zero(T), c) + @test all(c .=== BLAS.gemv!('N', one(T), A, b, zero(T), similar(c))) - c = copy(b) - muladd!(2one(T),A,b,one(T),c) - @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, one(T), copy(b))) + c = copy(b) + muladd!(2one(T),A,b,one(T),c) + @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, one(T), copy(b))) + end end end -end -@testset "Matrix * Matrix" begin - @testset "gemm" begin - for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), - view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)), - B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), - view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)) - C = similar(B); + @testset "Matrix * Matrix" begin + @testset "gemm" begin + for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), + view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)), + B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), + view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)) + C = similar(B); - C .= MulAdd(1.0,A,B,0.0,C) - @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, similar(C))) + C .= MulAdd(1.0,A,B,0.0,C) + @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, similar(C))) - C .= MulAdd(2.0,A,B,0.0,C) - @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 0.0, similar(C))) + C .= MulAdd(2.0,A,B,0.0,C) + @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 0.0, similar(C))) - C = copy(B) - C .= MulAdd(2.0,A,B,1.0,C) - @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, copy(B))) + C = copy(B) + C .= MulAdd(2.0,A,B,1.0,C) + @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, copy(B))) + end end - end - @testset "gemm Complex" begin - for T in (ComplexF64,), - A in (randn(T,5,5), view(randn(T,5,5),:,:)), - B in (randn(T,5,5), view(randn(T,5,5),:,:)) - C = similar(B); + @testset "gemm Complex" begin + for T in (ComplexF64,), + A in (randn(T,5,5), view(randn(T,5,5),:,:)), + B in (randn(T,5,5), view(randn(T,5,5),:,:)) + C = similar(B); - C .= MulAdd(one(T),A,B,zero(T),C) - @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, zero(T), similar(C))) + C .= MulAdd(one(T),A,B,zero(T),C) + @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, zero(T), similar(C))) - C .= MulAdd(2one(T),A,B,zero(T),C) - @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, zero(T), similar(C))) + C .= MulAdd(2one(T),A,B,zero(T),C) + @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, zero(T), similar(C))) - C = copy(B) - C .= MulAdd(one(T),A,B,one(T),C) - @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, one(T), copy(B))) + C = copy(B) + C .= MulAdd(one(T),A,B,one(T),C) + @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, one(T), copy(B))) + end end - end - @testset "gemm mixed array types" begin - (A, B, C) = (randn(5,5), randn(5,5), reshape(1.0:25.0,5,5)) - D = similar(B) - D .= MulAdd(1.0,A,B,1.0, C) - @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, Matrix{Float64}(C))) + @testset "gemm mixed array types" begin + (A, B, C) = (randn(5,5), randn(5,5), reshape(1.0:25.0,5,5)) + D = similar(B) + D .= MulAdd(1.0,A,B,1.0, C) + @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, Matrix{Float64}(C))) + end end -end -@testset "adjtrans" begin - @testset "gemv adjtrans" begin - for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), - view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)), - b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)), - Ac in (transpose(A), A') - c = similar(b); + @testset "adjtrans" begin + @testset "gemv adjtrans" begin + for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:), + view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)), + b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)), + Ac in (transpose(A), A') + c = similar(b); - c .= MulAdd(1.0,Ac,b,0.0,c) - @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 0.0, similar(c))) + c .= MulAdd(1.0,Ac,b,0.0,c) + @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 0.0, similar(c))) - c .= MulAdd(2.0,Ac,b,0.0,c) - @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 0.0, similar(c))) + c .= MulAdd(2.0,Ac,b,0.0,c) + @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 0.0, similar(c))) - c = copy(b) - c .= MulAdd(1.0,Ac,b,1.0,c) - @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 1.0, copy(b))) + c = copy(b) + c .= MulAdd(1.0,Ac,b,1.0,c) + @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 1.0, copy(b))) + end end - end - @testset "gemv adjtrans mixed types" begin - (A, b, c) = (randn(5,5), randn(5), 1.0:5.0) - for Ac in (transpose(A), A') - d = similar(b) - d .= MulAdd(1.0,Ac,b,1.0, c) - @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 1.0, Vector{Float64}(c))) + @testset "gemv adjtrans mixed types" begin + (A, b, c) = (randn(5,5), randn(5), 1.0:5.0) + for Ac in (transpose(A), A') + d = similar(b) + d .= MulAdd(1.0,Ac,b,1.0, c) + @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 1.0, Vector{Float64}(c))) - d .= MulAdd(3.0,Ac,b,2.0, c) - @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, Vector{Float64}(c))) + d .= MulAdd(3.0,Ac,b,2.0, c) + @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, Vector{Float64}(c))) + end end - end - @testset "gemv adjtrans Complex" begin - for T in (ComplexF64,), - A in (randn(T,5,5), view(randn(T,5,5),:,:)), - b in (randn(T,5), view(randn(T,5),:)), - (Ac,trans) in ((transpose(A),'T'), (A','C')) - c = similar(b); + @testset "gemv adjtrans Complex" begin + for T in (ComplexF64,), + A in (randn(T,5,5), view(randn(T,5,5),:,:)), + b in (randn(T,5), view(randn(T,5),:)), + (Ac,trans) in ((transpose(A),'T'), (A','C')) + c = similar(b); - c .= MulAdd(one(T),Ac,b,zero(T),c) - @test all(c .=== BLAS.gemv!(trans, one(T), A, b, zero(T), similar(c))) + c .= MulAdd(one(T),Ac,b,zero(T),c) + @test all(c .=== BLAS.gemv!(trans, one(T), A, b, zero(T), similar(c))) - c = copy(b) - c .= MulAdd(3one(T),Ac,b,2one(T),c) - @test all(c .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b))) - end + c = copy(b) + c .= MulAdd(3one(T),Ac,b,2one(T),c) + @test all(c .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b))) + end - C = randn(6,6) + im*randn(6,6) - V = view(C', 2:3, 3:4) - c = randn(2) + im*randn(2) - @test all(muladd!(1.0+0im,V,c,0.0+0im,similar(c,2)) .=== BLAS.gemv!('C', 1.0+0im, Matrix(V'), c, 0.0+0im, similar(c,2))) - @test all(muladd!(1.0+0im,V',c,0.0+0im,similar(c,2)) .=== BLAS.gemv!('N', 1.0+0im, Matrix(V'), c, 0.0+0im, similar(c,2))) - end + C = randn(6,6) + im*randn(6,6) + V = view(C', 2:3, 3:4) + c = randn(2) + im*randn(2) + @test all(muladd!(1.0+0im,V,c,0.0+0im,similar(c,2)) .=== BLAS.gemv!('C', 1.0+0im, Matrix(V'), c, 0.0+0im, similar(c,2))) + @test all(muladd!(1.0+0im,V',c,0.0+0im,similar(c,2)) .=== BLAS.gemv!('N', 1.0+0im, Matrix(V'), c, 0.0+0im, similar(c,2))) + end - @testset "gemm adjtrans" begin - for A in (randn(5,5), view(randn(5,5),:,:)), - B in (randn(5,5), view(randn(5,5),1:5,:)) - for Ac in (transpose(A), A') - C = copy(B) - C .= MulAdd(3.0, Ac, B, 2.0, C) - @test all(C .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B))) - end - for Bc in (transpose(B), B') - C = copy(B) - C .= MulAdd(3.0, A, Bc, 2.0, C) - @test all(C .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B))) - end - for Ac in (transpose(A), A'), Bc in (transpose(B), B') - C = copy(B) - d = similar(C) - d .= MulAdd(3.0, Ac, Bc, 2.0, C) - @test all(d .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B))) + @testset "gemm adjtrans" begin + for A in (randn(5,5), view(randn(5,5),:,:)), + B in (randn(5,5), view(randn(5,5),1:5,:)) + for Ac in (transpose(A), A') + C = copy(B) + C .= MulAdd(3.0, Ac, B, 2.0, C) + @test all(C .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B))) + end + for Bc in (transpose(B), B') + C = copy(B) + C .= MulAdd(3.0, A, Bc, 2.0, C) + @test all(C .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B))) + end + for Ac in (transpose(A), A'), Bc in (transpose(B), B') + C = copy(B) + d = similar(C) + d .= MulAdd(3.0, Ac, Bc, 2.0, C) + @test all(d .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B))) + end end end - end - @testset "symv adjtrans" begin - A = randn(100,100) - x = randn(100) - for As in (Symmetric(A), Hermitian(A), Symmetric(A)', Symmetric(view(A,:,:)',:L), view(Symmetric(A),:,:)) - @test all( muladd!(1.0,As,x,0.0,similar(x)) .=== BLAS.symv!('U', 1.0, A, x, 0.0, similar(x)) ) - end + @testset "symv adjtrans" begin + A = randn(100,100) + x = randn(100) + for As in (Symmetric(A), Hermitian(A), Symmetric(A)', Symmetric(view(A,:,:)',:L), view(Symmetric(A),:,:)) + @test all( muladd!(1.0,As,x,0.0,similar(x)) .=== BLAS.symv!('U', 1.0, A, x, 0.0, similar(x)) ) + end - for As in (Symmetric(A,:L), Hermitian(A,:L), Symmetric(A,:L)', Symmetric(view(A,:,:)',:U), view(Hermitian(A,:L),:,:)) - @test all( muladd!(1.0,As,x,0.0,similar(x)) .=== BLAS.symv!('L', 1.0, A, x, 0.0, similar(x)) ) - end + for As in (Symmetric(A,:L), Hermitian(A,:L), Symmetric(A,:L)', Symmetric(view(A,:,:)',:U), view(Hermitian(A,:L),:,:)) + @test all( muladd!(1.0,As,x,0.0,similar(x)) .=== BLAS.symv!('L', 1.0, A, x, 0.0, similar(x)) ) + end - T = ComplexF64 - A = randn(T,100,100) - x = randn(T,100) - for As in (Symmetric(A), transpose(Symmetric(A)), Symmetric(transpose(view(A,:,:)),:L), view(Symmetric(A),:,:)) - @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.symv!('U', one(T), A, x, zero(T), similar(x)) ) - end + T = ComplexF64 + A = randn(T,100,100) + x = randn(T,100) + for As in (Symmetric(A), transpose(Symmetric(A)), Symmetric(transpose(view(A,:,:)),:L), view(Symmetric(A),:,:)) + @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.symv!('U', one(T), A, x, zero(T), similar(x)) ) + end - for As in (Symmetric(A,:L), transpose(Symmetric(A,:L)), Symmetric(transpose(view(A,:,:)),:U), view(Symmetric(A,:L),:,:)) - @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.symv!('L', one(T), A, x, zero(T), similar(x)) ) - end + for As in (Symmetric(A,:L), transpose(Symmetric(A,:L)), Symmetric(transpose(view(A,:,:)),:U), view(Symmetric(A,:L),:,:)) + @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.symv!('L', one(T), A, x, zero(T), similar(x)) ) + end - for As in(Hermitian(A), Hermitian(A)', view(Hermitian(A),:,:)) - @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.hemv!('U', one(T), A, x, zero(T), similar(x)) ) - end + for As in(Hermitian(A), Hermitian(A)', view(Hermitian(A),:,:)) + @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.hemv!('U', one(T), A, x, zero(T), similar(x)) ) + end - for As in(Hermitian(A,:L), Hermitian(A,:L)', view(Hermitian(A,:L),:,:)) - @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.hemv!('L', one(T), A, x, zero(T), similar(x)) ) + for As in(Hermitian(A,:L), Hermitian(A,:L)', view(Hermitian(A,:L),:,:)) + @test all( muladd!(one(T),As,x,zero(T),similar(x)) .=== BLAS.hemv!('L', one(T), A, x, zero(T), similar(x)) ) + end end end -end -@testset "Mul" begin - @testset "Mixed types" begin - A = randn(5,6) - b = rand(Int,6) - c = Array{Float64}(undef, 5) - d = similar(c) - copyto!(d, MulAdd(1, A, b, 0.0, d)) - @test d == copyto!(similar(d), MulAdd(1, A, b, 0.0, d)) ≈ A*b - @test copyto!(similar(d), MulAdd(1, A, b, 1.0, d)) ≈ A*b + d + @testset "Mul" begin + @testset "Mixed types" begin + A = randn(5,6) + b = rand(Int,6) + c = Array{Float64}(undef, 5) + d = similar(c) + copyto!(d, MulAdd(1, A, b, 0.0, d)) + @test d == copyto!(similar(d), MulAdd(1, A, b, 0.0, d)) ≈ A*b + @test copyto!(similar(d), MulAdd(1, A, b, 1.0, d)) ≈ A*b + d - @test all((similar(d) .= MulAdd(1, A, b, 1.0, d)) .=== copyto!(similar(d), MulAdd(1, A, b, 1.0, d))) + @test all((similar(d) .= MulAdd(1, A, b, 1.0, d)) .=== copyto!(similar(d), MulAdd(1, A, b, 1.0, d))) - B = rand(Int,6,4) - C = Array{Float64}(undef, 5, 4) + B = rand(Int,6,4) + C = Array{Float64}(undef, 5, 4) - @test_throws DimensionMismatch muladd!(1.0,A,B,0.0,similar(C,4,4)) + @test_throws DimensionMismatch muladd!(1.0,A,B,0.0,similar(C,4,4)) - A = randn(Float64,20,22) - B = randn(ComplexF64,22,24) - C = similar(B,20,24) - @test copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)) ≈ A*B - end - - @testset "no allocation" begin - A = randn(5,5); x = randn(5); y = randn(5); c = similar(y); - muladd!(2.0, A, x, 3.0, y) - @test @allocated(muladd!(2.0, A, x, 3.0, y)) == 0 - VA = view(A,:,1:2) - vx = view(x,1:2) - vy = view(y,:) - muladd!(2.0, VA, vx, 3.0, vy) - @test @allocated(muladd!(2.0, VA, vx, 3.0, vy)) == 0 - end - - @testset "BigFloat" begin - A = BigFloat.(randn(5,5)) - x = BigFloat.(randn(5)) - @test A*x == muladd!(1.0,A,x,0.0,copy(x)) - @test_throws UndefRefError muladd!(1.0,A,x,0.0,similar(x)) - end + A = randn(Float64,20,22) + B = randn(ComplexF64,22,24) + C = similar(B,20,24) + @test copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)) ≈ A*B + end - @testset "Scalar * Vector" begin - A, x = [1 2; 3 4] , [[1,2],[3,4]] - @test muladd!(1.0,A,x,0.0, 1.0*x) == A*x - end + @testset "no allocation" begin + A = randn(5,5); x = randn(5); y = randn(5); c = similar(y); + muladd!(2.0, A, x, 3.0, y) + @test @allocated(muladd!(2.0, A, x, 3.0, y)) == 0 + VA = view(A,:,1:2) + vx = view(x,1:2) + vy = view(y,:) + muladd!(2.0, VA, vx, 3.0, vy) + @test @allocated(muladd!(2.0, VA, vx, 3.0, vy)) == 0 + end - @testset "adjoint" begin - A = randn(5,5) - V = view(A',1:3,1:3) - b = randn(3) - @test MemoryLayout(V) isa RowMajor - @test MemoryLayout(V') isa ColumnMajor - @test strides(V) == (5,1) - @test strides(V') == (1,5) - @test mul(V, b) ≈ V*b ≈ Matrix(V)*b - @test mul(V', b) ≈ V'b ≈ Matrix(V)'*b - end -end - -@testset "Lmul/Rmul" begin - @testset "tri Lmul" begin - @testset "Float * Float vector" begin - A = randn(Float64, 100, 100) - x = randn(Float64, 100) - - L = Lmul(UpperTriangular(A),x) - @test size(L) == (size(L,1),) == (100,) - @test axes(L) == (axes(L,1),) == (Base.OneTo(100),) - @test eltype(L) == Float64 - @test length(L) == 100 - - @test similar(L) isa Vector{Float64} - @test similar(L,Int) isa Vector{Int} - - @test all(copy(Lmul(UpperTriangular(A),x)) .=== - ArrayLayouts.lmul!(UpperTriangular(A),copy(x)) .=== - copyto!(similar(x),Lmul(UpperTriangular(A),x)) .=== - UpperTriangular(A)*x .=== - BLAS.trmv!('U', 'N', 'N', A, copy(x))) - @test all(copyto!(similar(x),Lmul(UnitUpperTriangular(A),x)) .=== - UnitUpperTriangular(A)*x .=== - BLAS.trmv!('U', 'N', 'U', A, copy(x))) - @test all(copyto!(similar(x),Lmul(LowerTriangular(A),x)) .=== - LowerTriangular(A)*x .=== - BLAS.trmv!('L', 'N', 'N', A, copy(x))) - @test all(ArrayLayouts.lmul!(UnitLowerTriangular(A),copy(x)) .=== - UnitLowerTriangular(A)x .=== - BLAS.trmv!('L', 'N', 'U', A, copy(x))) - - @test all(ArrayLayouts.lmul!(UpperTriangular(A)',copy(x)) .=== - ArrayLayouts.lmul!(LowerTriangular(A'),copy(x)) .=== - UpperTriangular(A)'x .=== - BLAS.trmv!('U', 'T', 'N', A, copy(x))) - @test all(ArrayLayouts.lmul!(UnitUpperTriangular(A)',copy(x)) .=== - ArrayLayouts.lmul!(UnitLowerTriangular(A'),copy(x)) .=== - UnitUpperTriangular(A)'*x .=== - BLAS.trmv!('U', 'T', 'U', A, copy(x))) - @test all(ArrayLayouts.lmul!(LowerTriangular(A)',copy(x)) .=== - ArrayLayouts.lmul!(UpperTriangular(A'),copy(x)) .=== - LowerTriangular(A)'*x .=== - BLAS.trmv!('L', 'T', 'N', A, copy(x))) - @test all(ArrayLayouts.lmul!(UnitLowerTriangular(A)',copy(x)) .=== - ArrayLayouts.lmul!(UnitUpperTriangular(A'),copy(x)) .=== - UnitLowerTriangular(A)'*x .=== - BLAS.trmv!('L', 'T', 'U', A, copy(x))) + @testset "BigFloat" begin + A = BigFloat.(randn(5,5)) + x = BigFloat.(randn(5)) + @test A*x == muladd!(1.0,A,x,0.0,copy(x)) + @test_throws UndefRefError muladd!(1.0,A,x,0.0,similar(x)) end - @testset "Float * Complex vector" begin - T = ComplexF64 - A = randn(T, 100, 100) - x = randn(T, 100) - - @test all((y = copy(x); y .= Lmul(UpperTriangular(A),y) ) .=== - UpperTriangular(A)x .=== - BLAS.trmv!('U', 'N', 'N', A, copy(x))) - @test all((y = copy(x); y .= Lmul(UnitUpperTriangular(A),y) ) .=== - UnitUpperTriangular(A)x .=== - BLAS.trmv!('U', 'N', 'U', A, copy(x))) - @test all((y = copy(x); y .= Lmul(LowerTriangular(A),y) ) .=== - LowerTriangular(A)x .=== - BLAS.trmv!('L', 'N', 'N', A, copy(x))) - @test all((y = copy(x); y .= Lmul(UnitLowerTriangular(A),y) ) .=== - UnitLowerTriangular(A)x .=== - BLAS.trmv!('L', 'N', 'U', A, copy(x))) - @test LowerTriangular(A') == UpperTriangular(A)' - - @test all((y = copy(x); y .= Lmul(transpose(UpperTriangular(A)),y) ) .=== - (similar(x) .= Lmul(LowerTriangular(transpose(A)),x)) .=== - transpose(UpperTriangular(A))x .=== - BLAS.trmv!('U', 'T', 'N', A, copy(x))) - @test all((y = copy(x); y .= Lmul(transpose(UnitUpperTriangular(A)),y) ) .=== - (similar(x) .= Lmul(UnitLowerTriangular(transpose(A)),x)) .=== - transpose(UnitUpperTriangular(A))x .=== - BLAS.trmv!('U', 'T', 'U', A, copy(x))) - @test all((y = copy(x); y .= Lmul(transpose(LowerTriangular(A)),y) ) .=== - (similar(x) .= Lmul(UpperTriangular(transpose(A)),x)) .=== - transpose(LowerTriangular(A))x .=== - BLAS.trmv!('L', 'T', 'N', A, copy(x))) - @test all((y = copy(x); y .= Lmul(transpose(UnitLowerTriangular(A)),y) ) .=== - (similar(x) .= Lmul(UnitUpperTriangular(transpose(A)),x)) .=== - transpose(UnitLowerTriangular(A))x .=== - BLAS.trmv!('L', 'T', 'U', A, copy(x))) - - @test all((y = copy(x); y .= Lmul(UpperTriangular(A)',y) ) .=== - (similar(x) .= Lmul(LowerTriangular(A'),x)) .=== - UpperTriangular(A)'x .=== - BLAS.trmv!('U', 'C', 'N', A, copy(x))) - @test all((y = copy(x); y .= Lmul(UnitUpperTriangular(A)',y) ) .=== - (similar(x) .= Lmul(UnitLowerTriangular(A'),x)) .=== - UnitUpperTriangular(A)'x .=== - BLAS.trmv!('U', 'C', 'U', A, copy(x))) - @test all((y = copy(x); y .= Lmul(LowerTriangular(A)',y) ) .=== - (similar(x) .= Lmul(UpperTriangular(A'),x)) .=== - LowerTriangular(A)'x .=== - BLAS.trmv!('L', 'C', 'N', A, copy(x))) - @test all((y = copy(x); y .= Lmul(UnitLowerTriangular(A)',y) ) .=== - (similar(x) .= Lmul(UnitUpperTriangular(A'),x)) .=== - UnitLowerTriangular(A)'x .=== - BLAS.trmv!('L', 'C', 'U', A, copy(x))) + @testset "Scalar * Vector" begin + A, x = [1 2; 3 4] , [[1,2],[3,4]] + @test muladd!(1.0,A,x,0.0, 1.0*x) == A*x + end + + @testset "adjoint" begin + A = randn(5,5) + V = view(A',1:3,1:3) + b = randn(3) + @test MemoryLayout(V) isa RowMajor + @test MemoryLayout(V') isa ColumnMajor + @test strides(V) == (5,1) + @test strides(V') == (1,5) + @test mul(V, b) ≈ V*b ≈ Matrix(V)*b + @test mul(V', b) ≈ V'b ≈ Matrix(V)'*b end + end - @testset "Float * Float Matrix" begin - A = randn(Float64, 100, 100) - x = randn(Float64, 100, 100) + @testset "Lmul/Rmul" begin + @testset "tri Lmul" begin + @testset "Float * Float vector" begin + A = randn(Float64, 100, 100) + x = randn(Float64, 100) - @test UpperTriangular(A)*x ≈ (similar(x) .= Lmul(UpperTriangular(A), x)) - end + L = Lmul(UpperTriangular(A),x) + @test size(L) == (size(L,1),) == (100,) + @test axes(L) == (axes(L,1),) == (Base.OneTo(100),) + @test eltype(L) == Float64 + @test length(L) == 100 - @testset "adjtrans" begin - for T in (Float64, ComplexF64) - A = randn(T,100,100) - b = randn(T,100) - - @test all(copy(Lmul(UpperTriangular(A)', b)) .=== UpperTriangular(A)'b) - @test all(copy(Lmul(UnitUpperTriangular(A)', b)) .=== UnitUpperTriangular(A)'b) - @test all(copy(Lmul(LowerTriangular(A)', b)) .=== LowerTriangular(A)'b) - @test all(copy(Lmul(UnitLowerTriangular(A)', b)) .=== UnitLowerTriangular(A)'b) - - @test all(copy(Lmul(transpose(UpperTriangular(A)), b)) .=== transpose(UpperTriangular(A))b) - @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), b)) .=== transpose(UnitUpperTriangular(A))b) - @test all(copy(Lmul(transpose(LowerTriangular(A)), b)) .=== transpose(LowerTriangular(A))b) - @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), b)) .=== transpose(UnitLowerTriangular(A))b) - - B = randn(T,100,100) - - @test all(copy(Lmul(UpperTriangular(A)', B)) .=== UpperTriangular(A)'B) - @test all(copy(Lmul(UnitUpperTriangular(A)', B)) .=== UnitUpperTriangular(A)'B) - @test all(copy(Lmul(LowerTriangular(A)', B)) .=== LowerTriangular(A)'B) - @test all(copy(Lmul(UnitLowerTriangular(A)', B)) .=== UnitLowerTriangular(A)'B) - - @test all(copy(Lmul(transpose(UpperTriangular(A)), B)) .=== transpose(UpperTriangular(A))B) - @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), B)) .=== transpose(UnitUpperTriangular(A))B) - @test all(copy(Lmul(transpose(LowerTriangular(A)), B)) .=== transpose(LowerTriangular(A))B) - @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), B)) .=== transpose(UnitLowerTriangular(A))B) + @test similar(L) isa Vector{Float64} + @test similar(L,Int) isa Vector{Int} + + @test all(copy(Lmul(UpperTriangular(A),x)) .=== + ArrayLayouts.lmul!(UpperTriangular(A),copy(x)) .=== + copyto!(similar(x),Lmul(UpperTriangular(A),x)) .=== + UpperTriangular(A)*x .=== + BLAS.trmv!('U', 'N', 'N', A, copy(x))) + @test all(copyto!(similar(x),Lmul(UnitUpperTriangular(A),x)) .=== + UnitUpperTriangular(A)*x .=== + BLAS.trmv!('U', 'N', 'U', A, copy(x))) + @test all(copyto!(similar(x),Lmul(LowerTriangular(A),x)) .=== + LowerTriangular(A)*x .=== + BLAS.trmv!('L', 'N', 'N', A, copy(x))) + @test all(ArrayLayouts.lmul!(UnitLowerTriangular(A),copy(x)) .=== + UnitLowerTriangular(A)x .=== + BLAS.trmv!('L', 'N', 'U', A, copy(x))) + + @test all(ArrayLayouts.lmul!(UpperTriangular(A)',copy(x)) .=== + ArrayLayouts.lmul!(LowerTriangular(A'),copy(x)) .=== + UpperTriangular(A)'x .=== + BLAS.trmv!('U', 'T', 'N', A, copy(x))) + @test all(ArrayLayouts.lmul!(UnitUpperTriangular(A)',copy(x)) .=== + ArrayLayouts.lmul!(UnitLowerTriangular(A'),copy(x)) .=== + UnitUpperTriangular(A)'*x .=== + BLAS.trmv!('U', 'T', 'U', A, copy(x))) + @test all(ArrayLayouts.lmul!(LowerTriangular(A)',copy(x)) .=== + ArrayLayouts.lmul!(UpperTriangular(A'),copy(x)) .=== + LowerTriangular(A)'*x .=== + BLAS.trmv!('L', 'T', 'N', A, copy(x))) + @test all(ArrayLayouts.lmul!(UnitLowerTriangular(A)',copy(x)) .=== + ArrayLayouts.lmul!(UnitUpperTriangular(A'),copy(x)) .=== + UnitLowerTriangular(A)'*x .=== + BLAS.trmv!('L', 'T', 'U', A, copy(x))) end - for T in (Float64, ComplexF64) - A = big.(randn(T,100,100)) - b = big.(randn(T,100)) + @testset "Float * Complex vector" begin + T = ComplexF64 + A = randn(T, 100, 100) + x = randn(T, 100) + + @test all((y = copy(x); y .= Lmul(UpperTriangular(A),y) ) .=== + UpperTriangular(A)x .=== + BLAS.trmv!('U', 'N', 'N', A, copy(x))) + @test all((y = copy(x); y .= Lmul(UnitUpperTriangular(A),y) ) .=== + UnitUpperTriangular(A)x .=== + BLAS.trmv!('U', 'N', 'U', A, copy(x))) + @test all((y = copy(x); y .= Lmul(LowerTriangular(A),y) ) .=== + LowerTriangular(A)x .=== + BLAS.trmv!('L', 'N', 'N', A, copy(x))) + @test all((y = copy(x); y .= Lmul(UnitLowerTriangular(A),y) ) .=== + UnitLowerTriangular(A)x .=== + BLAS.trmv!('L', 'N', 'U', A, copy(x))) + @test LowerTriangular(A') == UpperTriangular(A)' + + @test all((y = copy(x); y .= Lmul(transpose(UpperTriangular(A)),y) ) .=== + (similar(x) .= Lmul(LowerTriangular(transpose(A)),x)) .=== + transpose(UpperTriangular(A))x .=== + BLAS.trmv!('U', 'T', 'N', A, copy(x))) + @test all((y = copy(x); y .= Lmul(transpose(UnitUpperTriangular(A)),y) ) .=== + (similar(x) .= Lmul(UnitLowerTriangular(transpose(A)),x)) .=== + transpose(UnitUpperTriangular(A))x .=== + BLAS.trmv!('U', 'T', 'U', A, copy(x))) + @test all((y = copy(x); y .= Lmul(transpose(LowerTriangular(A)),y) ) .=== + (similar(x) .= Lmul(UpperTriangular(transpose(A)),x)) .=== + transpose(LowerTriangular(A))x .=== + BLAS.trmv!('L', 'T', 'N', A, copy(x))) + @test all((y = copy(x); y .= Lmul(transpose(UnitLowerTriangular(A)),y) ) .=== + (similar(x) .= Lmul(UnitUpperTriangular(transpose(A)),x)) .=== + transpose(UnitLowerTriangular(A))x .=== + BLAS.trmv!('L', 'T', 'U', A, copy(x))) + + @test all((y = copy(x); y .= Lmul(UpperTriangular(A)',y) ) .=== + (similar(x) .= Lmul(LowerTriangular(A'),x)) .=== + UpperTriangular(A)'x .=== + BLAS.trmv!('U', 'C', 'N', A, copy(x))) + @test all((y = copy(x); y .= Lmul(UnitUpperTriangular(A)',y) ) .=== + (similar(x) .= Lmul(UnitLowerTriangular(A'),x)) .=== + UnitUpperTriangular(A)'x .=== + BLAS.trmv!('U', 'C', 'U', A, copy(x))) + @test all((y = copy(x); y .= Lmul(LowerTriangular(A)',y) ) .=== + (similar(x) .= Lmul(UpperTriangular(A'),x)) .=== + LowerTriangular(A)'x .=== + BLAS.trmv!('L', 'C', 'N', A, copy(x))) + @test all((y = copy(x); y .= Lmul(UnitLowerTriangular(A)',y) ) .=== + (similar(x) .= Lmul(UnitUpperTriangular(A'),x)) .=== + UnitLowerTriangular(A)'x .=== + BLAS.trmv!('L', 'C', 'U', A, copy(x))) + end - @test all(copy(Lmul(UpperTriangular(A)', b)) ≈ UpperTriangular(A)'b) - @test all(copy(Lmul(UnitUpperTriangular(A)', b)) ≈ UnitUpperTriangular(A)'b) - @test all(copy(Lmul(LowerTriangular(A)', b)) ≈ LowerTriangular(A)'b) - @test all(copy(Lmul(UnitLowerTriangular(A)', b)) ≈ UnitLowerTriangular(A)'b) + @testset "Float * Float Matrix" begin + A = randn(Float64, 100, 100) + x = randn(Float64, 100, 100) - @test all(copy(Lmul(transpose(UpperTriangular(A)), b)) ≈ transpose(UpperTriangular(A))b) - @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), b)) ≈ transpose(UnitUpperTriangular(A))b) - @test all(copy(Lmul(transpose(LowerTriangular(A)), b)) ≈ transpose(LowerTriangular(A))b) - @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), b)) ≈ transpose(UnitLowerTriangular(A))b) + @test UpperTriangular(A)*x ≈ (similar(x) .= Lmul(UpperTriangular(A), x)) + end - B = big.(randn(T,100,100)) - - @test all(copy(Lmul(UpperTriangular(A)', B)) ≈ UpperTriangular(A)'B) - @test all(copy(Lmul(UnitUpperTriangular(A)', B)) ≈ UnitUpperTriangular(A)'B) - @test all(copy(Lmul(LowerTriangular(A)', B)) ≈ LowerTriangular(A)'B) - @test all(copy(Lmul(UnitLowerTriangular(A)', B)) ≈ UnitLowerTriangular(A)'B) - - @test all(copy(Lmul(transpose(UpperTriangular(A)), B)) ≈ transpose(UpperTriangular(A))B) - @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), B)) ≈ transpose(UnitUpperTriangular(A))B) - @test all(copy(Lmul(transpose(LowerTriangular(A)), B)) ≈ transpose(LowerTriangular(A))B) - @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), B)) ≈ transpose(UnitLowerTriangular(A))B) + @testset "adjtrans" begin + for T in (Float64, ComplexF64) + A = randn(T,100,100) + b = randn(T,100) + + @test all(copy(Lmul(UpperTriangular(A)', b)) .=== UpperTriangular(A)'b) + @test all(copy(Lmul(UnitUpperTriangular(A)', b)) .=== UnitUpperTriangular(A)'b) + @test all(copy(Lmul(LowerTriangular(A)', b)) .=== LowerTriangular(A)'b) + @test all(copy(Lmul(UnitLowerTriangular(A)', b)) .=== UnitLowerTriangular(A)'b) + + @test all(copy(Lmul(transpose(UpperTriangular(A)), b)) .=== transpose(UpperTriangular(A))b) + @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), b)) .=== transpose(UnitUpperTriangular(A))b) + @test all(copy(Lmul(transpose(LowerTriangular(A)), b)) .=== transpose(LowerTriangular(A))b) + @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), b)) .=== transpose(UnitLowerTriangular(A))b) + + B = randn(T,100,100) + + @test all(copy(Lmul(UpperTriangular(A)', B)) .=== UpperTriangular(A)'B) + @test all(copy(Lmul(UnitUpperTriangular(A)', B)) .=== UnitUpperTriangular(A)'B) + @test all(copy(Lmul(LowerTriangular(A)', B)) .=== LowerTriangular(A)'B) + @test all(copy(Lmul(UnitLowerTriangular(A)', B)) .=== UnitLowerTriangular(A)'B) + + @test all(copy(Lmul(transpose(UpperTriangular(A)), B)) .=== transpose(UpperTriangular(A))B) + @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), B)) .=== transpose(UnitUpperTriangular(A))B) + @test all(copy(Lmul(transpose(LowerTriangular(A)), B)) .=== transpose(LowerTriangular(A))B) + @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), B)) .=== transpose(UnitLowerTriangular(A))B) + end + + for T in (Float64, ComplexF64) + A = big.(randn(T,100,100)) + b = big.(randn(T,100)) + + @test all(copy(Lmul(UpperTriangular(A)', b)) ≈ UpperTriangular(A)'b) + @test all(copy(Lmul(UnitUpperTriangular(A)', b)) ≈ UnitUpperTriangular(A)'b) + @test all(copy(Lmul(LowerTriangular(A)', b)) ≈ LowerTriangular(A)'b) + @test all(copy(Lmul(UnitLowerTriangular(A)', b)) ≈ UnitLowerTriangular(A)'b) + + @test all(copy(Lmul(transpose(UpperTriangular(A)), b)) ≈ transpose(UpperTriangular(A))b) + @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), b)) ≈ transpose(UnitUpperTriangular(A))b) + @test all(copy(Lmul(transpose(LowerTriangular(A)), b)) ≈ transpose(LowerTriangular(A))b) + @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), b)) ≈ transpose(UnitLowerTriangular(A))b) + + B = big.(randn(T,100,100)) + + @test all(copy(Lmul(UpperTriangular(A)', B)) ≈ UpperTriangular(A)'B) + @test all(copy(Lmul(UnitUpperTriangular(A)', B)) ≈ UnitUpperTriangular(A)'B) + @test all(copy(Lmul(LowerTriangular(A)', B)) ≈ LowerTriangular(A)'B) + @test all(copy(Lmul(UnitLowerTriangular(A)', B)) ≈ UnitLowerTriangular(A)'B) + + @test all(copy(Lmul(transpose(UpperTriangular(A)), B)) ≈ transpose(UpperTriangular(A))B) + @test all(copy(Lmul(transpose(UnitUpperTriangular(A)), B)) ≈ transpose(UnitUpperTriangular(A))B) + @test all(copy(Lmul(transpose(LowerTriangular(A)), B)) ≈ transpose(LowerTriangular(A))B) + @test all(copy(Lmul(transpose(UnitLowerTriangular(A)), B)) ≈ transpose(UnitLowerTriangular(A))B) + end end end - end - @testset "tri Rmul" begin - for T in (Float64, ComplexF64) - A = randn(T, 100,100) - B = randn(T, 100,100) - R = Rmul(copy(A), UpperTriangular(B)) - @test size(R) == (size(R,1),size(R,2)) == (100,100) - @test axes(R) == (axes(R,1),axes(R,2)) == (Base.OneTo(100),Base.OneTo(100)) - @test eltype(R) == T - @test length(R) == 100^2 - - @test similar(R) isa Matrix{T} - @test similar(R,Int) isa Matrix{Int} - - R2 = deepcopy(R) - @test all(BLAS.trmm('R', 'U', 'N', 'N', one(T), B, A) .=== copyto!(similar(R2), R2) .=== materialize!(R)) - @test R.A ≠ A - @test all(BLAS.trmm('R', 'U', 'T', 'N', one(T), B, A) .=== copy(Rmul(A, transpose(UpperTriangular(B)))) .=== A*transpose(UpperTriangular(B))) - @test all(BLAS.trmm('R', 'U', 'N', 'U', one(T), B, A) .=== copy(Rmul(A, UnitUpperTriangular(B))).=== A*UnitUpperTriangular(B)) - @test all(BLAS.trmm('R', 'U', 'T', 'U', one(T), B, A) .=== copy(Rmul(A, transpose(UnitUpperTriangular(B)))) .=== A*transpose(UnitUpperTriangular(B))) - @test all(BLAS.trmm('R', 'L', 'N', 'N', one(T), B, A) .=== copy(Rmul(A, LowerTriangular(B))).=== A*LowerTriangular(B)) - @test all(BLAS.trmm('R', 'L', 'T', 'N', one(T), B, A) .=== copy(Rmul(A, transpose(LowerTriangular(B)))) .=== A*transpose(LowerTriangular(B))) - @test all(BLAS.trmm('R', 'L', 'N', 'U', one(T), B, A) .=== copy(Rmul(A, UnitLowerTriangular(B))).=== A*UnitLowerTriangular(B)) - @test all(BLAS.trmm('R', 'L', 'T', 'U', one(T), B, A) .=== copy(Rmul(A, transpose(UnitLowerTriangular(B)))) .=== A*transpose(UnitLowerTriangular(B))) - - if T <: Complex - @test all(BLAS.trmm('R', 'U', 'C', 'N', one(T), B, A) .=== copy(Rmul(A, UpperTriangular(B)')) .=== A*UpperTriangular(B)') - @test all(BLAS.trmm('R', 'U', 'C', 'U', one(T), B, A) .=== copy(Rmul(A, UnitUpperTriangular(B)')) .=== A*UnitUpperTriangular(B)') - @test all(BLAS.trmm('R', 'L', 'C', 'N', one(T), B, A) .=== copy(Rmul(A, LowerTriangular(B)')) .=== A*LowerTriangular(B)') - @test all(BLAS.trmm('R', 'L', 'C', 'U', one(T), B, A) .=== copy(Rmul(A, UnitLowerTriangular(B)')) .=== A*UnitLowerTriangular(B)') + @testset "tri Rmul" begin + for T in (Float64, ComplexF64) + A = randn(T, 100,100) + B = randn(T, 100,100) + R = Rmul(copy(A), UpperTriangular(B)) + @test size(R) == (size(R,1),size(R,2)) == (100,100) + @test axes(R) == (axes(R,1),axes(R,2)) == (Base.OneTo(100),Base.OneTo(100)) + @test eltype(R) == T + @test length(R) == 100^2 + + @test similar(R) isa Matrix{T} + @test similar(R,Int) isa Matrix{Int} + + R2 = deepcopy(R) + @test all(BLAS.trmm('R', 'U', 'N', 'N', one(T), B, A) .=== copyto!(similar(R2), R2) .=== materialize!(R)) + @test R.A ≠ A + @test all(BLAS.trmm('R', 'U', 'T', 'N', one(T), B, A) .=== copy(Rmul(A, transpose(UpperTriangular(B)))) .=== A*transpose(UpperTriangular(B))) + @test all(BLAS.trmm('R', 'U', 'N', 'U', one(T), B, A) .=== copy(Rmul(A, UnitUpperTriangular(B))).=== A*UnitUpperTriangular(B)) + @test all(BLAS.trmm('R', 'U', 'T', 'U', one(T), B, A) .=== copy(Rmul(A, transpose(UnitUpperTriangular(B)))) .=== A*transpose(UnitUpperTriangular(B))) + @test all(BLAS.trmm('R', 'L', 'N', 'N', one(T), B, A) .=== copy(Rmul(A, LowerTriangular(B))).=== A*LowerTriangular(B)) + @test all(BLAS.trmm('R', 'L', 'T', 'N', one(T), B, A) .=== copy(Rmul(A, transpose(LowerTriangular(B)))) .=== A*transpose(LowerTriangular(B))) + @test all(BLAS.trmm('R', 'L', 'N', 'U', one(T), B, A) .=== copy(Rmul(A, UnitLowerTriangular(B))).=== A*UnitLowerTriangular(B)) + @test all(BLAS.trmm('R', 'L', 'T', 'U', one(T), B, A) .=== copy(Rmul(A, transpose(UnitLowerTriangular(B)))) .=== A*transpose(UnitLowerTriangular(B))) + + if T <: Complex + @test all(BLAS.trmm('R', 'U', 'C', 'N', one(T), B, A) .=== copy(Rmul(A, UpperTriangular(B)')) .=== A*UpperTriangular(B)') + @test all(BLAS.trmm('R', 'U', 'C', 'U', one(T), B, A) .=== copy(Rmul(A, UnitUpperTriangular(B)')) .=== A*UnitUpperTriangular(B)') + @test all(BLAS.trmm('R', 'L', 'C', 'N', one(T), B, A) .=== copy(Rmul(A, LowerTriangular(B)')) .=== A*LowerTriangular(B)') + @test all(BLAS.trmm('R', 'L', 'C', 'U', one(T), B, A) .=== copy(Rmul(A, UnitLowerTriangular(B)')) .=== A*UnitLowerTriangular(B)') + end end + + T = Float64 + A = big.(randn(100,100)) + B = big.(randn(100,100)) + @test ArrayLayouts.rmul!(copy(A),UpperTriangular(B)) == A*UpperTriangular(B) end - T = Float64 - A = big.(randn(100,100)) - B = big.(randn(100,100)) - @test ArrayLayouts.rmul!(copy(A),UpperTriangular(B)) == A*UpperTriangular(B) + @testset "Diagonal and SymTridiagonal" begin + A = randn(5,5) + B = Diagonal(randn(5)) + @test MemoryLayout(B) == DiagonalLayout{DenseColumnMajor}() + + @test A*B == ArrayLayouts.rmul!(copy(A),B) + @test B*A == ArrayLayouts.lmul!(B,copy(A)) + end end - @testset "Diagonal and SymTridiagonal" begin + @testset "MulAdd" begin A = randn(5,5) - B = Diagonal(randn(5)) - @test MemoryLayout(B) == DiagonalLayout{DenseColumnMajor}() - - @test A*B == ArrayLayouts.rmul!(copy(A),B) - @test B*A == ArrayLayouts.lmul!(B,copy(A)) - end -end - -@testset "MulAdd" begin - A = randn(5,5) - B = randn(5,4) - C = randn(5,4) - b = randn(5) - c = randn(5) - - M = MulAdd(2.0,A,B,3.0,C) - @test size(M) == size(C) - @test size(M,1) == size(C,1) - @test size(M,2) == size(C,2) - @test_broken size(M,3) == size(C,3) - @test length(M) == length(C) - @test axes(M) == axes(C) - @test eltype(M) == Float64 - @test materialize(M) ≈ 2.0A*B + 3.0C - - @test_throws DimensionMismatch materialize(MulAdd(2.0,A,randn(3),1.0,B)) - @test_throws DimensionMismatch materialize(MulAdd([1,2],A,B,[1,2],C)) - @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(3,4))) - @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(5,5))) - - B = randn(5,5) - C = randn(5,5) - @test materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) == - muladd!(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(copy(C))) == 2.0Diagonal(A)*Diagonal(B) + 3.0*Diagonal(C) - @test_broken materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) isa Diagonal - - @test materialize(MulAdd(1.0, Eye(5), A, 3.0, C)) == muladd!(1.0, Eye(5), A, 3.0, copy(C)) == A + 3.0C - @test materialize(MulAdd(1.0, A, Eye(5), 3.0, C)) == muladd!(1.0, A, Eye(5), 3.0, copy(C)) == A + 3.0C - - @testset "Degenerate" begin - C = BigFloat.(randn(3,3)) - @test muladd!(1.0, BigFloat.(randn(3,0)), randn(0,3), 2.0, copy(C) ) == 2C - C = BigFloat.(randn(0,0)) - @test muladd!(1.0, BigFloat.(randn(0,3)), randn(3,0), 2.0, copy(C)) == C + B = randn(5,4) + C = randn(5,4) + b = randn(5) + c = randn(5) + + M = MulAdd(2.0,A,B,3.0,C) + @test size(M) == size(C) + @test size(M,1) == size(C,1) + @test size(M,2) == size(C,2) + @test_broken size(M,3) == size(C,3) + @test length(M) == length(C) + @test axes(M) == axes(C) + @test eltype(M) == Float64 + @test materialize(M) ≈ 2.0A*B + 3.0C + + @test_throws DimensionMismatch materialize(MulAdd(2.0,A,randn(3),1.0,B)) + @test_throws DimensionMismatch materialize(MulAdd([1,2],A,B,[1,2],C)) + @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(3,4))) + @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(5,5))) + + B = randn(5,5) + C = randn(5,5) + @test materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) == + muladd!(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(copy(C))) == 2.0Diagonal(A)*Diagonal(B) + 3.0*Diagonal(C) + @test_broken materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) isa Diagonal + + @test materialize(MulAdd(1.0, Eye(5), A, 3.0, C)) == muladd!(1.0, Eye(5), A, 3.0, copy(C)) == A + 3.0C + @test materialize(MulAdd(1.0, A, Eye(5), 3.0, C)) == muladd!(1.0, A, Eye(5), 3.0, copy(C)) == A + 3.0C + + @testset "Degenerate" begin + C = BigFloat.(randn(3,3)) + @test muladd!(1.0, BigFloat.(randn(3,0)), randn(0,3), 2.0, copy(C) ) == 2C + C = BigFloat.(randn(0,0)) + @test muladd!(1.0, BigFloat.(randn(0,3)), randn(3,0), 2.0, copy(C)) == C + end + + @testset "Generic muladd" begin + A = view(reshape(1:6,2,3),1:2,:) + @test IndexStyle(typeof(A)) == IndexCartesian() + b = randn(3) + c = similar(b,2) + @test muladd!(1.0, A, b, 0.0, c) == A*b + end end end \ No newline at end of file