diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index c95e52d..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,26 +0,0 @@
-# Documentation: http://docs.travis-ci.com/user/languages/julia
-language: julia
-notifications:
-  email: false
-julia:
-  - 1.6
-  - 1.8
-  - nightly
-os:
-  - linux
-arch:
-  - x64
-cache:
-  directories:
-    - ~/.julia/artifacts
-jobs:
-  fast_finish: true
-  allow_failures:
-    - julia: nightly
-after_success:
-  - |
-    julia -e '
-      using Pkg
-      Pkg.add("Coverage")
-      using Coverage
-      Codecov.submit(process_folder())'
diff --git a/src/RadialBasisFunctions.jl b/src/RadialBasisFunctions.jl
index 508712f..ebf9963 100644
--- a/src/RadialBasisFunctions.jl
+++ b/src/RadialBasisFunctions.jl
@@ -16,6 +16,7 @@ export AbstractPHS, PHS, PHS1, PHS3, PHS5, PHS7
 export IMQ
 export Gaussian
 export MonomialBasis
+export degree, dim
 
 include("utils.jl")
 export find_neighbors, reorder_points!
@@ -42,7 +43,7 @@ export ∂virtual
 
 include("operators/monomial/monomial.jl")
 
-include("operators/operator_combinations.jl")
+include("operators/operator_algebra.jl")
 
 include("interpolation.jl")
 export Interpolator
@@ -52,7 +53,7 @@ export Regrid, regrid
 
 # Some consts and aliases
 const Δ = ∇² # some people like this notation for the Laplacian
-const AVOID_NAN = 1e-16
+const AVOID_INF = 1e-16
 
 using PrecompileTools
 @setup_workload begin
diff --git a/src/basis/basis.jl b/src/basis/basis.jl
index 4ab9ee5..47858d6 100644
--- a/src/basis/basis.jl
+++ b/src/basis/basis.jl
@@ -29,6 +29,9 @@ function (ℒmon::ℒMonomialBasis{Dim,Deg})(x) where {Dim,Deg}
 end
 (m::ℒMonomialBasis)(b, x) = m.f(b, x)
 
+degree(::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg} = Deg
+dim(::ℒMonomialBasis{Dim,Deg}) where {Dim,Deg} = Dim
+
 include("polyharmonic_spline.jl")
 include("inverse_multiquadric.jl")
 include("gaussian.jl")
diff --git a/src/basis/polyharmonic_spline.jl b/src/basis/polyharmonic_spline.jl
index 19b18ce..d1b2b60 100644
--- a/src/basis/polyharmonic_spline.jl
+++ b/src/basis/polyharmonic_spline.jl
@@ -39,7 +39,7 @@ end
 
 (phs::PHS1)(x, xᵢ) = euclidean(x, xᵢ)
 function ∂(::PHS1, dim::Int)
-    ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + AVOID_NAN)
+    ∂ℒ(x, xᵢ) = (x[dim] - xᵢ[dim]) / (euclidean(x, xᵢ) + AVOID_INF)
     return ℒRadialBasisFunction(∂ℒ)
 end
 function ∇(::PHS1)
@@ -49,14 +49,14 @@ end
 function ∂²(::PHS1, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return (-(x[dim] - xᵢ[dim])^2 + sqeuclidean(x, xᵢ)) /
-               (euclidean(x, xᵢ)^3 + AVOID_NAN)
+               (euclidean(x, xᵢ)^3 + AVOID_INF)
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS1)
     function ∇²ℒ(x, xᵢ)
         return sum(
-            (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + AVOID_NAN)
+            (-(x .- xᵢ) .^ 2 .+ sqeuclidean(x, xᵢ)) / (euclidean(x, xᵢ)^3 + AVOID_INF)
         )
     end
     return ℒRadialBasisFunction(∇²ℒ)
@@ -87,14 +87,14 @@ end
 function ∂²(::PHS3, dim::Int)
     function ∂²ℒ(x, xᵢ)
         return 3 * (sqeuclidean(x, xᵢ) + (x[dim] - xᵢ[dim])^2) /
-               (euclidean(x, xᵢ) + AVOID_NAN)
+               (euclidean(x, xᵢ) + AVOID_INF)
     end
     return ℒRadialBasisFunction(∂²ℒ)
 end
 function ∇²(::PHS3)
     function ∇²ℒ(x, xᵢ)
         return sum(
-            3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + AVOID_NAN)
+            3 * (sqeuclidean(x, xᵢ) .+ (x .- xᵢ) .^ 2) / (euclidean(x, xᵢ) + AVOID_INF)
         )
     end
     return ℒRadialBasisFunction(∇²ℒ)
diff --git a/src/operators/directional.jl b/src/operators/directional.jl
index 85acdfe..9599416 100644
--- a/src/operators/directional.jl
+++ b/src/operators/directional.jl
@@ -14,12 +14,12 @@ end
 Builds a `RadialBasisOperator` where the operator is the directional derivative, `Directional`.
 """
 function directional(
-    data::AbstractVector{D},
+    data::AbstractVector,
     v::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
     ℒ = Directional(f, v)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
@@ -45,14 +45,12 @@ end
 
 function RadialBasisOperator(
     ℒ::Directional,
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = spzeros(eltype(D), Na, Nd)
+) where {T<:Int,B<:AbstractRadialBasis}
+    weights = _build_weights(ℒ, data, data, adjl, basis)
     return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
 end
 
@@ -64,12 +62,26 @@ function RadialBasisOperator(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = spzeros(eltype(TD), Na, Nd)
+    weights = _build_weights(ℒ, data, eval_points, adjl, basis)
     return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
 end
 
+function _build_weights(ℒ::Directional, data, eval_points, adjl, basis)
+    v = ℒ.v
+    N = length(first(data))
+    @assert length(v) == N || length(v) == length(data) "wrong size for v"
+    if length(v) == N
+        return mapreduce(+, enumerate(ℒ.ℒ)) do (i, ℒ)
+            _build_weights(ℒ, data, eval_points, adjl, basis) * v[i]
+        end
+    else
+        vv = ntuple(i -> getindex.(v, i), N)
+        return mapreduce(+, enumerate(ℒ.ℒ)) do (i, ℒ)
+            Diagonal(vv[i]) * _build_weights(ℒ, data, eval_points, adjl, basis)
+        end
+    end
+end
+
 function update_weights!(op::RadialBasisOperator{<:Directional})
     v = op.ℒ.v
     N = length(first(op.data))
@@ -88,7 +100,5 @@ function update_weights!(op::RadialBasisOperator{<:Directional})
     return nothing
 end
 
-Base.size(op::RadialBasisOperator{<:Directional}) = size(op.weights)
-
 # pretty printing
 print_op(op::Directional) = "Directional Gradient (∇f⋅v)"
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
index 77a6cf9..86fdedc 100644
--- a/src/operators/gradient.jl
+++ b/src/operators/gradient.jl
@@ -14,11 +14,11 @@ end
 Builds a `RadialBasisOperator` where the operator is the gradient, `Gradient`.
 """
 function gradient(
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,B<:AbstractRadialBasis,T<:Int}
+) where {B<:AbstractRadialBasis,T<:Int}
     f = ntuple(dim -> Base.Fix2(∂, dim), length(first(data)))
     ℒ = Gradient(f)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
diff --git a/src/operators/laplacian.jl b/src/operators/laplacian.jl
index 8fedacf..9e268f5 100644
--- a/src/operators/laplacian.jl
+++ b/src/operators/laplacian.jl
@@ -9,11 +9,11 @@ end
 
 # convienience constructors
 function laplacian(
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {T<:Int,B<:AbstractRadialBasis}
     ℒ = Laplacian(∇²)
     return RadialBasisOperator(ℒ, data, basis; k=k, adjl=adjl)
 end
diff --git a/src/operators/monomial/monomial.jl b/src/operators/monomial/monomial.jl
index 4e1425f..1068589 100644
--- a/src/operators/monomial/monomial.jl
+++ b/src/operators/monomial/monomial.jl
@@ -15,10 +15,9 @@ end
 function ∇²(m::MonomialBasis{Dim,Deg}) where {Dim,Deg}
     ∂² = ntuple(dim -> ∂(m, 2, dim), Dim)
     function basis!(b, x)
-        cache = ones(size(b))
-        b .= 0
+        cache = ones(eltype(x), size(b))
+        b .= zero(eltype(x))
         for ∂²! in ∂²
-            # use mapreduce here instead?
             ∂²!(cache, x)
             b .+= cache
         end
@@ -34,8 +33,8 @@ end
 
 function build_monomial_basis(ids::Vector{Vector{Vector{T}}}, c::Vector{T}) where {T<:Int}
     function basis!(db::AbstractVector{B}, x::AbstractVector) where {B}
-        db .= 1
-        # TODO flatten loop - why does it allocate here
+        db .= one(eltype(x))
+        # TODO optimize - allocations
         @views @inbounds for i in eachindex(ids), j in eachindex(ids[i])
             db[ids[i][j]] *= x[i]
         end
diff --git a/src/operators/operator_algebra.jl b/src/operators/operator_algebra.jl
new file mode 100644
index 0000000..98d8fef
--- /dev/null
+++ b/src/operators/operator_algebra.jl
@@ -0,0 +1,45 @@
+for op in (:+, :-)
+    @eval function Base.$op(a::ℒRadialBasisFunction, b::ℒRadialBasisFunction)
+        additive_ℒRBF(x, xᵢ) = Base.$op(a(x, xᵢ), b(x, xᵢ))
+        return ℒRadialBasisFunction(additive_ℒRBF)
+    end
+end
+
+for op in (:+, :-)
+    @eval function Base.$op(
+        a::ℒMonomialBasis{Dim,Deg}, b::ℒMonomialBasis{Dim,Deg}
+    ) where {Dim,Deg}
+        function additive_ℒMon(m, x)
+            m .= Base.$op.(a(x), b(x))
+            return nothing
+        end
+        return ℒMonomialBasis(Dim, Deg, additive_ℒMon)
+    end
+end
+
+for op in (:+, :-)
+    @eval function Base.$op(op1::RadialBasisOperator, op2::RadialBasisOperator)
+        _check_compatible(op1, op2)
+        k = _update_stencil(op1, op2)
+        ℒ(x) = Base.$op(op1.ℒ(x), op2.ℒ(x))
+        return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k, adjl=op1.adjl)
+    end
+end
+
+function _check_compatible(op1::RadialBasisOperator, op2::RadialBasisOperator)
+    if !all(op1.data .≈ op2.data)
+        throw(
+            ArgumentError("Can not add operators that were not built with the same data.")
+        )
+    end
+    if !all(op1.adjl .≈ op2.adjl)
+        throw(ArgumentError("Can not add operators that do not have the same stencils."))
+    end
+end
+
+function _update_stencil(op1::RadialBasisOperator, op2::RadialBasisOperator)
+    k1 = length(first((op1.adjl)))
+    k2 = length(first((op2.adjl)))
+    k = k1 > k2 ? k1 : k2
+    return k
+end
diff --git a/src/operators/operator_combinations.jl b/src/operators/operator_combinations.jl
deleted file mode 100644
index 6d5076d..0000000
--- a/src/operators/operator_combinations.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-for op in (:+, :-, :*, :/)
-    @eval function Base.$op(a::ℒRadialBasisFunction, b::ℒRadialBasisFunction)
-        return additive_ℒRBF(x, xᵢ) = Base.$op(a(x, xᵢ), b(x, xᵢ))
-    end
-end
-
-for op in (:+, :-, :*, :/)
-    @eval function Base.$op(a::ℒMonomialBasis, b::ℒMonomialBasis)
-        function additive_ℒMon(m, x)
-            cache = ones(size(m))
-            m .= 0
-            a(cache, x)
-            m .+= cache
-            b(cache, x)
-            @. m = Base.$op(m, cache)
-            return nothing
-        end
-    end
-end
-
-for op in (:+, :-, :*, :/)
-    @eval function Base.$op(op1::RadialBasisOperator, op2::RadialBasisOperator)
-        if !all(op1.data .≈ op2.data)
-            throw(
-                ArgumentError(
-                    "Can not add operators that were not built with the same data."
-                ),
-            )
-        end
-        if !all(op1.adjl .≈ op2.adjl)
-            throw(
-                ArgumentError("Can not add operators that do not have the same stencils.")
-            )
-        end
-        k1 = length(first((op1.adjl)))
-        k2 = length(first((op2.adjl)))
-        k = k1 > k2 ? k1 : k2
-        ℒ = Base.$op(op1.ℒ, op2.ℒ)
-        return RadialBasisOperator(ℒ, op1.data, op1.basis; k=k, adjl=op1.adjl)
-    end
-end
-
-for op in (:+, :-, :*, :/)
-    @eval function Base.$op(a::ScalarValuedOperator, b::ScalarValuedOperator)
-        order = (a.order..., b.order)
-        dim = (a.dim..., b.dim)
-        f(x) = Base.$op(a.ℒ(x), a.ℒ(x))
-        return Partial(f, order, dim)
-    end
-end
diff --git a/src/operators/operators.jl b/src/operators/operators.jl
index 031100b..10adfb8 100644
--- a/src/operators/operators.jl
+++ b/src/operators/operators.jl
@@ -2,6 +2,8 @@ abstract type AbstractOperator end
 abstract type ScalarValuedOperator <: AbstractOperator end
 abstract type VectorValuedOperator <: AbstractOperator end
 
+(op::AbstractOperator)(x) = op.ℒ(x)
+
 """
     struct RadialBasisOperator
 
@@ -16,24 +18,30 @@ struct RadialBasisOperator{L,W,D,C,A,B<:AbstractRadialBasis}
     basis::B
     valid_cache::Base.RefValue{Bool}
     function RadialBasisOperator(
-        ℒ::L, weights::W, data::D, eval_points::C, adjl::A, basis::B
+        ℒ::L,
+        weights::W,
+        data::D,
+        eval_points::C,
+        adjl::A,
+        basis::B,
+        cache_status::Bool=false,
     ) where {L,W,D,C,A,B<:AbstractRadialBasis}
-        return new{L,W,D,C,A,B}(ℒ, weights, data, eval_points, adjl, basis, Ref(false))
+        return new{L,W,D,C,A,B}(
+            ℒ, weights, data, eval_points, adjl, basis, Ref(cache_status)
+        )
     end
 end
 
 # convienience constructors
 function RadialBasisOperator(
     ℒ,
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = spzeros(eltype(D), Na, Nd)
-    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
+) where {T<:Int,B<:AbstractRadialBasis}
+    weights = _build_weights(ℒ, data, data, adjl, basis)
+    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis, true)
 end
 
 function RadialBasisOperator(
@@ -44,23 +52,19 @@ function RadialBasisOperator(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = spzeros(eltype(TD), Na, Nd)
-    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
+    weights = _build_weights(ℒ, data, eval_points, adjl, basis)
+    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis, true)
 end
 
 function RadialBasisOperator(
     ℒ::VectorValuedOperator,
-    data::AbstractVector{D},
+    data::AbstractVector,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
-    TD = eltype(D)
-    N = length(adjl)
-    weights = ntuple(_ -> _allocate_weights(TD, N, N, k), length(ℒ.ℒ))
-    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis)
+) where {T<:Int,B<:AbstractRadialBasis}
+    weights = ntuple(i -> _build_weights(ℒ.ℒ[i], data, data, adjl, basis), length(ℒ.ℒ))
+    return RadialBasisOperator(ℒ, weights, data, data, adjl, basis, true)
 end
 
 function RadialBasisOperator(
@@ -71,15 +75,16 @@ function RadialBasisOperator(
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, eval_points, k),
 ) where {TD,TE,T<:Int,B<:AbstractRadialBasis}
-    Na = length(adjl)
-    Nd = length(data)
-    weights = ntuple(_ -> _allocate_weights(eltype(TD), Na, Nd, k), length(ℒ.ℒ))
-    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis)
+    weights = ntuple(length(ℒ.ℒ)) do i
+        _build_weights(ℒ.ℒ[i], data, eval_points, adjl, basis)
+    end
+    return RadialBasisOperator(ℒ, weights, data, eval_points, adjl, basis, true)
 end
 
 # extend Base methods
 Base.length(op::RadialBasisOperator) = length(op.adjl)
 Base.size(op::RadialBasisOperator) = size(op.weights)
+Base.size(op::RadialBasisOperator, dim::Int) = size(op.weights, dim)
 function Base.size(op::RadialBasisOperator{<:VectorValuedOperator})
     return ntuple(i -> size(op.weights[i]), embeddim(op))
 end
diff --git a/src/operators/partial.jl b/src/operators/partial.jl
index 7208a8f..d9c0775 100644
--- a/src/operators/partial.jl
+++ b/src/operators/partial.jl
@@ -16,13 +16,13 @@ end
 Builds a `RadialBasisOperator` where the operator is the partial derivative, `Partial`, of `order` with respect to `dim`.
 """
 function partial(
-    data::AbstractVector{D},
+    data::AbstractVector,
     order::T,
     dim::T,
     basis::B=PHS(3; poly_deg=2);
     k::T=autoselect_k(data, basis),
     adjl=find_neighbors(data, k),
-) where {D<:AbstractArray,T<:Int,B<:AbstractRadialBasis}
+) where {T<:Int,B<:AbstractRadialBasis}
     f = let o = order, dim = dim
         x -> ∂(x, o, dim)
     end
diff --git a/src/operators/regridding.jl b/src/operators/regridding.jl
index 01639ff..11dae54 100644
--- a/src/operators/regridding.jl
+++ b/src/operators/regridding.jl
@@ -7,6 +7,7 @@ struct Regrid
     ℒ
     Regrid() = new(identity)
 end
+(op::Regrid)(x) = op.ℒ(x)
 
 # convienience constructors
 """
diff --git a/src/solve.jl b/src/solve.jl
index cb54f99..fadf64a 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -1,6 +1,12 @@
-function _build_weights(ℒ, op; nchunks=Threads.nthreads())
+function _build_weights(ℒ, op)
     data = op.data
+    eval_points = op.eval_points
+    adjl = op.adjl
     basis = op.basis
+    return _build_weights(ℒ, data, eval_points, adjl, basis)
+end
+
+function _build_weights(ℒ, data, eval_points, adjl, basis)
     dim = length(first(data)) # dimension of data
 
     # build monomial basis and operator
@@ -8,15 +14,11 @@ function _build_weights(ℒ, op; nchunks=Threads.nthreads())
     ℒmon = ℒ(mon)
     ℒrbf = ℒ(basis)
 
-    return _build_weights(op, ℒrbf, ℒmon, mon; nchunks=nchunks)
+    return _build_weights(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon)
 end
 
-function _build_weights(op, ℒrbf, ℒmon, mon; nchunks=Threads.nthreads())
-    data = op.data
-    eval_points = op.eval_points
-    adjl = op.adjl
-    basis = op.basis
-
+function _build_weights(data, eval_points, adjl, basis, ℒrbf, ℒmon, mon)
+    nchunks = Threads.nthreads()
     TD = eltype(first(data))
     dim = length(first(data)) # dimension of data
     nmon = binomial(dim + basis.poly_deg, basis.poly_deg)
@@ -66,8 +68,8 @@ function _build_stencil!(
 end
 
 function _build_collocation_matrix!(
-    A::Symmetric, data::AbstractVector{D}, basis::B, mon::MonomialBasis{Dim,Deg}, k::K
-) where {D<:AbstractArray,B<:AbstractRadialBasis,K<:Int,Dim,Deg}
+    A::Symmetric, data::AbstractVector, basis::B, mon::MonomialBasis{Dim,Deg}, k::K
+) where {B<:AbstractRadialBasis,K<:Int,Dim,Deg}
     # radial basis section
     AA = parent(A)
     N = size(A, 2)
diff --git a/test/operators/operator_algebra.jl b/test/operators/operator_algebra.jl
new file mode 100644
index 0000000..810d428
--- /dev/null
+++ b/test/operators/operator_algebra.jl
@@ -0,0 +1,24 @@
+using RadialBasisFunctions
+using StaticArrays
+using LinearAlgebra
+using Statistics
+using HaltonSequences
+
+mean_percent_error(test, correct) = mean(abs.((test .- correct) ./ correct)) * 100
+
+f(x) = 2 * x[1] + 3 * x[2]
+df_dx(x) = 2
+df_dy(x) = 3
+
+N = 1000
+x = SVector{2}.(HaltonPoint(2)[1:N])
+y = f.(x)
+
+dx = partial(x, 1, 1)
+dy = partial(x, 1, 2)
+
+dxdy = dx + dy
+@test mean_percent_error(dxdy(y), df_dx.(x) .+ df_dy.(x)) < 1e-6
+
+dxdy = dx - dy
+@test mean_percent_error(dxdy(y), df_dx.(x) .- df_dy.(x)) < 1e-6
diff --git a/test/runtests.jl b/test/runtests.jl
index 9314673..759eb83 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -44,6 +44,10 @@ end
     include("operators/virtual.jl")
 end
 
+@safetestset "Operator Algebra" begin
+    include("operators/operator_algebra.jl")
+end
+
 @safetestset "Stencil" begin
     include("solve.jl")
 end