diff --git a/src/NCMAlgorithm.jl b/src/NCMAlgorithm.jl
index 1ef7063..38ebf26 100644
--- a/src/NCMAlgorithm.jl
+++ b/src/NCMAlgorithm.jl
@@ -52,14 +52,15 @@ default_tol(::Type{<:Rational}) = 0
 default_tol(::Type{<:Integer}) = 0
 
 """
-    supports_mask(alg)
+    supports_mask(algtype)
 
 Trait for whether an algorithm supports a fixed-element mask. When `true`, the mask is enforced
 during the solve (for alternating projections, this means projecting onto the fixed-element
 subspace each iteration). When `false`, passing a mask throws an informative error. Default is
 `false`.
 """
-supports_mask(::NCMAlgorithm) = false
+supports_mask(::Type{T}) where {T <: NCMAlgorithm} = false
+supports_mask(alg::T) where {T <: NCMAlgorithm} = supports_mask(T)
 
 """
     default_iters(alg, A)
@@ -100,7 +101,7 @@ If `false`, then a copy using the upper or lower matrix is used instead.
 supports_symmetric(::NCMAlgorithm) = false
 
 """
-    supports_parameterless_construction(alg)
+    supports_parameterless_construction(algtype)
 
 Trait for if an algorithm can be constructed without any parameters (default is `false`).
 """
diff --git a/src/NCMSolution.jl b/src/NCMSolution.jl
index 87af951..0077752 100644
--- a/src/NCMSolution.jl
+++ b/src/NCMSolution.jl
@@ -83,16 +83,17 @@ function CommonSolve.solve!(solver::NCMSolver, args...; kwargs...)
     sol = solve!(solver, solver.alg, args...; kwargs...)
 
     if sol.solver.ensure_pd && !isposdef(sol.X)
+        project_psd!(sol.X, solver.min_eigenvalue)
+
         # Strict PD and exact fixed-element feasibility cannot both be guaranteed: repairing
         # definiteness perturbs every entry, so re-apply the mask afterwards. The fixed elements
         # (and unit diagonal) take precedence - the result is PD up to O(√eps).
-        project_psd!(sol.X, sqrt(eps(eltype(sol.X))))
-
         if sol.solver.mask !== nothing
             project_fixed!(sol.X, sol.solver.A_orig, sol.solver.mask)
+            project_unit!(sol.X)
+        else
+            cov2cor!(sol.X)
         end
-
-        project_unit!(sol.X)
     end
 
     return sol
diff --git a/src/NCMSolver.jl b/src/NCMSolver.jl
index a22bc69..02a67d3 100644
--- a/src/NCMSolver.jl
+++ b/src/NCMSolver.jl
@@ -14,6 +14,7 @@ Common interface for solving NCM problems. Algorithm-specific cache is stored in
 - `maxiters`: The number of iterations allowed. Defaults to `size(A,1)`
 - `ensure_pd`: Checks (and corrects) that the resulting matrix is positive definite.
   Defaults to `false`.
+- `min_eigenvalue`: The minimum eigenvalue to enforce when `ensure_pd` == true.
 - `verbose`: Whether to print extra information. Defaults to `false`.
 - `mask`: The fixed-element mask, or `nothing` if unmasked.
 - `A_orig`: The original values of A to be used when a mask is supplied.
@@ -28,6 +29,7 @@ mutable struct NCMSolver{TA, P, Talg, Tc, Ttol, Tm}
     reltol::Ttol    # relative tolerance for convergence
     maxiters::Int   # maximum number of iterations
     ensure_pd::Bool # ensures that the resulting matrix is positive definite
+    min_eigenvalue::Union{Nothing, Real} # the minimum eigenvalue to enforce
     verbose::Bool   # whether to print extra information
     mask::Tm        # fixed-element mask, or nothing
     A_orig::TA      # a copy of A, or an alias of A if no mask is given
@@ -38,7 +40,7 @@ end
 
 Get the default algorithm type for a given input matrix.
 """
-default_algtype(prob::NCMProblem) = prob.mask === nothing ? Newton : AlternatingProjections
+default_algtype(prob::NCMProblem) = prob.mask === nothing ? Newton : AcceleratedAP
 
 """
     init(prob, alg, args...; kwargs...)
@@ -84,6 +86,7 @@ function CommonSolve.init(
         convert_f16::Bool = false,
         force_f16::Bool = false,
         ensure_pd::Bool = false,
+        min_eigenvalue = nothing,
         verbose::Bool = false,
         kwargs...
     )
@@ -113,6 +116,8 @@ function CommonSolve.init(
     A = prob.A
     p = prob.p
 
+    T = eltype(A)
+
     A = if alias_A
         verbose && println("Aliasing A")
         A
@@ -159,7 +164,7 @@ function CommonSolve.init(
         end
     end
 
-    if eltype(A) === Float16 && !supports_float16(alg)
+    if T === Float16 && !supports_float16(alg)
         if convert_f16
             verbose &&
                 println(
@@ -188,15 +193,35 @@ function CommonSolve.init(
     A_orig = mask === nothing ? A : copy(A)
 
     # Guard against type mismatch for user-specified reltol/abstol
-    reltol = real(eltype(A))(reltol)
-    abstol = real(eltype(A))(abstol)
+    reltol = real(T)(reltol)
+    reltol = max(reltol, sqrt(eps(T)))
+    abstol = real(T)(abstol)
+    abstol = max(abstol, eps(T))
+
+    min_eigenvalue = if min_eigenvalue === nothing
+        if ensure_pd
+            if mask === nothing
+                # no mask, can default to sqrt(eps(T))
+                sqrt(eps(T))
+            else
+                # be more conservative about the min eigenvalue when there is a mask
+                sqrt(sqrt(eps(T)))
+            end
+        else
+            # no checks for PD -> min_eigenvalue is not used
+            nothing
+        end
+    else
+        # user explicitly set min_eigenvalue. Just ensure that it is Real
+        real(T)(min_eigenvalue)
+    end
 
     cacheval = init_cacheval(alg, A; maxiters = maxiters, abstol = abstol, reltol = reltol, verbose = verbose)
     isfresh = true
     Tc = typeof(cacheval)
 
     solver = NCMSolver{typeof(A), typeof(p), typeof(alg), Tc, typeof(reltol), Union{Nothing, typeof(mask)}}(
-        A, p, alg, cacheval, isfresh, abstol, reltol, maxiters, ensure_pd, verbose, mask, A_orig
+        A, p, alg, cacheval, isfresh, abstol, reltol, maxiters, ensure_pd, min_eigenvalue, verbose, mask, A_orig
     )
 
     return solver
diff --git a/src/NearestCorrelationMatrix.jl b/src/NearestCorrelationMatrix.jl
index 79dc638..fa23525 100644
--- a/src/NearestCorrelationMatrix.jl
+++ b/src/NearestCorrelationMatrix.jl
@@ -20,6 +20,7 @@ include("simple_interface.jl")
 include("algorithms/Newton.jl")
 include("algorithms/DirectProjection.jl")
 include("algorithms/AlternatingProjections.jl")
+include("algorithms/AcceleratedAP.jl")
 include("algorithms/JuMPAlgorithm.jl")
 
 export
@@ -39,9 +40,10 @@ export
     nearest_cor,
     nearest_cor!,
     # algorithms
-    Newton,
-    DirectProjection,
+    AcceleratedAP,
     AlternatingProjections,
-    JuMPAlgorithm
+    DirectProjection,
+    JuMPAlgorithm,
+    Newton
 
 end
diff --git a/src/algorithms/AlternatingProjectionsAA.jl b/src/algorithms/AcceleratedAP.jl
similarity index 56%
rename from src/algorithms/AlternatingProjectionsAA.jl
rename to src/algorithms/AcceleratedAP.jl
index 66eeed7..5523dde 100644
--- a/src/algorithms/AlternatingProjectionsAA.jl
+++ b/src/algorithms/AcceleratedAP.jl
@@ -1,33 +1,40 @@
-struct AlternatingProjectionsAA{A, K} <: NCMAlgorithm
+"""
+    AcceleratedAP(; tau=0, m=2)
+
+The alternating projections algorithm with Anderson acceleration applied. Should converge
+in roughly half the number of steps of the standard alternating projections algorithm.
+"""
+struct AcceleratedAP{A, K} <: NCMAlgorithm
     tau::Real
     m::Int
     args::A
     kwargs::K
 end
 
-function AlternatingProjectionsAA(args...; tau::Real = 0, m::Int = 2, kwargs...)
-    return AlternatingProjectionsAA(tau, m, args, kwargs)
+function AcceleratedAP(args...; tau::Real = 0, m::Int = 2, kwargs...)
+    return AcceleratedAP(tau, m, args, kwargs)
 end
 
-default_iters(::AlternatingProjectionsAA, A) = clamp(size(A, 1), 20, 200)
-modifies_in_place(::AlternatingProjectionsAA) = true
-supports_float16(::AlternatingProjectionsAA) = true
-supports_symmetric(::AlternatingProjectionsAA) = false
-supports_parameterless_construction(::Type{AlternatingProjectionsAA}) = true
+default_iters(::AcceleratedAP, A) = clamp(size(A, 1), 20, 200)
+modifies_in_place(::AcceleratedAP) = true
+supports_float16(::AcceleratedAP) = true
+supports_symmetric(::AcceleratedAP) = false
+supports_parameterless_construction(::Type{<:AcceleratedAP}) = true
+supports_mask(::Type{<:AcceleratedAP}) = true
 
-function autotune(::Type{AlternatingProjectionsAA}, prob::NCMProblem)
-    return AlternatingProjectionsAA(; tau = eps(eltype(prob.A)), m = 2)
+function autotune(::Type{<:AcceleratedAP}, prob::NCMProblem)
+    return AcceleratedAP(; tau = sqrt(eps(eltype(prob.A))), m = 2)
 end
 
-function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kwargs...)
+function CommonSolve.solve!(solver::NCMSolver, alg::AcceleratedAP; kwargs...)
     A = solver.A
     n = size(A, 1)
-    size(A, 2) == n || throw(DimensionMismatch("Input matrix A must be square."))
 
     T = eltype(A)
     m = alg.m
-    tol = solver.reltol
-    maxiter = solver.maxiters
+    tau = convert(T, alg.tau)
+    mask = solver.mask
+    A_orig = solver.A_orig
 
     # Initialize working matrices
     X = copy(A)
@@ -35,6 +42,7 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kw
     S = zeros(T, n, n)
     R = similar(A)
     G = similar(A)
+    scratch = similar(A)
 
     # Pre-allocate memory for Anderson Acceleration history
     vec_dim = n * n
@@ -47,36 +55,26 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kw
     m_eff = 0
 
     iter = 0
-    converged = false
-    rel_err = 0.0
+    resid = Inf
 
-    while iter < maxiter
+    while iter < solver.maxiters
         iter += 1
 
-        # R = Y - S
-        @. R = Y - S
-
-        # X = P_S(R) : Project onto Positive Semidefinite Cone S+
-        X .= project_s(R)
-
-        # S = X - R : Update Dykstra correction
-        @. S = X - R
-
-        # G = P_U(X) : Project onto Unit Diagonal U
+        R .= Y .- S
+        project_psd!(X, R, tau, scratch)
+        S .= X .- R
         copyto!(G, X)
-        for i in 1:n
-            G[i, i] = one(T)
+
+        if mask !== nothing
+            project_fixed!(G, A_orig, mask)
         end
 
-        # Relative residual error check
-        rel_err = norm(X .- G, 2) / max(one(T), norm(X, 2))
+        # project unit after projecting fixed to ensure that the unit diagonal is preserved
+        project_unit!(G)
 
-        if solver.verbose
-            println("Iter $iter: rel_err = $rel_err")
-        end
+        resid = norm(X .- G) / norm(X)
 
-        if rel_err <= tol
-            converged = true
+        if resid <= solver.reltol
             break
         end
 
@@ -131,11 +129,10 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kw
         end
     end
 
-    # Ensure output X strictly satisfies unit diagonal and exact symmetry
-    for i in 1:n
-        X[i, i] = one(T)
+    if mask !== nothing
+        project_fixed!(X, A_orig, mask)
     end
-    X .= (X .+ X') ./ 2
+    project_unit!(X)
 
-    return build_ncm_solution(alg, X, rel_err, solver; iters = iter)
+    return build_ncm_solution(alg, X, resid, solver; iters = iter)
 end
diff --git a/src/algorithms/AlternatingProjections.jl b/src/algorithms/AlternatingProjections.jl
index 64a3d46..e7e9702 100644
--- a/src/algorithms/AlternatingProjections.jl
+++ b/src/algorithms/AlternatingProjections.jl
@@ -32,11 +32,11 @@ default_iters(::AlternatingProjections, A) = clamp(size(A, 1), 20, 200)
 modifies_in_place(::AlternatingProjections) = true
 supports_float16(::AlternatingProjections) = true
 supports_symmetric(::AlternatingProjections) = false
-supports_parameterless_construction(::Type{AlternatingProjections}) = true
-supports_mask(::AlternatingProjections) = true
+supports_parameterless_construction(::Type{<:AlternatingProjections}) = true
+supports_mask(::Type{<:AlternatingProjections}) = true
 
-function autotune(::Type{AlternatingProjections}, prob::NCMProblem)
-    return AlternatingProjections(; tau = eps(eltype(prob.A)))
+function autotune(::Type{<:AlternatingProjections}, prob::NCMProblem)
+    return AlternatingProjections(; tau = sqrt(eps(eltype(prob.A))))
 end
 
 function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwargs...)
@@ -54,7 +54,9 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwar
     iter = 0
     resid = Inf
 
-    while iter < solver.maxiters && resid ≥ solver.reltol
+    while iter < solver.maxiters
+        iter += 1
+
         R .= Y .- ΔS
         project_psd!(X, R, tau, scratch)
         ΔS .= X .- R
@@ -68,7 +70,10 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwar
         project_unit!(Y)
 
         resid = norm(Y .- X) / norm(Y)
-        iter += 1
+
+        if resid ≤ solver.reltol
+            break
+        end
     end
 
     return build_ncm_solution(alg, Y, resid, solver; iters = iter)
diff --git a/src/algorithms/DirectProjection.jl b/src/algorithms/DirectProjection.jl
index e8a96bf..7a2b383 100644
--- a/src/algorithms/DirectProjection.jl
+++ b/src/algorithms/DirectProjection.jl
@@ -1,5 +1,5 @@
 """
-    DirectProjection(; tau=eps())
+    DirectProjection(args...; tau=0, kwargs...)
 
 Single step projection of the input matrix into the set of correlation matrices. Useful when
 a "close" correlation matrix is needed without concern for it being the most optimal.
@@ -11,6 +11,7 @@ struct DirectProjection{A, K} <: NCMAlgorithm
     tau::Real
     args::A
     kwargs::K
+
 end
 
 function DirectProjection(args...; tau::Real = 0, kwargs...)
@@ -24,7 +25,7 @@ supports_parameterless_construction(::Type{DirectProjection}) = true
 
 autotune(::Type{DirectProjection}, prob::NCMProblem) = _autotune(DirectProjection, prob.A)
 
-function _autotune(::Type{DirectProjection}, A::AbstractMatrix{Float64})
+function _autotune(::Type{DirectProjection}, ::AbstractMatrix{Float64})
     return DirectProjection(; tau = 1.0e-12)
 end
 
@@ -61,12 +62,14 @@ function _autotune(::Type{DirectProjection}, A::AbstractMatrix{Float16})
 end
 
 function CommonSolve.solve!(solver::NCMSolver, alg::DirectProjection; kwargs...)
-    X = solver.A
-    T = eltype(X)
-    tau = max(T(alg.tau), zero(T))
+    A = solver.A
+    X = copy(A)
+    tau = convert(eltype(X), alg.tau)
 
     project_psd!(X, tau)
     cov2cor!(X)
 
-    return build_ncm_solution(alg, X, nothing, solver; iters = 1)
+    resid = norm(X .- A) / norm(X)
+
+    return build_ncm_solution(alg, X, resid, solver; iters = 1)
 end
diff --git a/src/algorithms/Newton.jl b/src/algorithms/Newton.jl
index 9b19994..cc7e31f 100644
--- a/src/algorithms/Newton.jl
+++ b/src/algorithms/Newton.jl
@@ -33,7 +33,23 @@ function Newton(
 end
 
 autotune(::Type{Newton}, prob::NCMProblem) = _autotune(Newton, prob.A)
-_autotune(::Type{Newton}, A::AbstractMatrix{Float64}) = Newton(; tau = 1.0e-12)
+function _autotune(::Type{Newton}, A::AbstractMatrix{Float64})
+    n = size(A, 1)
+
+    tau = if n ≤ 50
+        1.0e-12
+    elseif n ≤ 100
+        1.0e-11
+    elseif n ≤ 500
+        1.0e-10
+    elseif n ≤ 1000
+        1.0e-8
+    else
+        1.0e-6
+    end
+
+    return Newton(; tau = tau)
+end
 
 function _autotune(::Type{Newton}, A::AbstractMatrix{Float32})
     n = size(A, 1)
diff --git a/test/runtests.jl b/test/runtests.jl
index 3c9e63e..b38aa9b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,7 @@ include("datadeps_registration.jl")
 # Algorithm Robustnes
 @safetestset "Constructors" include("test_constructors.jl")
 @safetestset "Convergence" include("test_convergence.jl")
+@safetestset "Fixed Element Masking" include("test_masking.jl")
 @safetestset "Real World Data" include("test_real_world.jl")
 
 # Extension Packages
diff --git a/test/test_masking.jl b/test/test_masking.jl
new file mode 100644
index 0000000..845f1cd
--- /dev/null
+++ b/test/test_masking.jl
@@ -0,0 +1,34 @@
+using Test
+using LinearAlgebra
+using InteractiveUtils
+using NearestCorrelationMatrix
+import NearestCorrelationMatrix as NCM
+
+include("Datasets.jl")
+using .Datasets
+
+include("CustomTestMacros.jl")
+using .CustomTestMacros
+
+masking_algs = filter(alg -> NCM.supports_mask(alg), subtypes(NCMAlgorithm))
+supported_types = (Float64, Float32, Float16)
+
+for algtype in masking_algs, T in supported_types
+    @testset "$(NCM.alg_name(algtype)) - $T" begin
+        A, m = usgs13()
+        A = convert(Matrix{T}, A)
+        X = copy(A)
+        prob = NCMProblem(X; mask = m)
+        alg = autotune(algtype, prob)
+        solver = init(prob, alg)
+        @test solver.mask !== nothing
+        @test Base.mightalias(solver.A_orig, A) == false
+        sol = solve!(solver)
+        @test isapprox(sol.X, sol.X')
+        @test isapprox(diag(sol.X), ones(T, size(sol.X, 1)))
+        @test all(x -> prevfloat(-one(T)) <= x <= nextfloat(one(T)), sol.X)
+        evals = eigvals(sol.X)
+        @test all(>=(-sqrt(sqrt(eps(T)))), evals)
+        @test sol.X[m] == A[m]
+    end
+end
diff --git a/test/test_real_world.jl b/test/test_real_world.jl
index d9339af..d07bd4c 100644
--- a/test/test_real_world.jl
+++ b/test/test_real_world.jl
@@ -1,5 +1,7 @@
 using Test
+using InteractiveUtils
 using NearestCorrelationMatrix
+import NearestCorrelationMatrix as NCM
 using NearestCorrelationMatrix.Internals: iscorrelation
 
 include("Datasets.jl")
@@ -8,106 +10,119 @@ using .Datasets
 include("CustomTestMacros.jl")
 using .CustomTestMacros
 
-@testset "BCCD16" begin
-    A = bccd16()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "BEYU11" begin
-    A = beyu11()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "BHWI01" begin
-    A = bhwi01()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "COR1399" begin
-    A = cor1399()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "COR3120" begin
-    A = cor3120()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "FING97" begin
-    # unmasked
-    A, _ = fing97()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-
-    A, m = fing97()
-    @test !iscorrelation(A)
-    B = nearest_cor(A, AlternatingProjections(); mask = m)
-    @test_iscorrelation B
-    @test B[m] == A[m]
-end
-
-@testset "HIGH02" begin
-    A = high02()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "MMB13" begin
-    A = mmb13()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "TEC03" begin
-    A = tec03()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "TYDA99R1" begin
-    A = tyda99r1()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "TYDA99R2" begin
-    A = tyda99r2()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "TYDA99R3" begin
-    A = tyda99r3()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-end
-
-@testset "USGS13" begin
-    A, _ = usgs13()
-    @test !iscorrelation(A)
-    nearest_cor!(A)
-    @test_iscorrelation A
-
-    A, m = usgs13()
-    @test !iscorrelation(A)
-    B = nearest_cor(A, AlternatingProjections(tau = 1.0e-6); mask = m)
-    @test_iscorrelation B
-    @test B[m] == A[m]
+internal_algtypes = setdiff(subtypes(NCMAlgorithm), (JuMPAlgorithm,))
+fast_algtypes = (Newton, DirectProjection)
+
+for alg in internal_algtypes
+    @testset "$(NCM.alg_name(alg))" begin
+        @testset "BEYU11" begin
+            A = beyu11()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "BHWI01" begin
+            A = bhwi01()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "COR1399" begin
+            A = cor1399()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        if alg ∈ fast_algtypes
+            @testset "BCCD16" begin
+                A = bccd16()
+                @test !iscorrelation(A)
+                nearest_cor!(A, alg)
+                @test_iscorrelation A
+            end
+
+            @testset "COR3120" begin
+                A = cor3120()
+                @test !iscorrelation(A)
+                nearest_cor!(A, alg)
+                @test_iscorrelation A
+            end
+        end
+
+        @testset "FING97" begin
+            # unmasked
+            A, _ = fing97()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+
+            if NCM.supports_mask(alg)
+                A, m = fing97()
+                @test !iscorrelation(A)
+                B = nearest_cor(A, alg; mask = m)
+                @test_iscorrelation B
+                @test B[m] == A[m]
+            end
+        end
+
+        @testset "HIGH02" begin
+            A = high02()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "MMB13" begin
+            A = mmb13()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "TEC03" begin
+            A = tec03()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "TYDA99R1" begin
+            A = tyda99r1()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "TYDA99R2" begin
+            A = tyda99r2()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "TYDA99R3" begin
+            A = tyda99r3()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+        end
+
+        @testset "USGS13" begin
+            A, _ = usgs13()
+            @test !iscorrelation(A)
+            nearest_cor!(A, alg)
+            @test_iscorrelation A
+
+            if NCM.supports_mask(alg)
+                A, m = usgs13()
+                @test !iscorrelation(A)
+                B = nearest_cor(A, alg; mask = m)
+                @test_iscorrelation B
+                @test B[m] == A[m]
+            end
+        end
+    end
 end