diff --git a/test/runtests.jl b/test/runtests.jl
index c1499fb..1b748cd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -45,264 +45,272 @@ for m in [1, 3]
 end
 =#
 
+@testset "KroneckerTools" verbose = true begin
+    @testset "ExtendedMul" begin
+        ma = 2
+        na = 4
+        a = randn(ma, na)
+        mb = 4
+        nb = 3
+        b = randn(mb, nb)
+        c = randn(ma, nb)
+        A_mul_B!(c, 1, a, 1, ma, na, b, 1, nb)
+        @test c ≈ a * b
+        mb = 2
+        nb = 3
+        b = randn(mb, nb)
+        c = randn(na, nb)
+        At_mul_B!(c, 1, a, 1, na, ma, b, 1, nb)
+        @test c ≈ transpose(a) * b
+        nb = 4
+        mb = 3
+        b = randn(mb, nb)
+        c = randn(ma, mb)
+        A_mul_Bt!(c, 1, a, 1, ma, na, b, 1, mb)
+        @test c ≈ a * transpose(b)
+        mb = 4
+        nb = 2
+        b = randn(mb, nb)
+        c = randn(na, mb)
+        At_mul_Bt!(c, 1, a, 1, na, ma, b, 1, mb)
+        @test c ≈ transpose(a) * transpose(b)
+    end
 
-ma = 2
-na = 4
-a = randn(ma, na)
-mb = 4
-nb = 3
-b = randn(mb, nb)
-c = randn(ma, nb)
-A_mul_B!(c, 1, a, 1, ma, na, b, 1, nb) 
-@test c ≈ a*b
-mb = 2
-nb = 3
-b = randn(mb, nb)
-c = randn(na, nb)
-At_mul_B!(c, 1, a, 1, na, ma, b, 1, nb) 
-@test c ≈ transpose(a)*b
-nb = 4
-mb = 3
-b = randn(mb, nb)
-c = randn(ma, mb)
-A_mul_Bt!(c, 1, a, 1, ma, na, b, 1, mb) 
-@test c ≈ a*transpose(b)
-mb = 4
-nb = 2
-b = randn(mb, nb)
-c = randn(na, mb)
-At_mul_Bt!(c, 1, a, 1, na, ma, b, 1, mb) 
-@test c ≈ transpose(a)*transpose(b)
-
-order=2
-ma = 2
-na = 4
-q = 2
-p = na^(order-1)
-a = rand(ma, na)
-b_orig = rand(q*na^order)
-c = rand(ma*p*q)
-b = copy(b_orig)
-kron_mul_elem!(c, a, b, p, q)
-@test c ≈ kron(kron(I(p),a),I(q))*b_orig
-
+    @testset "Main kronecker operations" begin
+        order = 2
+        ma = 2
+        na = 4
+        q = 2
+        p = na^(order - 1)
+        a = rand(ma, na)
+        b_orig = rand(q * na^order)
+        c = rand(ma * p * q)
+        b = copy(b_orig)
+        kron_mul_elem!(c, a, b, p, q)
+        @test c ≈ kron(kron(I(p), a), I(q)) * b_orig
 
-ma = 2
-na = 2
-a = randn(ma,na)
-mc = 3
-order = 4
-c = randn(mc,mc)
-b = randn(na,mc^order)
-b_orig = copy(b)
-d = randn(ma,mc^order)
-w1 = Vector{Float64}(undef, ma*mc^order)
-w2 = Vector{Float64}(undef, ma*mc^order)
-a_mul_b_kron_c!(d, a, b, c, order, w1, w2)
-cc = c
-for i = 2:order
-    global cc = kron(cc,c)
-end
-@test d ≈ a*b_orig*cc
+        ma = 2
+        na = 2
+        a = randn(ma, na)
+        mc = 3
+        order = 4
+        c = randn(mc, mc)
+        b = randn(na, mc^order)
+        b_orig = copy(b)
+        d = randn(ma, mc^order)
+        w1 = Vector{Float64}(undef, ma * mc^order)
+        w2 = Vector{Float64}(undef, ma * mc^order)
+        a_mul_b_kron_c!(d, a, b, c, order, w1, w2)
+        cc = c
+        for i = 2:order
+            cc = kron(cc, c)
+        end
+        @test d ≈ a * b_orig * cc
 
-b = copy(b_orig)
-a_mul_kron_b!(d,b,c,order)
-cc = c
-for i = 2:order
-    global cc = kron(cc,c)
-end
-@test d ≈ b_orig*cc
+        b = copy(b_orig)
+        a_mul_kron_b!(d, b, c, order)
+        cc = c
+        for i = 2:order
+            cc = kron(cc, c)
+        end
+        @test d ≈ b_orig * cc
 
-order = 3
-ma = 2
-na = 4
-a = randn(ma,na)
-mb1 = 2
-nb1 = 4
-b1 = randn(mb1,nb1)
-mb2 = 2
-nb2 = 2
-b2 = randn(mb2,nb2)
-b = [b1, b2]
-c = randn(ma,nb1*nb2)
-work = zeros(16)
-a_mul_kron_b!(c,a,b,work)
-@test c ≈ a*kron(b[1],b[2])
+        order = 3
+        ma = 2
+        na = 4
+        a = randn(ma, na)
+        mb1 = 2
+        nb1 = 4
+        b1 = randn(mb1, nb1)
+        mb2 = 2
+        nb2 = 2
+        b2 = randn(mb2, nb2)
+        b = [b1, b2]
+        c = randn(ma, nb1 * nb2)
+        work = zeros(16)
+        a_mul_kron_b!(c, a, b, work)
+        @test c ≈ a * kron(b[1], b[2])
 
-order = 3
-ma = 2
-na = 4
-a = randn(ma,na)
-mb = 4
-nb = 8
-b = randn(mb,nb)
-c = randn(2,2)
-d = randn(2,2)
-work1 = zeros(mb*nb)
-work2 = zeros(mb*nb)
-e = randn(ma,8)
-a_mul_b_kron_c_d!(e, a, b, c, d, order, work1, work2)
-@test e ≈ a*b*kron(c,kron(d,d))
+        order = 3
+        ma = 2
+        na = 4
+        a = randn(ma, na)
+        mb = 4
+        nb = 8
+        b = randn(mb, nb)
+        c = randn(2, 2)
+        d = randn(2, 2)
+        work1 = zeros(mb * nb)
+        work2 = zeros(mb * nb)
+        e = randn(ma, 8)
+        a_mul_b_kron_c_d!(e, a, b, c, d, order, work1, work2)
+        @test e ≈ a * b * kron(c, kron(d, d))
 
-order = 2
-ma = 2
-na = 4
-a = randn(ma,na)
-mb = 4
-nb = 8
-b = randn(mb,nb)
-c = randn(ma*ma*nb)
-d = randn(na*na*mb)
-work1 = rand(na*na*mb)
-work2 = rand(na*na*mb)
-@time kron_at_kron_b_mul_c!(d, a, order, b, c, work1, work2)
-@time kron_at_kron_b_mul_c!(d, a, order, b, c, work1, work2)
+        order = 2
+        ma = 2
+        na = 4
+        a = randn(ma, na)
+        mb = 4
+        nb = 8
+        b = randn(mb, nb)
+        c = randn(ma * ma * nb)
+        d = randn(na * na * mb)
+        work1 = rand(na * na * mb)
+        work2 = rand(na * na * mb)
+        kron_at_kron_b_mul_c!(d, a, order, b, c, work1, work2)
+        kron_at_kron_b_mul_c!(d, a, order, b, c, work1, work2)
 
-@test d ≈ kron(kron(a',a'),b)*c
+        @test d ≈ kron(kron(a', a'), b) * c
+    end
 
-println("test1")
-order = 2
-ma = 2
-na = 4
-q = 2
-a = rand(ma, na)
-b = rand(q*na^order)
-c = rand(q*ma^order)
-work1 = rand(q*max(ma, na)^order)
-work2 = similar(work1)
-q=2
-kron_a_mul_b!(c, a, order, b, q, work1, work2)
-@time kron_a_mul_b!(c, a, order, b, q, work1, work2)
-@test c ≈ kron(kron(a,a),Matrix{Float64}(I(q)))*b
+    @testset "kron_a_mul_b" begin
+        order = 2
+        ma = 2
+        na = 4
+        q = 2
+        a = rand(ma, na)
+        b = rand(q * na^order)
+        c = rand(q * ma^order)
+        work1 = rand(q * max(ma, na)^order)
+        work2 = similar(work1)
+        q = 2
+        kron_a_mul_b!(c, a, order, b, q, work1, work2)
+        kron_a_mul_b!(c, a, order, b, q, work1, work2)
+        @test c ≈ kron(kron(a, a), Matrix{Float64}(I(q))) * b
 
-println("test2")
-b = rand(q*ma^order)
-c = rand(q*na^order)
-kron_at_mul_b!(c, a, order, b, q, work1, work2)
-@time kron_at_mul_b!(c, a, order, b, q, work1, work2)
-@test c ≈ kron(kron(a',a'),Matrix{Float64}(I(q)))*b
+        # test2
+        b = rand(q * ma^order)
+        c = rand(q * na^order)
+        kron_at_mul_b!(c, a, order, b, q, work1, work2)
+        kron_at_mul_b!(c, a, order, b, q, work1, work2)
+        @test c ≈ kron(kron(a', a'), Matrix{Float64}(I(q))) * b
+    end
 
-println("test3")
-order = 2
-mc = 2
-nc = 3
-c = rand(mc,nc)
-ma = 2
-na = 4
-a = rand(ma,na)
-nb = nc^order
-b = rand(na,nb)
-d = rand(ma,mc^order)
-work1 = rand(ma*max(mc, nc)^order)
-work2 = similar(work1)
-a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
-@time a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
-@test d ≈ a*b*kron(c',c')
+    @testset "a_mul_b_kron_c" begin
+        order = 2
+        mc = 2
+        nc = 3
+        c = rand(mc, nc)
+        ma = 2
+        na = 4
+        a = rand(ma, na)
+        nb = nc^order
+        b = rand(na, nb)
+        d = rand(ma, mc^order)
+        work1 = rand(ma * max(mc, nc)^order)
+        work2 = similar(work1)
+        a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
+        a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
+        @test d ≈ a * b * kron(c', c')
 
-println("test4")
-nb = mc^order
-b = rand(ma,nb)
-d = rand(na,nc^order)
-work1 = rand(na*max(mc, nc)^order)
-work2 = similar(work1)
-at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
-@time at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
-@test d ≈ a'*b*kron(c,c)
+        # test4
+        nb = mc^order
+        b = rand(ma, nb)
+        d = rand(na, nc^order)
+        work1 = rand(na * max(mc, nc)^order)
+        work2 = similar(work1)
+        at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
+        at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
+        @test d ≈ a' * b * kron(c, c)
 
-println("test1a")
-order = 2
-ma = 4
-na = 4
-q = 2
-a = QuasiUpperTriangular(triu(rand(ma, na)))
-b = rand(q*na^order)
-c = rand(q*ma^order)
-work1 = rand(q*max(ma, na)^order)
-work2 = similar(work1)
-kron_a_mul_b!(c, a, order, b, q, work1, work2)
-@time kron_a_mul_b!(c, a, order, b, q, work1, work2)
-@test c ≈ kron(kron(a,a),Matrix{Float64}(I(q)))*b
+        # Doubtful usefulness, probably matching different combination of sizes
+        order = 2
+        mc = 2
+        nc = 3
+        c = rand(mc, nc)
+        ma = 4
+        na = 4
+        a = rand(ma, na)
+        nb = nc^order
+        b = rand(na, nb)
+        d = rand(ma, mc^order)
+        work1 = rand(ma * max(mc, nc)^order)
+        work2 = similar(work1)
+        a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
+        a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
+        @test d ≈ a * b * kron(c', c')
 
-println("test2a")
-b = rand(q*ma^order)
-c = rand(q*na^order)
-kron_at_mul_b!(c, a, order, b, q, work1, work2)
-@time kron_at_mul_b!(c, a, order, b, q, work1, work2)
-@test c ≈ kron(kron(a',a'),Matrix{Float64}(I(q)))*b
+        #test4a
+        nb = mc^order
+        b = rand(ma, nb)
+        d = rand(na, nc^order)
+        work1 = rand(na * max(mc, nc)^order)
+        work2 = similar(work1)
+        at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
+        at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
+        @test d ≈ a' * b * kron(c, c)
+    end
 
-println("test3a")
-order = 2
-mc = 2
-nc = 3
-c = rand(mc,nc)
-ma = 4
-na = 4
-a = rand(ma,na)
-nb = nc^order
-b = rand(na,nb)
-d = rand(ma,mc^order)
-work1 = rand(ma*max(mc, nc)^order)
-work2 = similar(work1)
-a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
-@time a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
-@test d ≈ a*b*kron(c',c')
+    @testset "QuasiUpperTriangular" begin
+        order = 2
+        ma = 4
+        na = 4
+        q = 2
+        a = QuasiUpperTriangular(triu(rand(ma, na)))
+        b = rand(q * na^order)
+        c = rand(q * ma^order)
+        work1 = rand(q * max(ma, na)^order)
+        work2 = similar(work1)
+        kron_a_mul_b!(c, a, order, b, q, work1, work2)
+        kron_a_mul_b!(c, a, order, b, q, work1, work2)
+        @test c ≈ kron(kron(a, a), Matrix{Float64}(I(q))) * b
 
-println("test4a")
-nb = mc^order
-b = rand(ma,nb)
-d = rand(na,nc^order)
-work1 = rand(na*max(mc, nc)^order)
-work2 = similar(work1)
-at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
-@time at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
-@test d ≈ a'*b*kron(c,c)
+        # test2a
+        b = rand(q * ma^order)
+        c = rand(q * na^order)
+        kron_at_mul_b!(c, a, order, b, q, work1, work2)
+        kron_at_mul_b!(c, a, order, b, q, work1, work2)
+        @test c ≈ kron(kron(a', a'), Matrix{Float64}(I(q))) * b
 
-println("test1b")
-order = 2
-ma = 4
-na = 4
-q = 2
-a = QuasiUpperTriangular(triu(rand(ma, na)))
-b = view(rand(q*na^order),:)
-c = view(rand(q*ma^order),:)
-work1 = view(rand(q*max(ma, na)^order),:)
-work2 = view(rand(q*max(ma, na)^order),:)
-kron_a_mul_b!(c, a, order, b, q, work1, work2)
-@time kron_a_mul_b!(c, a, order, b, q, work1, work2)
-@test c ≈ kron(kron(a,a),Matrix{Float64}(I(q)))*b
+        # Testing the same thing with views
+        order = 2
+        ma = 4
+        na = 4
+        q = 2
+        a = QuasiUpperTriangular(triu(rand(ma, na)))
+        b = view(rand(q * na^order), :)
+        c = view(rand(q * ma^order), :)
+        work1 = view(rand(q * max(ma, na)^order), :)
+        work2 = view(rand(q * max(ma, na)^order), :)
+        kron_a_mul_b!(c, a, order, b, q, work1, work2)
+        kron_a_mul_b!(c, a, order, b, q, work1, work2)
+        @test c ≈ kron(kron(a, a), Matrix{Float64}(I(q))) * b
 
-println("test2b")
-b = view(rand(q*ma^order),:)
-c = view(rand(q*na^order),:)
-work1 = view(rand(q*max(ma, na)^order),:)
-work2 = view(rand(q*max(ma, na)^order),:)
-kron_at_mul_b!(c, a, order, b, q, work1, work2)
-@time kron_at_mul_b!(c, a, order, b, q, work1, work2)
-@test c ≈ kron(kron(a',a'),Matrix{Float64}(I(q)))*b
+        # test2b
+        b = view(rand(q * ma^order), :)
+        c = view(rand(q * na^order), :)
+        work1 = view(rand(q * max(ma, na)^order), :)
+        work2 = view(rand(q * max(ma, na)^order), :)
+        kron_at_mul_b!(c, a, order, b, q, work1, work2)
+        kron_at_mul_b!(c, a, order, b, q, work1, work2)
+        @test c ≈ kron(kron(a', a'), Matrix{Float64}(I(q))) * b
+    end
 
-println("test3b")
-order = 2
-mc = 2
-nc = 3
-c = view(rand(mc,nc),:,:)
-ma = 4
-na = 4
-a = rand(ma,na)
-nb = nc^order
-b = view(rand(na,nb),:,:)
-d = view(rand(ma,mc^order),:,:)
-work1 = view(rand(ma*max(mc, nc)^order),:)
-work2 = view(rand(ma*max(mc, nc)^order),:)
-a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
-@time a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
-@test d ≈ a*b*kron(c',c')
+    @testset "Views" begin
+        order = 2
+        mc = 2
+        nc = 3
+        c = view(rand(mc, nc), :, :)
+        ma = 4
+        na = 4
+        a = rand(ma, na)
+        nb = nc^order
+        b = view(rand(na, nb), :, :)
+        d = view(rand(ma, mc^order), :, :)
+        work1 = view(rand(ma * max(mc, nc)^order), :)
+        work2 = view(rand(ma * max(mc, nc)^order), :)
+        a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
+        a_mul_b_kron_ct!(d, a, b, c, order, work1, work2)
+        @test d ≈ a * b * kron(c', c')
 
-println("test4b")
-nb = mc^order
-b = view(rand(ma,nb),:,:)
-d = view(rand(na,nc^order),:,:)
-work1 = rand(na*max(mc, nc)^order)
-work2 = similar(work1)
-at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
-@time at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
-@test d ≈ a'*b*kron(c,c)
+        # test3b
+        nb = mc^order
+        b = view(rand(ma, nb), :, :)
+        d = view(rand(na, nc^order), :, :)
+        work1 = rand(na * max(mc, nc)^order)
+        work2 = similar(work1)
+        at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
+        at_mul_b_kron_c!(d, a, b, c, order, work1, work2)
+        @test d ≈ a' * b * kron(c, c)
+    end
+end