diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl
index 8ab3b3ec6..83d53fb0b 100644
--- a/src/factorizations/matrixalgebrakit.jl
+++ b/src/factorizations/matrixalgebrakit.jl
@@ -64,6 +64,38 @@ end
 
 MAK.zero!(t::AbstractTensorMap) = zerovector!(t)
 
+# Default algorithm
+# -----------------
+for f in [
+        :lq_full, :lq_compact, :lq_null,
+        :qr_full, :qr_compact, :qr_null,
+        :schur_full, :schur_vals,
+        :eig_full, :eig_vals, :eig_trunc, :eig_trunc_no_error,
+        :eigh_full, :eigh_vals, :eigh_trunc, :eigh_trunc_no_error,
+        :svd_full, :svd_compact, :svd_trunc, :svd_trunc_no_error, :svd_vals,
+        :left_polar, :right_polar,
+        :left_orth, :right_orth, :left_null, :right_null,
+        :project_hermitian, :project_antihermitian, :project_isometric,
+    ]
+    f! = Symbol(f, :!)
+    @eval MAK.$f!(t::AbstractTensorMap, alg::DefaultAlgorithm) =
+        MAK.$f!(t, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
+    @eval MAK.$f!(t::AbstractTensorMap, out, alg::DefaultAlgorithm) =
+        MAK.$f!(t, out, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
+
+    # disambiguate
+    @eval MAK.$f!(t::AdjointTensorMap, alg::DefaultAlgorithm) =
+        MAK.$f!(t, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
+    @eval MAK.$f!(t::AdjointTensorMap, out, alg::DefaultAlgorithm) =
+        MAK.$f!(t, out, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
+
+    @eval MAK.$f!(t::DiagonalTensorMap, alg::DefaultAlgorithm) =
+        MAK.$f!(t, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
+    @eval MAK.$f!(t::DiagonalTensorMap, out, alg::DefaultAlgorithm) =
+        MAK.$f!(t, out, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
+end
+
+
 # Singular value decomposition
 # ----------------------------
 function MAK.initialize_output(::typeof(svd_full!), t::AbstractTensorMap, ::AbstractAlgorithm)
diff --git a/test/cuda/factorizations.jl b/test/cuda/factorizations.jl
index fdeca843d..bdb141235 100644
--- a/test/cuda/factorizations.jl
+++ b/test/cuda/factorizations.jl
@@ -2,7 +2,7 @@ using Adapt, CUDA, CUDA.cuRAND, cuTENSOR
 using Test, TestExtras
 using TensorKit
 using LinearAlgebra: LinearAlgebra
-using MatrixAlgebraKit: diagview
+using MatrixAlgebraKit: DefaultAlgorithm, diagview
 const CUDAExt = Base.get_extension(TensorKit, :TensorKitCUDAExt)
 @assert !isnothing(CUDAExt) "Failed to load TensorKit - CUDA extension"
 const CuTensorMap = getglobal(CUDAExt, :CuTensorMap)
@@ -35,21 +35,41 @@ for V in spacelist
                 @test Q * R ≈ t
                 @test isunitary(Q)
 
+                Q, R = @constinferred qr_full(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isunitary(Q)
+
                 Q, R = @constinferred qr_compact(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
 
+                Q, R = @constinferred qr_compact(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+
                 Q, R = @constinferred left_orth(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
 
+                Q, R = @constinferred left_orth(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+
                 N = @constinferred qr_null(t)
                 @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
 
+                N = @constinferred qr_null(t, DefaultAlgorithm())
+                @test isisometric(N)
+                @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
+
                 N = @constinferred left_null(t)
                 @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
+
+                N = @constinferred left_null(t, DefaultAlgorithm())
+                @test isisometric(N)
+                @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
             end
 
             # empty tensor
@@ -61,19 +81,38 @@ for V in spacelist
                 @test isunitary(Q)
                 @test dim(R) == dim(t) == 0
 
+                Q, R = @constinferred qr_full(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isunitary(Q)
+                @test dim(R) == dim(t) == 0
+
                 Q, R = @constinferred qr_compact(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
                 @test dim(Q) == dim(R) == dim(t)
 
+                Q, R = @constinferred qr_compact(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+                @test dim(Q) == dim(R) == dim(t)
+
                 Q, R = @constinferred left_orth(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
                 @test dim(Q) == dim(R) == dim(t)
 
+                Q, R = @constinferred left_orth(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+                @test dim(Q) == dim(R) == dim(t)
+
                 N = @constinferred qr_null(t)
                 @test isunitary(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
+
+                N = @constinferred qr_null(t, DefaultAlgorithm())
+                @test isunitary(N)
+                @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
             end
         end
 
@@ -90,17 +129,33 @@ for V in spacelist
                 @test L * Q ≈ t
                 @test isunitary(Q)
 
+                L, Q = @constinferred lq_full(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isunitary(Q)
+
                 L, Q = @constinferred lq_compact(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
 
+                L, Q = @constinferred lq_compact(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+
                 L, Q = @constinferred right_orth(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
 
+                L, Q = @constinferred right_orth(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+
                 Nᴴ = @constinferred lq_null(t)
                 @test isisometric(Nᴴ; side = :right)
                 @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
+
+                Nᴴ = @constinferred lq_null(t, DefaultAlgorithm())
+                @test isisometric(Nᴴ; side = :right)
+                @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
             end
 
             for T in eltypes
@@ -112,19 +167,38 @@ for V in spacelist
                 @test isunitary(Q)
                 @test dim(L) == dim(t) == 0
 
+                L, Q = @constinferred lq_full(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isunitary(Q)
+                @test dim(L) == dim(t) == 0
+
                 L, Q = @constinferred lq_compact(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
                 @test dim(Q) == dim(L) == dim(t)
 
+                L, Q = @constinferred lq_compact(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+                @test dim(Q) == dim(L) == dim(t)
+
                 L, Q = @constinferred right_orth(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
                 @test dim(Q) == dim(L) == dim(t)
 
+                L, Q = @constinferred right_orth(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+                @test dim(Q) == dim(L) == dim(t)
+
                 Nᴴ = @constinferred lq_null(t)
                 @test isunitary(Nᴴ)
                 @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
+
+                Nᴴ = @constinferred lq_null(t, DefaultAlgorithm())
+                @test isunitary(Nᴴ)
+                @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
             end
         end
 
@@ -143,6 +217,11 @@ for V in spacelist
                 @test isisometric(w)
                 @test isposdef(p)
 
+                w, p = @constinferred left_polar(t, DefaultAlgorithm())
+                @test w * p ≈ t
+                @test isisometric(w)
+                @test isposdef(p)
+
                 w, p = @constinferred left_orth(t; alg = :polar)
                 @test w * p ≈ t
                 @test isisometric(w)
@@ -162,6 +241,11 @@ for V in spacelist
                 @test isisometric(wᴴ; side = :right)
                 @test isposdef(p)
 
+                p, wᴴ = @constinferred right_polar(t, DefaultAlgorithm())
+                @test p * wᴴ ≈ t
+                @test isisometric(wᴴ; side = :right)
+                @test isposdef(p)
+
                 p, wᴴ = @constinferred right_orth(t; alg = :polar)
                 @test p * wᴴ ≈ t
                 @test isisometric(wᴴ; side = :right)
@@ -182,16 +266,31 @@ for V in spacelist
                 @test isunitary(u)
                 @test isunitary(vᴴ)
 
+                u, s, vᴴ = @constinferred svd_full(t, DefaultAlgorithm())
+                @test u * s * vᴴ ≈ t
+                @test isunitary(u)
+                @test isunitary(vᴴ)
+
                 u, s, vᴴ = @constinferred svd_compact(t)
                 @test u * s * vᴴ ≈ t
                 @test isisometric(u)
                 @test isposdef(s)
                 @test isisometric(vᴴ; side = :right)
 
+                u, s, vᴴ = @constinferred svd_compact(t, DefaultAlgorithm())
+                @test u * s * vᴴ ≈ t
+                @test isisometric(u)
+                @test isposdef(s)
+                @test isisometric(vᴴ; side = :right)
+
                 s′ = @constinferred svd_vals(t)
                 @test parent(s′) ≈ parent(diagview(s))
                 @test s′ isa TensorKit.SectorVector
 
+                s′ = @constinferred svd_vals(t, DefaultAlgorithm())
+                @test parent(s′) ≈ parent(diagview(s))
+                @test s′ isa TensorKit.SectorVector
+
                 s2 = @constinferred DiagonalTensorMap(s′)
                 @test s2 ≈ s
 
@@ -230,9 +329,18 @@ for V in spacelist
                 @test isunitary(U)
                 @test isunitary(Vᴴ)
 
+                U, S, Vᴴ = @constinferred svd_full(t, DefaultAlgorithm())
+                @test U * S * Vᴴ ≈ t
+                @test isunitary(U)
+                @test isunitary(Vᴴ)
+
                 U, S, Vᴴ = @constinferred svd_compact(t)
                 @test U * S * Vᴴ ≈ t
                 @test dim(U) == dim(S) == dim(Vᴴ) == dim(t) == 0
+
+                U, S, Vᴴ = @constinferred svd_compact(t, DefaultAlgorithm())
+                @test U * S * Vᴴ ≈ t
+                @test dim(U) == dim(S) == dim(Vᴴ) == dim(t) == 0
             end
         end
 
@@ -253,6 +361,12 @@ for V in spacelist
                 @test isisometric(U)
                 @test isisometric(Vᴴ; side = :right)
 
+                U, S, Vᴴ, ϵ = @constinferred svd_trunc(t, DefaultAlgorithm(; trunc = notrunc()))
+                @test U * S * Vᴴ ≈ t
+                @test ϵ ≈ 0
+                @test isisometric(U)
+                @test isisometric(Vᴴ; side = :right)
+
                 # dimension of S is a float for IsingBimodule
                 nvals = round(Int, dim(domain(S)) / 2)
                 trunc = truncrank(nvals)
@@ -316,10 +430,17 @@ for V in spacelist
                 d, v = @constinferred eig_full(t)
                 @test t * v ≈ v * d
 
+                d, v = @constinferred eig_full(t, DefaultAlgorithm())
+                @test t * v ≈ v * d
+
                 d′ = @constinferred eig_vals(t)
                 @test parent(d′) ≈ parent(diagview(d))
                 @test d′ isa TensorKit.SectorVector
 
+                d′ = @constinferred eig_vals(t, DefaultAlgorithm())
+                @test parent(d′) ≈ parent(diagview(d))
+                @test d′ isa TensorKit.SectorVector
+
                 d2 = @constinferred DiagonalTensorMap(d′)
                 @test d2 ≈ d
 
@@ -332,12 +453,21 @@ for V in spacelist
                 @test t * v ≈ v * d
                 #test_dim_isapprox(domain(d), nvals)
 
+                d, v = @constinferred eig_trunc(t, DefaultAlgorithm(; trunc = truncrank(nvals)))
+                @test t * v ≈ v * d
+                #test_dim_isapprox(domain(d), nvals)
+
                 t2 = @constinferred project_hermitian(t)
                 D, V = eigen(t2)
                 @test isisometric(V)
                 D̃, Ṽ = @constinferred eigh_full(t2)
                 @test D ≈ D̃
                 @test V ≈ Ṽ
+
+                D̃, Ṽ = @constinferred eigh_full(t2, DefaultAlgorithm())
+                @test D ≈ D̃
+                @test V ≈ Ṽ
+
                 λ = minimum(real, parent(diagview(D)))
                 @test cond(Ṽ) ≈ one(real(T))
                 @test isposdef(t2) == isposdef(λ)
@@ -352,6 +482,10 @@ for V in spacelist
                 @test parent(d′) ≈ parent(diagview(d))
                 @test d′ isa TensorKit.SectorVector
 
+                d′ = @constinferred eigh_vals(t2, DefaultAlgorithm())
+                @test parent(d′) ≈ parent(diagview(d))
+                @test d′ isa TensorKit.SectorVector
+
                 λ = minimum(real, parent(diagview(d)))
                 @test cond(v) ≈ one(real(T))
                 @test isposdef(t2) == isposdef(λ)
@@ -361,6 +495,10 @@ for V in spacelist
                 d, v = @constinferred eigh_trunc(t2; trunc = truncrank(nvals))
                 @test t2 * v ≈ v * d
                 #test_dim_isapprox(domain(d), nvals)
+
+                d, v = @constinferred eigh_trunc(t2, DefaultAlgorithm(; trunc = truncrank(nvals)))
+                @test t2 * v ≈ v * d
+                #test_dim_isapprox(domain(d), nvals)
             end
         end
 
@@ -423,6 +561,11 @@ for V in spacelist
                 th′ = @constinferred project_hermitian(t)
                 @test ishermitian(th′)
                 @test th′ ≈ th
+
+                th′ = @constinferred project_hermitian(t, DefaultAlgorithm())
+                @test ishermitian(th′)
+                @test th′ ≈ th
+
                 @test t == tc
                 th_approx = th + noisefactor * ta
                 @test !ishermitian(th_approx) || (T <: Real && t isa DiagonalTensorMap)
@@ -431,6 +574,11 @@ for V in spacelist
                 ta′ = project_antihermitian(t)
                 @test isantihermitian(ta′)
                 @test ta′ ≈ ta
+
+                ta′ = @constinferred project_antihermitian(t, DefaultAlgorithm())
+                @test isantihermitian(ta′)
+                @test ta′ ≈ ta
+
                 @test t == tc
                 ta_approx = ta + noisefactor * th
                 @test !isantihermitian(ta_approx)
@@ -448,6 +596,10 @@ for V in spacelist
                     )
                 t2 = project_isometric(t)
                 @test isisometric(t2)
+                t2′ = @constinferred project_isometric(t, DefaultAlgorithm())
+                @test isisometric(t2′)
+                @test t2′ * ((t2′)' * t) ≈ t
+
                 t3 = project_isometric(t2)
                 @test t3 ≈ t2 # stability of the projection
                 @test t2 * (t2' * t) ≈ t
diff --git a/test/factorizations/eig.jl b/test/factorizations/eig.jl
index 84909a030..55ff622cd 100644
--- a/test/factorizations/eig.jl
+++ b/test/factorizations/eig.jl
@@ -1,7 +1,7 @@
 using Test, TestExtras
 using TensorKit
 using LinearAlgebra: LinearAlgebra
-using MatrixAlgebraKit: diagview
+using MatrixAlgebraKit: DefaultAlgorithm, diagview
 
 
 spacelist = factorization_spacelist(fast_tests)
@@ -29,10 +29,17 @@ for V in spacelist
                 d, v = @constinferred eig_full(t)
                 @test t * v ≈ v * d
 
+                d, v = @constinferred eig_full(t, DefaultAlgorithm())
+                @test t * v ≈ v * d
+
                 d′ = @constinferred eig_vals(t)
                 @test d′ ≈ diagview(d)
                 @test d′ isa TensorKit.SectorVector
 
+                d′ = @constinferred eig_vals(t, DefaultAlgorithm())
+                @test d′ ≈ diagview(d)
+                @test d′ isa TensorKit.SectorVector
+
                 d2 = @constinferred DiagonalTensorMap(d′)
                 @test d2 ≈ d
 
@@ -45,12 +52,21 @@ for V in spacelist
                 @test t * v ≈ v * d
                 @test abs(dim(domain(d)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
 
+                d, v = @constinferred eig_trunc(t, DefaultAlgorithm(; trunc = truncrank(nvals)))
+                @test t * v ≈ v * d
+                @test abs(dim(domain(d)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
+
                 t2 = @constinferred project_hermitian(t)
                 D, V = eigen(t2)
                 @test isisometric(V)
                 D̃, Ṽ = @constinferred eigh_full(t2)
                 @test D ≈ D̃
                 @test V ≈ Ṽ
+
+                D̃, Ṽ = @constinferred eigh_full(t2, DefaultAlgorithm())
+                @test D ≈ D̃
+                @test V ≈ Ṽ
+
                 λ = minimum(real, diagview(D))
                 @test cond(Ṽ) ≈ one(real(T))
                 @test isposdef(t2) == isposdef(λ)
@@ -65,6 +81,10 @@ for V in spacelist
                 @test d′ ≈ diagview(d)
                 @test d′ isa TensorKit.SectorVector
 
+                d′ = @constinferred eigh_vals(t2, DefaultAlgorithm())
+                @test d′ ≈ diagview(d)
+                @test d′ isa TensorKit.SectorVector
+
                 λ = minimum(real, diagview(d))
                 @test cond(v) ≈ one(real(T))
                 @test isposdef(t2) == isposdef(λ)
@@ -74,6 +94,10 @@ for V in spacelist
                 d, v = @constinferred eigh_trunc(t2; trunc = truncrank(nvals))
                 @test t2 * v ≈ v * d
                 @test abs(dim(domain(d)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
+
+                d, v = @constinferred eigh_trunc(t2, DefaultAlgorithm(; trunc = truncrank(nvals)))
+                @test t2 * v ≈ v * d
+                @test abs(dim(domain(d)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1)
             end
         end
     end
diff --git a/test/factorizations/ortho.jl b/test/factorizations/ortho.jl
index 768f9f33e..170d596a2 100644
--- a/test/factorizations/ortho.jl
+++ b/test/factorizations/ortho.jl
@@ -1,7 +1,7 @@
 using Test, TestExtras
 using TensorKit
 using LinearAlgebra: LinearAlgebra
-using MatrixAlgebraKit: diagview
+using MatrixAlgebraKit: DefaultAlgorithm, diagview
 
 
 spacelist = factorization_spacelist(fast_tests)
@@ -32,21 +32,41 @@ for V in spacelist
                 @test Q * R ≈ t
                 @test isunitary(Q)
 
+                Q, R = @constinferred qr_full(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isunitary(Q)
+
                 Q, R = @constinferred qr_compact(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
 
+                Q, R = @constinferred qr_compact(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+
                 Q, R = @constinferred left_orth(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
 
+                Q, R = @constinferred left_orth(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+
                 N = @constinferred qr_null(t)
                 @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
 
+                N = @constinferred qr_null(t, DefaultAlgorithm())
+                @test isisometric(N)
+                @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
+
                 N = @constinferred left_null(t)
                 @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
+
+                N = @constinferred left_null(t, DefaultAlgorithm())
+                @test isisometric(N)
+                @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
             end
 
             # empty tensor
@@ -58,19 +78,38 @@ for V in spacelist
                 @test isunitary(Q)
                 @test dim(R) == dim(t) == 0
 
+                Q, R = @constinferred qr_full(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isunitary(Q)
+                @test dim(R) == dim(t) == 0
+
                 Q, R = @constinferred qr_compact(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
                 @test dim(Q) == dim(R) == dim(t)
 
+                Q, R = @constinferred qr_compact(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+                @test dim(Q) == dim(R) == dim(t)
+
                 Q, R = @constinferred left_orth(t)
                 @test Q * R ≈ t
                 @test isisometric(Q)
                 @test dim(Q) == dim(R) == dim(t)
 
+                Q, R = @constinferred left_orth(t, DefaultAlgorithm())
+                @test Q * R ≈ t
+                @test isisometric(Q)
+                @test dim(Q) == dim(R) == dim(t)
+
                 N = @constinferred qr_null(t)
                 @test isunitary(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
+
+                N = @constinferred qr_null(t, DefaultAlgorithm())
+                @test isunitary(N)
+                @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
             end
         end
 
@@ -87,17 +126,33 @@ for V in spacelist
                 @test L * Q ≈ t
                 @test isunitary(Q)
 
+                L, Q = @constinferred lq_full(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isunitary(Q)
+
                 L, Q = @constinferred lq_compact(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
 
+                L, Q = @constinferred lq_compact(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+
                 L, Q = @constinferred right_orth(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
 
+                L, Q = @constinferred right_orth(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+
                 Nᴴ = @constinferred lq_null(t)
                 @test isisometric(Nᴴ; side = :right)
                 @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
+
+                Nᴴ = @constinferred lq_null(t, DefaultAlgorithm())
+                @test isisometric(Nᴴ; side = :right)
+                @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
             end
 
             for T in eltypes
@@ -109,19 +164,38 @@ for V in spacelist
                 @test isunitary(Q)
                 @test dim(L) == dim(t) == 0
 
+                L, Q = @constinferred lq_full(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isunitary(Q)
+                @test dim(L) == dim(t) == 0
+
                 L, Q = @constinferred lq_compact(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
                 @test dim(Q) == dim(L) == dim(t)
 
+                L, Q = @constinferred lq_compact(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+                @test dim(Q) == dim(L) == dim(t)
+
                 L, Q = @constinferred right_orth(t)
                 @test L * Q ≈ t
                 @test isisometric(Q; side = :right)
                 @test dim(Q) == dim(L) == dim(t)
 
+                L, Q = @constinferred right_orth(t, DefaultAlgorithm())
+                @test L * Q ≈ t
+                @test isisometric(Q; side = :right)
+                @test dim(Q) == dim(L) == dim(t)
+
                 Nᴴ = @constinferred lq_null(t)
                 @test isunitary(Nᴴ)
                 @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
+
+                Nᴴ = @constinferred lq_null(t, DefaultAlgorithm())
+                @test isunitary(Nᴴ)
+                @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
             end
         end
     end
diff --git a/test/factorizations/projections.jl b/test/factorizations/projections.jl
index 209bf4ef6..fd32a804b 100644
--- a/test/factorizations/projections.jl
+++ b/test/factorizations/projections.jl
@@ -1,7 +1,7 @@
 using Test, TestExtras
 using TensorKit
 using LinearAlgebra: LinearAlgebra
-using MatrixAlgebraKit: diagview
+using MatrixAlgebraKit: DefaultAlgorithm, diagview
 
 
 spacelist = factorization_spacelist(fast_tests)
@@ -35,6 +35,11 @@ for V in spacelist
                 th′ = @constinferred project_hermitian(t)
                 @test ishermitian(th′)
                 @test th′ ≈ th
+
+                th′ = @constinferred project_hermitian(t, DefaultAlgorithm())
+                @test ishermitian(th′)
+                @test th′ ≈ th
+
                 @test t == tc
                 th_approx = th + noisefactor * ta
                 @test !ishermitian(th_approx) || (T <: Real && t isa DiagonalTensorMap)
@@ -43,6 +48,11 @@ for V in spacelist
                 ta′ = project_antihermitian(t)
                 @test isantihermitian(ta′)
                 @test ta′ ≈ ta
+
+                ta′ = @constinferred project_antihermitian(t, DefaultAlgorithm())
+                @test isantihermitian(ta′)
+                @test ta′ ≈ ta
+
                 @test t == tc
                 ta_approx = ta + noisefactor * th
                 @test !isantihermitian(ta_approx)
@@ -64,6 +74,10 @@ for V in spacelist
                     )
                 t2 = project_isometric(t)
                 @test isisometric(t2)
+                t2′ = @constinferred project_isometric(t, DefaultAlgorithm())
+                @test isisometric(t2′)
+                @test t2′ * ((t2′)' * t) ≈ t
+
                 t3 = project_isometric(t2)
                 @test t3 ≈ t2 # stability of the projection
                 @test t2 * (t2' * t) ≈ t
diff --git a/test/factorizations/svd.jl b/test/factorizations/svd.jl
index cec5e1dda..0678db827 100644
--- a/test/factorizations/svd.jl
+++ b/test/factorizations/svd.jl
@@ -1,7 +1,7 @@
 using Test, TestExtras
 using TensorKit
 using LinearAlgebra: LinearAlgebra
-using MatrixAlgebraKit: defaulttol, diagview
+using MatrixAlgebraKit: DefaultAlgorithm, defaulttol, diagview
 
 spacelist = factorization_spacelist(fast_tests)
 
@@ -71,6 +71,10 @@ for V in spacelist
                 @test isisometric(w)
                 @test isposdef(p)
 
+                w′, p′ = @constinferred left_polar(t, DefaultAlgorithm())
+                @test w ≈ w′
+                @test p ≈ p′
+
                 w, p = @constinferred left_orth(t; alg = :polar)
                 @test w * p ≈ t
                 @test isisometric(w)
@@ -90,6 +94,10 @@ for V in spacelist
                 @test isisometric(wᴴ; side = :right)
                 @test isposdef(p)
 
+                p′, wᴴ′ = @constinferred right_polar(t, DefaultAlgorithm())
+                @test p′ ≈ p
+                @test wᴴ′ ≈ wᴴ
+
                 p, wᴴ = @constinferred right_orth(t; alg = :polar)
                 @test p * wᴴ ≈ t
                 @test isisometric(wᴴ; side = :right)
@@ -110,16 +118,30 @@ for V in spacelist
                 @test isunitary(u)
                 @test isunitary(vᴴ)
 
+                u′, s′, vᴴ′ = @constinferred svd_full(t, DefaultAlgorithm())
+                @test u ≈ u′
+                @test s ≈ s′
+                @test vᴴ ≈ vᴴ′
+
                 u, s, vᴴ = @constinferred svd_compact(t)
                 @test u * s * vᴴ ≈ t
                 @test isisometric(u)
                 @test isposdef(s)
                 @test isisometric(vᴴ; side = :right)
 
+                u′, s′, vᴴ′ = @constinferred svd_compact(t, DefaultAlgorithm())
+                @test u ≈ u′
+                @test s ≈ s′
+                @test vᴴ ≈ vᴴ′
+
                 s′ = @constinferred svd_vals(t)
                 @test s′ ≈ diagview(s)
                 @test s′ isa TensorKit.SectorVector
 
+                s′ = @constinferred svd_vals(t, DefaultAlgorithm())
+                @test s′ ≈ diagview(s)
+                @test s′ isa TensorKit.SectorVector
+
                 s2 = @constinferred DiagonalTensorMap(s′)
                 @test s2 ≈ s
 
@@ -157,9 +179,18 @@ for V in spacelist
                 @test isunitary(U)
                 @test isunitary(Vᴴ)
 
+                U, S, Vᴴ = @constinferred svd_full(t, DefaultAlgorithm())
+                @test U * S * Vᴴ ≈ t
+                @test isunitary(U)
+                @test isunitary(Vᴴ)
+
                 U, S, Vᴴ = @constinferred svd_compact(t)
                 @test U * S * Vᴴ ≈ t
                 @test dim(U) == dim(S) == dim(Vᴴ) == dim(t) == 0
+
+                U, S, Vᴴ = @constinferred svd_compact(t, DefaultAlgorithm())
+                @test U * S * Vᴴ ≈ t
+                @test dim(U) == dim(S) == dim(Vᴴ) == dim(t) == 0
             end
         end
 
@@ -180,6 +211,12 @@ for V in spacelist
                 @test isisometric(U)
                 @test isisometric(Vᴴ; side = :right)
 
+                U, S, Vᴴ, ϵ = @constinferred svd_trunc(t, DefaultAlgorithm(; trunc = notrunc()))
+                @test U * S * Vᴴ ≈ t
+                @test ϵ ≈ 0
+                @test isisometric(U)
+                @test isisometric(Vᴴ; side = :right)
+
                 # when rank of t is already smaller than truncrank
                 t_rank = ceil(Int, min(dim(codomain(t)), dim(domain(t))))
                 U, S, Vᴴ, ϵ = @constinferred svd_trunc(t; trunc = truncrank(t_rank + 1))
@@ -188,6 +225,12 @@ for V in spacelist
                 @test isisometric(U)
                 @test isisometric(Vᴴ; side = :right)
 
+                U, S, Vᴴ, ϵ = @constinferred svd_trunc(t, DefaultAlgorithm(; trunc = truncrank(t_rank + 1)))
+                @test U * S * Vᴴ ≈ t
+                @test ϵ ≈ 0
+                @test isisometric(U)
+                @test isisometric(Vᴴ; side = :right)
+
                 # dimension of S is a float for IsingBimodule
                 nvals = round(Int, dim(domain(S)) / 2)
                 trunc = truncrank(nvals)