Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ Quaternions = "0.7"
Random = "1.6"
SparseArrays = "1.6"
StableRNGs = "1"
StaticArrays = "1"
Test = "1.6"
julia = "1.6"

Expand All @@ -33,7 +34,8 @@ Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Aqua", "Infinities", "Quaternions", "Random", "StableRNGs", "SparseArrays", "Test"]
test = ["Aqua", "Infinities", "Quaternions", "Random", "StableRNGs", "SparseArrays", "StaticArrays", "Test"]
10 changes: 8 additions & 2 deletions src/ArrayLayouts.jl
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,12 @@ else
LinearAlgebra.UnitLowerTriangular{T,S}}
end

@static if VERSION ≥ v"1.8.0"
import Base: LazyString
else
const LazyString = string
end

# Originally defined in FillArrays
_copy_oftype(A::AbstractArray, ::Type{S}) where {S} = eltype(A) == S ? copy(A) : AbstractArray{S}(A)
_copy_oftype(A::AbstractRange, ::Type{S}) where {S} = eltype(A) == S ? copy(A) : map(S, A)
Expand Down Expand Up @@ -205,7 +211,7 @@ getindex(A::AdjOrTrans{<:Any,<:LayoutVector}, kr::Integer, jr::AbstractVector) =
function *(a::Transpose{T, <:LayoutVector{T}}, b::Zeros{T, 1}) where T<:Real
la, lb = length(a), length(b)
if la ≠ lb
throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
throw(DimensionMismatch(LazyString("dot product arguments have lengths ", la, " and ", lb)))
end
return zero(T)
end
Expand Down Expand Up @@ -334,7 +340,7 @@ end
@inline function reflectorApply!(x::AbstractVector, τ::Number, A::AbstractVecOrMat)
m,n = size(A,1),size(A,2)
if length(x) != m
throw(DimensionMismatch("reflector has length $(length(x)), which must match the first dimension of matrix A, $m"))
throw(DimensionMismatch(LazyString("reflector has length ", length(x), ", which must match the first dimension of matrix A, ", m)))
end
m == 0 && return A
@inbounds begin
Expand Down
16 changes: 8 additions & 8 deletions src/factorizations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -152,8 +152,8 @@ function copyto!(dest::AbstractArray, M::Ldiv{<:AbstractQLayout})
ldiv!(A,dest)
end

materialize!(M::Lmul{LAY}) where LAY<:AbstractQLayout = error("Overload materialize!(::Lmul{$(LAY)})")
materialize!(M::Rmul{LAY}) where LAY<:AbstractQLayout = error("Overload materialize!(::Rmul{$(LAY)})")
materialize!(M::Lmul{LAY}) where LAY<:AbstractQLayout = error(LazyString("Overload materialize!(::Lmul{", LAY, "})"))
materialize!(M::Rmul{LAY}) where LAY<:AbstractQLayout = error(LazyString("Overload materialize!(::Rmul{", LAY, "})"))

materialize!(M::Ldiv{<:AbstractQLayout}) = materialize!(Lmul(M.A',M.B))

Expand All @@ -179,7 +179,7 @@ function materialize!(M::Lmul{<:QRPackedQLayout})
mA, nA = size(A.factors)
mB, nB = size(B,1), size(B,2)
if mA != mB
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
throw(DimensionMismatch(LazyString("matrix A has dimensions (", mA, ",", nA, ") but B has dimensions (", mB, ", ", nB, ")")))
end
Afactors = A.factors
@inbounds begin
Expand Down Expand Up @@ -217,7 +217,7 @@ function materialize!(M::Lmul{<:AdjQRPackedQLayout})
mA, nA = size(A.factors)
mB, nB = size(B,1), size(B,2)
if mA != mB
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
throw(DimensionMismatch(LazyString("matrix A has dimensions (", mA, ",", nA, ") but B has dimensions (", mB, ", ", nB, ")")))
end
Afactors = A.factors
@inbounds begin
Expand Down Expand Up @@ -248,7 +248,7 @@ function materialize!(M::Rmul{<:Any,<:QRPackedQLayout})
mQ, nQ = size(Q.factors)
mA, nA = size(A,1), size(A,2)
if nA != mQ
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
throw(DimensionMismatch(LazyString("matrix A has dimensions (", mA, ",", nA, ") but matrix Q has dimensions (", mQ, ", ", nQ, ")")))
end
Qfactors = Q.factors
@inbounds begin
Expand Down Expand Up @@ -284,7 +284,7 @@ function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
mQ, nQ = size(Q.factors)
mA, nA = size(A,1), size(A,2)
if nA != mQ
throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
throw(DimensionMismatch(LazyString("matrix A has dimensions (", mA, ",", nA, ") but matrix Q has dimensions (", mQ, ", ", nQ, ")")))
end
Qfactors = Q.factors
@inbounds begin
Expand All @@ -309,10 +309,10 @@ end
__qr(layout, lengths, A; kwds...) = invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_qr(layout, axes, A; kwds...) = __qr(layout, map(length, axes), A; kwds...)
_qr(layout, axes, A, pivot::P; kwds...) where P = invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
_qr!(layout, axes, A, args...; kwds...) = error(LazyString("Overload _qr!(::", typeof(layout), ", axes, A)"))
_lu(layout, axes, A; kwds...) = invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_lu(layout, axes, A, pivot::P; kwds...) where P = invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
_lu!(layout, axes, A, args...; kwds...) = error(LazyString("Overload _lu!(::", typeof(layout), ", axes, A)"))
_cholesky(layout, axes, A, ::CNoPivot=CNoPivot(); check::Bool = true) = cholesky!(cholcopy(A); check = check)
_cholesky(layout, axes, A, ::CRowMaximum; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), CRowMaximum(); tol = tol, check = check)
_factorize(layout, axes, A) = qr(A) # Default to QR
Expand Down
10 changes: 5 additions & 5 deletions src/ldiv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,10 +56,10 @@ _getindex(::Type{Tuple{I,J}}, L::Ldiv, (k,j)::Tuple{Colon,J}) where {I,J} = Ldiv
_getindex(::Type{Tuple{I,J}}, L::Ldiv, (k,j)::Tuple{I,J}) where {I,J} = L[:,j][k]

check_ldiv_axes(A, B) =
axes(A,1) == axes(B,1) || throw(DimensionMismatch("First axis of A, $(axes(A,1)), and first axis of B, $(axes(B,1)) must match"))
axes(A,1) == axes(B,1) || throw(DimensionMismatch(LazyString("First axis of A, ", axes(A,1), ", and first axis of B, ", axes(B,1), " must match")))

check_rdiv_axes(A, B) =
axes(A,2) == axes(B,2) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and second axis of B, $(axes(B,2)) must match"))
axes(A,2) == axes(B,2) || throw(DimensionMismatch(LazyString("Second axis of A, ", axes(A,2), ", and second axis of B, ", axes(B,2), " must match")))



Expand All @@ -73,12 +73,12 @@ end
Rdiv(instantiate(L.A), instantiate(L.B))
end

__ldiv!(::Mat, ::Mat, B) where Mat = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(Mat))),$(typeof(MemoryLayout(B)))})")
__ldiv!(::Mat, ::Mat, B::LayoutArray) where Mat = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(Mat))),$(typeof(MemoryLayout(B)))})")
__ldiv!(::Mat, ::Mat, B) where Mat = error(LazyString("Overload materialize!(::Ldiv{", typeof(MemoryLayout(Mat)), ",", typeof(MemoryLayout(B)), "})"))
__ldiv!(::Mat, ::Mat, B::LayoutArray) where Mat = error(LazyString("Overload materialize!(::Ldiv{", typeof(MemoryLayout(Mat)), ",", typeof(MemoryLayout(B)), "})"))
__ldiv!(_, F, B) = LinearAlgebra.ldiv!(F, B)
@inline _ldiv!(A, B) = __ldiv!(A, factorize(A), B)
@inline _ldiv!(A::Factorization, B) = LinearAlgebra.ldiv!(A, B)
@inline _ldiv!(A::Factorization, B::LayoutArray) = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(A))),$(typeof(MemoryLayout(B)))})")
@inline _ldiv!(A::Factorization, B::LayoutArray) = error(LazyString("Overload materialize!(::Ldiv{", typeof(MemoryLayout(A)), ",", typeof(MemoryLayout(B)), "})"))

@inline _ldiv!(dest, A, B; kwds...) = ldiv!(dest, factorize(A), B; kwds...)
@inline _ldiv!(dest, A::Factorization, B; kwds...) = LinearAlgebra.ldiv!(dest, A, B; kwds...)
Expand Down
4 changes: 2 additions & 2 deletions src/memorylayout.jl
Original file line number Diff line number Diff line change
Expand Up @@ -585,8 +585,8 @@ 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::Bidiagonal) = D.uplo == 'U' ? D.ev : throw(ArgumentError(LazyString(D, " is lower-bidiagonal")))
subdiagonaldata(D::Bidiagonal) = D.uplo == 'L' ? D.ev : throw(ArgumentError(LazyString(D, " is upper-bidiagonal")))

supdiagonaldata(D::SymTridiagonal) = D.ev
subdiagonaldata(D::SymTridiagonal) = D.ev
Expand Down
18 changes: 14 additions & 4 deletions src/mul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,18 +92,28 @@ check_mul_axes(A) = nothing
_check_mul_axes(::Number, ::Number) = nothing
_check_mul_axes(::Number, _) = nothing
_check_mul_axes(_, ::Number) = nothing
_check_mul_axes(A, B) = axes(A, 2) == axes(B, 1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))
_check_mul_axes(A, B) = axes(A, 2) == axes(B, 1) || throw_mul_axes_err(axes(A,2), axes(B,1))
@noinline function throw_mul_axes_err(axA2, axB1)
throw(
DimensionMismatch(
LazyString("second axis of A, ", axA2, ", and first axis of B, ", axB1, ", must match")))
end
@noinline function throw_mul_axes_err(axA2::Base.OneTo, axB1::Base.OneTo)
throw(
DimensionMismatch(
LazyString("second dimension of A, ", length(axA2), ", does not match length of x, ", length(axB1))))
end
# we need to special case AbstractQ as it allows non-compatiple multiplication
const FlexibleLeftQs = Union{QRCompactWYQ,QRPackedQ,HessenbergQ}
_check_mul_axes(::FlexibleLeftQs, ::Number) = nothing
_check_mul_axes(Q::FlexibleLeftQs, B) =
axes(Q.factors, 1) == axes(B, 1) || axes(Q.factors, 2) == axes(B, 1) ||
throw(DimensionMismatch("First axis of B, $(axes(B,1)) must match either axes of A, $(axes(Q.factors))"))
throw(DimensionMismatch(LazyString("First axis of B, ", axes(B,1), " must match either axes of A, ", axes(Q.factors))))
_check_mul_axes(::Number, ::AdjointQtype{<:Any,<:FlexibleLeftQs}) = nothing
function _check_mul_axes(A, adjQ::AdjointQtype{<:Any,<:FlexibleLeftQs})
Q = parent(adjQ)
axes(A, 2) == axes(Q.factors, 1) || axes(A, 2) == axes(Q.factors, 2) ||
throw(DimensionMismatch("Second axis of A, $(axes(A,2)) must match either axes of B, $(axes(Q.factors))"))
throw(DimensionMismatch(LazyString("Second axis of A, ", axes(A,2), " must match either axes of B, ", axes(Q.factors))))
end
_check_mul_axes(Q::FlexibleLeftQs, adjQ::AdjointQtype{<:Any,<:FlexibleLeftQs}) =
invoke(_check_mul_axes, Tuple{Any,Any}, Q, adjQ)
Expand All @@ -115,7 +125,7 @@ end
# we need to special case AbstractQ as it allows non-compatiple multiplication
function check_mul_axes(A::Union{QRCompactWYQ,QRPackedQ}, B, C...)
axes(A.factors, 1) == axes(B, 1) || axes(A.factors, 2) == axes(B, 1) ||
throw(DimensionMismatch("First axis of B, $(axes(B,1)) must match either axes of A, $(axes(A))"))
throw(DimensionMismatch(LazyString("First axis of B, ", axes(B,1), " must match either axes of A, ", axes(A))))
check_mul_axes(B, C...)
end

Expand Down
4 changes: 2 additions & 2 deletions src/muladd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,8 @@ similar(M::MulAdd) = similar(M, eltype(M))
function checkdimensions(M::MulAdd)
@boundscheck check_mul_axes(M.α, M.A, M.B)
@boundscheck check_mul_axes(M.β, M.C)
@boundscheck axes(M.A,1) == axes(M.C,1) || throw(DimensionMismatch("First axis of A, $(axes(M.A,1)), and first axis of C, $(axes(M.C,1)) must match"))
@boundscheck axes(M.B,2) == axes(M.C,2) || throw(DimensionMismatch("Second axis of B, $(axes(M.B,2)), and second axis of C, $(axes(M.C,2)) must match"))
@boundscheck axes(M.A,1) == axes(M.C,1) || throw(DimensionMismatch(LazyString("First axis of A, ", axes(M.A,1), ", and first axis of C, ", axes(M.C,1), " must match")))
@boundscheck axes(M.B,2) == axes(M.C,2) || throw(DimensionMismatch(LazyString("Second axis of B, ", axes(M.B,2), ", and second axis of C, ", axes(M.C,2), " must match")))
end
@propagate_inbounds function instantiate(M::MulAdd)
checkdimensions(M)
Expand Down
32 changes: 11 additions & 21 deletions src/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ function materialize!(M::Lmul{<:TriangularLayout{'U','N'}})
A,B = M.A,M.B
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"))
throw(DimensionMismatch(LazyString("right hand side B needs first dimension of size ", size(A,1), ", has size ", m)))
end
Adata = triangulardata(A)
for j = rowsupport(B)
Expand All @@ -46,7 +46,7 @@ function materialize!(M::Lmul{<:TriangularLayout{'U','U'}})
A,B = M.A,M.B
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"))
throw(DimensionMismatch(LazyString("right hand side B needs first dimension of size ", size(A,1), ", has size ", m)))
end
Adata = triangulardata(A)
for j = rowsupport(B)
Expand All @@ -66,7 +66,7 @@ function materialize!(M::Lmul{<:TriangularLayout{'L','N'}})
A,B = M.A,M.B
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"))
throw(DimensionMismatch(LazyString("right hand side B needs first dimension of size ", size(A,1), ", has size ", m)))
end
Adata = triangulardata(A)
for j = 1:n
Expand All @@ -84,7 +84,7 @@ function materialize!(M::Lmul{<:TriangularLayout{'L','U'}})
A,B = M.A,M.B
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"))
throw(DimensionMismatch(LazyString("right hand side B needs first dimension of size ", size(A,1), ", has size ", m)))
end
Adata = triangulardata(A)
for j = 1:n
Expand Down Expand Up @@ -264,10 +264,7 @@ end
function materialize!(M::MatLdivVec{<:TriangularLayout{'U','N'}})
A,b = M.A,M.B
require_one_based_indexing(A, b)
n = size(A, 2)
if !(n == length(b))
throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
end
check_mul_axes(A, b)
data = triangulardata(A)
@inbounds for j in reverse(colsupport(b,1))
iszero(data[j,j]) && throw(SingularException(j))
Expand All @@ -282,10 +279,7 @@ end
function materialize!(M::MatLdivVec{<:TriangularLayout{'U','U'}})
A,b = M.A,M.B
require_one_based_indexing(A, b)
n = size(A, 2)
if !(n == length(b))
throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
end
check_mul_axes(A, b)
data = triangulardata(A)
@inbounds for j in reverse(colsupport(b,1))
iszero(data[j,j]) && throw(SingularException(j))
Expand All @@ -300,11 +294,9 @@ end
function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N'}})
A,b = M.A,M.B
require_one_based_indexing(A, b)
n = size(A, 2)
if !(n == length(b))
throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
end
check_mul_axes(A, b)
data = triangulardata(A)
n = size(A, 2)
@inbounds for j in 1:n
iszero(data[j,j]) && throw(SingularException(j))
bj = b[j] = data[j,j] \ b[j]
Expand All @@ -318,11 +310,9 @@ end
function materialize!(M::MatLdivVec{<:TriangularLayout{'L','U'}})
A,b = M.A,M.B
require_one_based_indexing(A, b)
n = size(A, 2)
if !(n == length(b))
throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
end
check_mul_axes(A, b)
data = triangulardata(A)
n = size(A, 2)
@inbounds for j in 1:n
iszero(data[j,j]) && throw(SingularException(j))
bj = b[j]
Expand Down Expand Up @@ -379,7 +369,7 @@ function materialize!(M::MatLdivVec{<:BidiagonalLayout})
require_one_based_indexing(A, b)
N = size(A, 2)
if N != length(b)
throw(DimensionMismatch("second dimension of A, $N, does not match one of the length of b, $(length(b))"))
throw(DimensionMismatch(LazyString("second dimension of A, ", N, ", does not match one of the length of b, ", length(b))))
end

if N == 0
Expand Down
20 changes: 20 additions & 0 deletions test/test_ldiv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ module TestLdiv

using ArrayLayouts, LinearAlgebra, FillArrays, Test
import ArrayLayouts: ApplyBroadcastStyle, QRCompactWYQLayout, QRCompactWYLayout, QRPackedQLayout, QRPackedLayout
using StaticArrays

@testset "Ldiv" begin
@testset "Float64 \\ *" begin
Expand Down Expand Up @@ -295,6 +296,25 @@ import ArrayLayouts: ApplyBroadcastStyle, QRCompactWYQLayout, QRCompactWYLayout,
@test ArrayLayouts.rdiv!(similar(B), B, D) == B / D
@test_broken ArrayLayouts.rdiv!(similar(B), D, B) == D / B
end

@testset "error paths" begin
v = rand(Int,3)
A = rand(2,2)
S = SMatrix{2,2}(A)
errA(U) = DimensionMismatch("second dimension of A, $(size(U,2)), does not match length of x, $(length(v))")
errS(U) = DimensionMismatch("second axis of A, $(axes(U,2)), and first axis of B, $(axes(v,1)), must match")
for (M, errf) in ((A, errA), (S, errS))
U = UpperTriangular(M)
err = errf(U)
@test_throws err ArrayLayouts.materialize!(ArrayLayouts.Ldiv(U, v))
UU = UnitUpperTriangular(M)
@test_throws err ArrayLayouts.materialize!(ArrayLayouts.Ldiv(UU, v))
L = LowerTriangular(M)
@test_throws err ArrayLayouts.materialize!(ArrayLayouts.Ldiv(L, v))
UL = UnitLowerTriangular(M)
@test_throws err ArrayLayouts.materialize!(ArrayLayouts.Ldiv(UL, v))
end
end
end

end
13 changes: 13 additions & 0 deletions test/test_muladd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -845,6 +845,19 @@ Random.seed!(0)
end
end

@testset "Error paths" begin
if VERSION >= v"1.10.0"
Q = qr(rand(2,2), ColumnNorm()).Q
else
Q = qr(rand(2,2), Val(true)).Q
end
v = rand(Float32, 3)
@test_throws DimensionMismatch ArrayLayouts.materialize!(ArrayLayouts.Rmul(v, Q))
@test_throws DimensionMismatch ArrayLayouts.materialize!(ArrayLayouts.Rmul(v, Q'))
@test_throws DimensionMismatch ArrayLayouts.materialize!(ArrayLayouts.Lmul(Q, v))
@test_throws DimensionMismatch ArrayLayouts.materialize!(ArrayLayouts.Lmul(Q', v))
end

@testset "dual" begin
a = randn(5)
X = randn(5,6)
Expand Down