diff --git a/examples/calderon.jl b/examples/calderon.jl
index 4e2d6084..5e2631c2 100644
--- a/examples/calderon.jl
+++ b/examples/calderon.jl
@@ -1,24 +1,45 @@
-# If you want a scatter plot of the spectra, define `plotresults = true`
-# prior to running this script.
-
 using CompScienceMeshes, BEAST
 using LinearAlgebra
 
 Γ = readmesh(joinpath(dirname(pathof(BEAST)),"../examples/sphere2.in"))
-println("Mesh with $(numvertices(Γ)) vertices and $(numcells(Γ)) cells.")
 X = raviartthomas(Γ)
-Y = BEAST.buffachristiansen2(Γ)
+Y = BEAST.buffachristiansen(Γ)
 
-κ = ω = 1.0; γ = κ*im
-T = MWSingleLayer3D(γ)
+κ = 1.0;
+T = Maxwell3D.singlelayer(wavenumber=κ)
 N = NCross()
 
-Txx = assemble(T,X,X); println("primal discretisation assembled.")
+E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
+e = (n × E) × n
+
+b = assemble(e, X)
+
+A = assemble(T,X,X); println("primal discretisation assembled.")
+
 Tyy = assemble(T,Y,Y); println("dual discretisation assembled.")
 Nxy = Matrix(assemble(N,X,Y)); println("duality form assembled.")
-
 iNxy = inv(Nxy); println("duality form inverted.")
-A = iNxy' * Tyy * iNxy * Txx
+P = iNxy' * Tyy * iNxy
+
+iT1 = BEAST.GMRESSolver(A; restart=1_500, abstol=1e-8, reltol=1e-8, maxiter=1_500)
+iT2 = BEAST.GMRESSolver(A; restart=1_500, abstol=1e-8, reltol=1e-8, maxiter=1_500, left_preconditioner=P)
+
+u1, ch1 = BEAST.solve(iT1, bx)
+u2, ch2 = BEAST.solve(iT2, bx)
+
+using Plots
+Plots.plot(title="log10 residual error vs iteration count")
+Plots.plot!(log10.(ch1[:resnorm]), label="A=b", l=2)
+Plots.plot!(log10.(ch2[:resnorm]), label="P*A=P*b", l=2)
+
+ss1 = svdvals(Txx)
+ss2 = svdvals(P*Txx)
+
+Plots.plot(title="singular value spectrum")
+Plots.plot!(ss1, label="A", l=2)
+Plots.plot!(ss2, label="P*A", l=2)
+
+@show norm(u1-u2) / norm(u1)
 
 @show cond(Txx)
 @show cond(Tyy)
diff --git a/src/bases/global/gwpglobal.jl b/src/bases/global/gwpglobal.jl
index 3270e4da..f302632f 100644
--- a/src/bases/global/gwpglobal.jl
+++ b/src/bases/global/gwpglobal.jl
@@ -147,7 +147,7 @@ function gwpcurl(mesh, edges=nothing; order)
             _addshapes!(fns, cell, gids, lids, inv(α'))
         end
 
-        order < 2 && continue
+        order < 1 && continue
         face_ch = chart(mesh, cell)
         gids = Ne + (cell-1)*nf .+ (1:nf)
         lids = localindices(localspace, refdom, Val{2}, 1)
diff --git a/src/operator.jl b/src/operator.jl
index b54b2a74..22dbb102 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -89,7 +89,7 @@ function assemble(operator::AbstractOperator, test_functions, trial_functions;
 end
 
 
-function assemble(A::Matrix, testfns, trialfns)
+function assemble(A::AbstractMatrix, testfns, trialfns)
     @assert numfunctions(testfns) == size(A,1)
     @assert numfunctions(trialfns) == size(A,2)
     return A
@@ -137,15 +137,15 @@ function allocatestorage(operator::AbstractOperator, test_functions, trial_funct
 end
 
 
-function allocatestorage(operator::LinearCombinationOfOperators,
-        test_functions::SpaceTimeBasis, trial_functions::SpaceTimeBasis,
-        storage_policy::Type{Val{:bandedstorage}},
-        long_delays_policy::Type{LongDelays{:ignore}})
+# function allocatestorage(operator::LinearCombinationOfOperators,
+#         test_functions::SpaceTimeBasis, trial_functions::SpaceTimeBasis,
+#         storage_policy::Type{Val{:bandedstorage}},
+#         long_delays_policy::Type{LongDelays{:ignore}})
 
-    # This works when are terms in the LC can share storage
-    return allocatestorage(operator.ops[end], test_functions, trial_functions,
-        storage_policy, long_delays_policy)
-end
+#     # This works when are terms in the LC can share storage
+#     return allocatestorage(operator.ops[end], test_functions, trial_functions,
+#         storage_policy, long_delays_policy)
+# end
 
 
 function allocatestorage(operator::LinearCombinationOfOperators,
diff --git a/src/utils/variational.jl b/src/utils/variational.jl
index aa0a5158..7bb9930a 100644
--- a/src/utils/variational.jl
+++ b/src/utils/variational.jl
@@ -286,6 +286,11 @@ function getindex(A, v::HilbertVector, u::HilbertVector)
     BilForm(v.space, u.space, terms)
 end
 
+function getindex(A::AbstractMatrix, v::HilbertVector, u::HilbertVector)
+    terms = [ BilTerm(v.idx, u.idx, v.opstack, u.opstack, 1, A) ]
+    BilForm(v.space, u.space, terms)
+end
+
 
 function getindex(A::BEAST.BlockDiagonalOperator, V::Vector{HilbertVector}, U::Vector{HilbertVector})
     op = A.op