Internal type promotion fixes (#112)

Author: dkarrasch

src/LinearMaps.jl

@@ -48,15 +48,20 @@ Base.ndims(::LinearMap) = 2
 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
 Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
 
+# check dimension consistency for multiplication A*B
+_iscompatible((A, B)) = size(A, 2) == size(B, 1)
+function check_dim_mul(A, B)
+    _iscompatible((A, B)) ||
+        throw(DimensionMismatch("second dimension of left factor, $(size(A, 2)), " *
+            "does not match first dimension of right factor, $(size(B, 1))"))
+    return nothing
+end
 # check dimension consistency for multiplication C = A*B
 function check_dim_mul(C, A, B)
-    # @info "checked vector dimensions" # uncomment for testing
     mA, nA = size(A) # A always has two dimensions
     mB, nB = size(B, 1), size(B, 2)
-    (mB == nA) ||
-        throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)"))
-    (size(C, 1) != mA || size(C, 2) != nB) &&
-        throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)"))
+    (mB == nA && size(C, 1) == mA && size(C, 2) == nB) ||
+        throw(DimensionMismatch("A has size ($mA,$nA), B has size ($mB,$nB), C has size $(size(C))"))
     return nothing
 end
 
@@ -93,10 +98,8 @@ julia> A*x
 ```
 """
 function Base.:(*)(A::LinearMap, x::AbstractVector)
-    m, n = size(A)
-    n == length(x) || throw(DimensionMismatch("linear map has dimensions ($m,$n), " *
-        "vector has length $(length(x))"))
-    return mul!(similar(x, promote_type(eltype(A), eltype(x)), m), A, x)
+    check_dim_mul(A, x)
+    return mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
 end
 
 """

src/blockmap.jl

@@ -4,8 +4,9 @@ struct BlockMap{T,As<:Tuple{Vararg{LinearMap}},Rs<:Tuple{Vararg{Int}},Rranges<:T
     rowranges::Rranges
     colranges::Cranges
     function BlockMap{T,R,S}(maps::R, rows::S) where {T, R<:Tuple{Vararg{LinearMap}}, S<:Tuple{Vararg{Int}}}
-        for A in maps
-            promote_type(T, eltype(A)) == T || throw(InexactError())
+        for n in eachindex(maps)
+            A = maps[n]
+            @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in BlockMap constructor"
         end
         rowranges, colranges = rowcolranges(maps, rows)
         return new{T,R,S,typeof(rowranges),typeof(colranges)}(maps, rows, rowranges, colranges)
@@ -391,8 +392,9 @@ struct BlockDiagonalMap{T,As<:Tuple{Vararg{LinearMap}},Ranges<:Tuple{Vararg{Unit
     rowranges::Ranges
     colranges::Ranges
     function BlockDiagonalMap{T,As}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}}
-        for A in maps
-            promote_type(T, eltype(A)) == T || throw(InexactError())
+        for n in eachindex(maps)
+            A = maps[n]
+            @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in BlockDiagonalMap constructor"
         end
         # row ranges
         inds = vcat(1, size.(maps, 1)...)

src/composition.jl

@@ -1,15 +1,13 @@
-# helper function
-check_dim_mul(A, B) = size(A, 2) == size(B, 1)
-
 struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
     maps::As # stored in order of application to vector
     function CompositeMap{T, As}(maps::As) where {T, As}
         N = length(maps)
         for n in 2:N
-            check_dim_mul(maps[n], maps[n-1]) || throw(DimensionMismatch("CompositeMap"))
+            check_dim_mul(maps[n], maps[n-1])
         end
-        for n in 1:N
-            promote_type(T, eltype(maps[n])) == T || throw(InexactError())
+        for n in eachindex(maps)
+            A = maps[n]
+            @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in CompositeMap constructor"
         end
         new{T, As}(maps)
     end

src/kronecker.jl

@@ -1,8 +1,9 @@
 struct KroneckerMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
     maps::As
     function KroneckerMap{T, As}(maps::As) where {T, As}
-        for A in maps
-            promote_type(T, eltype(A)) == T || throw(InexactError())
+        for n in eachindex(maps)
+            A = maps[n]
+            @assert promote_type(T, eltype(A)) == T  "eltype $(eltype(A)) cannot be promoted to $T in KroneckerMap constructor"
         end
         return new{T,As}(maps)
     end
@@ -151,7 +152,7 @@ end
 function _unsafe_mul!(y::AbstractVecOrMat, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
     require_one_based_indexing(y)
     B, A = L.maps
-    if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps))
+    if length(A.maps) == length(B.maps) && all(_iscompatible, zip(A.maps, B.maps))
         _unsafe_mul!(y, kron(map(*, A.maps, B.maps)...), x)
     else
         _unsafe_mul!(y, LinearMap(A)*B, x)
@@ -166,7 +167,7 @@ function _unsafe_mul!(y::AbstractVecOrMat, L::CompositeMap{T,<:Tuple{Vararg{Kron
     Bs = map(AB -> AB.maps[2], L.maps)
     As1, As2 = Base.front(As), Base.tail(As)
     Bs1, Bs2 = Base.front(Bs), Base.tail(Bs)
-    apply = all(A -> check_dim_mul(A...), zip(As1, As2)) && all(A -> check_dim_mul(A...), zip(Bs1, Bs2))
+    apply = all(_iscompatible, zip(As1, As2)) && all(_iscompatible, zip(Bs1, Bs2))
     if apply
         _unsafe_mul!(y, kron(prod(As), prod(Bs)), x)
     else
@@ -181,9 +182,10 @@ end
 struct KroneckerSumMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T}
     maps::As
     function KroneckerSumMap{T, As}(maps::As) where {T, As}
-        for A in maps
+        for n in eachindex(maps)
+            A = maps[n]
             size(A, 1) == size(A, 2) || throw(ArgumentError("operators need to be square in Kronecker sums"))
-            promote_type(T, eltype(A)) == T || throw(InexactError())
+            @assert promote_type(T, eltype(A)) == T  "eltype $(eltype(A)) cannot be promoted to $T in KroneckerSumMap constructor"
         end
         return new{T,As}(maps)
     end

src/linearcombination.jl

@@ -3,9 +3,10 @@ struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
     function LinearCombination{T, As}(maps::As) where {T, As}
         N = length(maps)
         sz = size(maps[1])
-        for n in 1:N
-            size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination"))
-            promote_type(T, eltype(maps[n])) == T || throw(InexactError())
+        for n in eachindex(maps)
+            A = maps[n]
+            size(A) == sz || throw(DimensionMismatch("LinearCombination"))
+            @assert promote_type(T, eltype(A)) == T  "eltype $(eltype(A)) cannot be promoted to $T in LinearCombination constructor"
         end
         new{T, As}(maps)
     end

src/scaledmap.jl

@@ -1,15 +1,15 @@
 struct ScaledMap{T, S<:RealOrComplex, A<:LinearMap} <: LinearMap{T}
     λ::S
     lmap::A
-    function ScaledMap{T,S,A}(λ::S, lmap::A) where {T, S <: RealOrComplex, A <: LinearMap}
-        Base.promote_op(*, S, eltype(lmap)) == T || throw(InexactError())
+    function ScaledMap{T}(λ::S, lmap::A) where {T, S <: RealOrComplex, A <: LinearMap}
+        @assert Base.promote_op(*, S, eltype(lmap)) == T "target type $T cannot hold products of $S and $(eltype(lmap)) objects"
         new{T,S,A}(λ, lmap)
     end
 end
 
 # constructor
-ScaledMap{T}(λ::S, lmap::A) where {T,S<:RealOrComplex,A<:LinearMap} =
-    ScaledMap{Base.promote_op(*, S, eltype(lmap)),S,A}(λ, lmap)
+ScaledMap(λ::S, lmap::A) where {S<:RealOrComplex,A<:LinearMap} =
+    ScaledMap{Base.promote_op(*, S, eltype(lmap))}(λ, lmap)
 
 # show
 function Base.show(io::IO, A::ScaledMap{T}) where {T}
@@ -32,14 +32,8 @@ Base.:(==)(A::ScaledMap, B::ScaledMap) =
     (eltype(A) == eltype(B) && A.lmap == B.lmap) && A.λ == B.λ
 
 # scalar multiplication and division
-function Base.:(*)(α::RealOrComplex, A::LinearMap)
-    T = Base.promote_op(*, typeof(α), eltype(A))
-    return ScaledMap{T}(α, A)
-end
-function Base.:(*)(A::LinearMap, α::RealOrComplex)
-    T = Base.promote_op(*, typeof(α), eltype(A))
-    return ScaledMap{T}(α, A)
-end
+Base.:(*)(α::RealOrComplex, A::LinearMap) = ScaledMap(α, A)
+Base.:(*)(A::LinearMap, α::RealOrComplex) = ScaledMap(α, A)
 
 Base.:(*)(α::Number, A::ScaledMap) = (α * A.λ) * A.lmap
 Base.:(*)(A::ScaledMap, α::Number) = A.lmap * (A.λ * α)

test/blockmap.jl

@@ -9,6 +9,9 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             L = @inferred hcat(LinearMap(A11), LinearMap(A12))
             @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
             @test L isa LinearMaps.BlockMap{elty}
+            if elty <: Complex
+                @test_throws AssertionError LinearMaps.BlockMap{Float64}((LinearMap(A11), LinearMap(A12)), (2,))
+            end
             A = [A11 A12]
             x = rand(10+n2)
             @test size(L) == size(A)
@@ -182,6 +185,9 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             M3 = randn(elty, m, n+3); L3 = LinearMap(M3)
 
             # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse:
+            if elty <: Complex
+                @test_throws AssertionError LinearMaps.BlockDiagonalMap{Float64}((L1, L2, L3, L2, L1))
+            end
             Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...))
             @test (@which blockdiag(sparse.((M1, M2, M3, M2, M1))...)).module != LinearMaps
             @test (@which cat(M1, M2, M3, M2, M1; dims=(1,2))).module != LinearMaps

test/composition.jl

@@ -6,6 +6,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     FCM = LinearMaps.CompositeMap{ComplexF64}((FC,))
     L = LowerTriangular(ones(10,10))
     @test_throws DimensionMismatch F * LinearMap(rand(2,2))
+    @test_throws AssertionError LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10))))
     A = 2 * rand(ComplexF64, (10, 10)) .- 1
     B = rand(size(A)...)
     H = LinearMap(Hermitian(A'A))

test/kronecker.jl

@@ -8,6 +8,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
         LA = LinearMap(A)
         LB = LinearMap(B)
         LK = @inferred kron(LA, LB)
+        @test_throws AssertionError LinearMaps.KroneckerMap{Float64}((LA, LB))
         @test @inferred size(LK) == size(K)
         @test LinearMaps.MulStyle(LK) === LinearMaps.ThreeArg()
         for i in (1, 2)

test/linearcombination.jl

@@ -10,6 +10,7 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     @test run(b, samples=3).allocs == 0
     n = 10
     L = sum(fill(CS!, n))
+    @test_throws AssertionError LinearMaps.LinearCombination{Float64}((CS!, CS!))
     @test mul!(u, L, v) ≈ n * cumsum(v)
     b = @benchmarkable mul!($u, $L, $v, 2, 2)
     @test run(b, samples=5).allocs <= 1