--- a/src/LazyTensors/lazy_tensor_operations.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -1,5 +1,15 @@
 using Sbplib.RegionIndices
 
+export LazyTensorMappingApplication
+export LazyTensorMappingTranspose
+export TensorMappingComposition
+export LazyLinearMap
+export IdentityMapping
+export InflatedTensorMapping
+export LazyOuterProduct
+export ⊗
+export SizeMismatch
+
 """
     LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R}
 
@@ -9,12 +19,16 @@
 With a mapping `m` and a vector `v` the LazyTensorMappingApplication object can be created by `m*v`.
 The actual result will be calcualted when indexing into `m*v`.
 """
-struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{T,R,D}, AA<:AbstractArray{T,D}} <: LazyArray{T,R}
+struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{<:Any,R,D}, AA<:AbstractArray{<:Any,D}} <: LazyArray{T,R}
     t::TM
     o::AA
+
+    function LazyTensorMappingApplication(t::TensorMapping{<:Any,R,D}, o::AbstractArray{<:Any,D}) where {R,D}
+        T = promote_type(eltype(t), eltype(o))
+        return new{T,R,D,typeof(t), typeof(o)}(t,o)
+    end
 end
 # TODO: Do boundschecking on creation!
-export LazyTensorMappingApplication
 
 Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Any,R}) where {T,R} = apply(ta.t, ta.o, I...)
 Base.getindex(ta::LazyTensorMappingApplication{T,1}, I::CartesianIndex{1}) where {T} = apply(ta.t, ta.o, I.I...) # Would otherwise be caught in the previous method.
@@ -43,15 +57,14 @@
 struct LazyTensorMappingTranspose{T,R,D, TM<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R}
     tm::TM
 end
-export LazyTensorMappingTranspose
 
 # # TBD: Should this be implemented on a type by type basis or through a trait to provide earlier errors?
 # Jonatan 2020-09-25: Is the problem that you can take the transpose of any TensorMapping even if it doesn't implement `apply_transpose`?
 Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm)
 Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm
 
-apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...)
-apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...)
+apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...)
+apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...)
 
 range_size(tmt::LazyTensorMappingTranspose) = domain_size(tmt.tm)
 domain_size(tmt::LazyTensorMappingTranspose) = range_size(tmt.tm)
@@ -67,8 +80,8 @@
 end
 # TODO: Boundschecking in constructor.
 
-apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...)
-apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...)
+apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...)
+apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...)
 
 range_size(tmBinOp::LazyTensorMappingBinaryOperation) = range_size(tmBinOp.tm1)
 domain_size(tmBinOp::LazyTensorMappingBinaryOperation) = domain_size(tmBinOp.tm1)
@@ -90,16 +103,15 @@
         return new{T,R,K,D, typeof(t1), typeof(t2)}(t1,t2)
     end
 end
-export TensorMappingComposition
 
 range_size(tm::TensorMappingComposition) = range_size(tm.t1)
 domain_size(tm::TensorMappingComposition) = domain_size(tm.t2)
 
-function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,K,D}
+function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,K,D}
     apply(c.t1, c.t2*v, I...)
 end
 
-function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,K,D}
+function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,K,D}
     apply_transpose(c.t2, c.t1'*v, I...)
 end
 
@@ -127,12 +139,11 @@
         return new{T,R,D,RD,AA}(A,range_indicies,domain_indicies)
     end
 end
-export LazyLinearMap
 
 range_size(llm::LazyLinearMap) = size(llm.A)[[llm.range_indicies...]]
 domain_size(llm::LazyLinearMap) = size(llm.A)[[llm.domain_indicies...]]
 
-function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D}
+function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D}
     view_index = ntuple(i->:,ndims(llm.A))
     for i ∈ 1:R
         view_index = Base.setindex(view_index, Int(I[i]), llm.range_indicies[i])
@@ -141,7 +152,7 @@
     return sum(A_view.*v)
 end
 
-function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D}
+function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D}
     apply(LazyLinearMap(llm.A, llm.domain_indicies, llm.range_indicies), v, I...)
 end
 
@@ -155,7 +166,6 @@
 struct IdentityMapping{T,D} <: TensorMapping{T,D,D}
     size::NTuple{D,Int}
 end
-export IdentityMapping
 
 IdentityMapping{T}(size::NTuple{D,Int}) where {T,D} = IdentityMapping{T,D}(size)
 IdentityMapping{T}(size::Vararg{Int,D}) where {T,D} = IdentityMapping{T,D}(size)
@@ -164,8 +174,8 @@
 range_size(tmi::IdentityMapping) = tmi.size
 domain_size(tmi::IdentityMapping) = tmi.size
 
-apply(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...]
-apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...]
+apply(tmi::IdentityMapping{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = v[I...]
+apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = v[I...]
 
 """
     Base.:∘(tm, tmi)
@@ -212,7 +222,6 @@
         return new{T,R,D,D_before,R_middle,D_middle,D_after, typeof(tm)}(before, tm, after)
     end
 end
-export InflatedTensorMapping
 """
     InflatedTensorMapping(before, tm, after)
     InflatedTensorMapping(before,tm)
@@ -258,7 +267,7 @@
     )
 end
 
-function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D}
+function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D}
     dim_before = range_dim(itm.before)
     dim_domain = domain_dim(itm.tm)
     dim_range = range_dim(itm.tm)
@@ -270,7 +279,7 @@
     return apply(itm.tm, v_inner, inner_index...)
 end
 
-function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D}
+function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D}
     dim_before = range_dim(itm.before)
     dim_domain = domain_dim(itm.tm)
     dim_range = range_dim(itm.tm)
@@ -398,7 +407,6 @@
 ```
 """
 function LazyOuterProduct end
-export LazyOuterProduct
 
 function LazyOuterProduct(tm1::TensorMapping{T}, tm2::TensorMapping{T}) where T
     itm1 = InflatedTensorMapping(tm1, IdentityMapping{T}(range_size(tm2)))
@@ -414,7 +422,6 @@
 LazyOuterProduct(tms::Vararg{TensorMapping}) = foldl(LazyOuterProduct, tms)
 
 ⊗(a::TensorMapping, b::TensorMapping) = LazyOuterProduct(a,b)
-export ⊗
 
 
 function check_domain_size(tm::TensorMapping, sz)
@@ -427,7 +434,6 @@
     tm::TensorMapping
     sz
 end
-export SizeMismatch
 
 function Base.showerror(io::IO, err::SizeMismatch)
     print(io, "SizeMismatch: ")

--- a/src/LazyTensors/tensor_mapping.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/src/LazyTensors/tensor_mapping.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -1,10 +1,16 @@
+export TensorMapping
+export apply
+export apply_transpose
+export range_dim, domain_dim
+export range_size, domain_size
+
 """
     TensorMapping{T,R,D}
 
 Describes a mapping of a `D` dimension tensor to an `R` dimension tensor.
 The action of the mapping is implemented through the method
 ```julia
-    apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+    apply(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg) where {R,D,T}
 ```
 
 The size of the range and domain that the operator works with should be returned by
@@ -21,23 +27,20 @@
 ```
 """
 abstract type TensorMapping{T,R,D} end
-export TensorMapping
 
 """
-    apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
+    apply(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg) where {R,D,T}
 
 Return the result of the mapping for a given index.
 """
 function apply end
-export apply
 
 """
-    apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg) where {R,D,T}
+    apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg) where {R,D,T}
 
 Return the result of the transposed mapping for a given index.
 """
 function apply_transpose end
-export apply_transpose
 
 """
     range_dim(::TensorMapping)
@@ -51,7 +54,6 @@
 """
 domain_dim(::TensorMapping{T,R,D}) where {T,R,D} = D
 
-export range_dim, domain_dim
 
 """
     range_size(M::TensorMapping)
@@ -67,7 +69,6 @@
 """
 function domain_size end
 
-export range_size, domain_size
 
 """
     eltype(::TensorMapping{T})

--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -62,27 +62,27 @@
 LazyTensors.range_size(op::BoundaryOperator) = ()
 LazyTensors.domain_size(op::BoundaryOperator) = (op.size,)
 
-function LazyTensors.apply(op::BoundaryOperator{T,Lower}, v::AbstractVector{T}) where T
+function LazyTensors.apply(op::BoundaryOperator{<:Any,Lower}, v::AbstractVector)
     apply_stencil(op.stencil,v,1)
 end
 
-function LazyTensors.apply(op::BoundaryOperator{T,Upper}, v::AbstractVector{T}) where T
+function LazyTensors.apply(op::BoundaryOperator{<:Any,Upper}, v::AbstractVector)
     apply_stencil_backwards(op.stencil,v,op.size)
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T
+function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Lower}, v::AbstractArray{<:Any,0}, i::Index{Lower})
     return op.stencil[Int(i)-1]*v[]
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T
+function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Upper}, v::AbstractArray{<:Any,0}, i::Index{Upper})
     return op.stencil[op.size[1] - Int(i)]*v[]
 end
 
 # Catch all combinations of Lower, Upper and Interior not caught by the two previous methods.
-function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i::Index) where T
-    return zero(T)
+function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i::Index)
+    return zero(eltype(v))
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T
+function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i)
     return LazyTensors.apply_transpose_with_region(op, v, closure_size(op), op.size[1], i)
 end

--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -28,19 +28,19 @@
 LazyTensors.domain_size(op::ConstantInteriorScalingOperator) = (op.size,)
 
 # TBD: @inbounds in apply methods?
-function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Lower})
     return op.closure_weights[Int(i)]*v[Int(i)]
 end
 
-function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Interior})
     return op.interior_weight*v[Int(i)]
 end
 
-function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T
+function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Upper})
     return op.closure_weights[op.size[1]-Int(i)+1]*v[Int(i)]
 end
 
-function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T
+function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i)
     return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i)
 end
 

--- a/src/SbpOperators/volumeops/volume_operator.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -43,18 +43,18 @@
 LazyTensors.range_size(op::VolumeOperator) = op.size
 LazyTensors.domain_size(op::VolumeOperator) = op.size
 
-function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Lower})
     return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i))
 end
 
-function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Interior})
     return apply_stencil(op.inner_stencil, v, Int(i))
 end
 
-function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T
+function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Upper})
     return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i))
 end
 
-function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T
+function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i)
     return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i)
 end

--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -80,6 +80,28 @@
     v = [[1 2];[3 4]]
     @test m*v == [[2 4];[6 8]]
     @test (m*v)[2,1] == 6
+
+    @testset "Promotion" begin
+        m = ScalingOperator{Int,1}(2,(3,))
+        v = [1.,2.,3.]
+        @test m*v isa AbstractVector{Float64}
+        @test m*v == [2.,4.,6.]
+
+        m = ScalingOperator{Int,2}(2,(2,2))
+        v = [[1. 2.];[3. 4.]]
+        @test m*v == [[2. 4.];[6. 8.]]
+        @test (m*v)[2,1] == 6.
+
+        m = ScalingOperator{ComplexF64,1}(2. +2. *im,(3,))
+        v = [1.,2.,3.]
+        @test m*v isa AbstractVector{ComplexF64}
+        @test m*v == [2. + 2. *im, 4. + 4. *im, 6. + 6. *im]
+
+        m = ScalingOperator{ComplexF64,1}(1,(3,))
+        v = [2. + 2. *im, 4. + 4. *im, 6. + 6. *im]
+        @test m*v isa AbstractVector{ComplexF64}
+        @test m*v == [2. + 2. *im, 4. + 4. *im, 6. + 6. *im]
+    end
 end
 
 @testset "TensorMapping binary operations" begin
@@ -107,6 +129,8 @@
 
     @test range_size(A+B) == range_size(A) == range_size(B)
     @test domain_size(A+B) == domain_size(A) == domain_size(B)
+
+    @test ((A+B)*ComplexF64[1.1,1.2,1.3])[3] isa ComplexF64
 end
 
 
@@ -129,6 +153,9 @@
 
     v = rand(2)
     @test (Ã∘B̃)'*v ≈ B'*A'*v rtol=1e-14
+
+    @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64
+    @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64
 end
 
 @testset "LazyLinearMap" begin
@@ -194,6 +221,10 @@
         @test I*v == v
         @test I'*v == v
 
+        v = rand(ComplexF64,sz...)
+        @test I*v == v
+        @test I'*v == v
+
         @test range_size(I) == sz
         @test domain_size(I) == sz
     end
@@ -322,6 +353,16 @@
             end
         end
 
+        @testset "application to other type" begin
+            tm = InflatedTensorMapping(I(3,2), A, I(4))
+
+            v = rand(ComplexF64, domain_size(tm)...)
+            @test (tm*v)[1,2,3,1] isa ComplexF64
+
+            v = rand(ComplexF64, domain_size(tm')...)
+            @test (tm'*v)[1,2,2,1] isa ComplexF64
+        end
+
         @testset "Inference of application" begin
             struct ScalingOperator{T,D} <: TensorMapping{T,D,D}
                 λ::T

--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -72,6 +72,14 @@
             @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2]
             @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
             @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
+
+            v = evalOn(g_1D, x->1. +x*im)
+            @test (op_l*v)[] isa ComplexF64
+
+            u = fill(1. +im)
+            @test (op_l'*u)[1] isa ComplexF64
+            @test (op_l'*u)[5] isa ComplexF64
+            @test (op_l'*u)[11] isa ComplexF64
         end
 
         @testset "2D" begin

--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -24,6 +24,9 @@
     @test a*v == [.1,.2,.5,.5,.5,.2,.1]
     @test a'*v == [.1,.2,.5,.5,.5,.2,.1]
 
+    @test (a*rand(ComplexF64, domain_size(a)... ))[1] isa ComplexF64
+    @test (a'*rand(ComplexF64, domain_size(a')...))[1] isa ComplexF64
+
     @test range_size(a) == (7,)
     @test domain_size(a) == (7,)
 

--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Mon Mar 14 09:46:22 2022 +0100
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Mon Mar 14 10:14:38 2022 +0100
@@ -89,6 +89,8 @@
     @testset "Application" begin
         @test op_x*v ≈ rx rtol = 1e-14
         @test op_y*v ≈ ry rtol = 1e-14
+
+        @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64
     end
 
     @testset "Regions" begin