diff --git a/src/cr.jl b/src/cr.jl
index 65e43c63e..fe51955b7 100644
--- a/src/cr.jl
+++ b/src/cr.jl
@@ -335,7 +335,7 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema
 
       kaxpy!(n, α, p, x)
       xNorm = knorm(n, x)
-      xNorm ≈ radius && (on_boundary = true)
+      xNorm ≈ radius > 0 && (on_boundary = true)
       kaxpy!(n, -α, Mq, r)  # residual
       if MisI
         rNorm² = kdotr(n, r, r)
@@ -401,11 +401,11 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema
 
     # Termination status
     tired               && (status = "maximum number of iterations exceeded")
-    on_boundary         && (status = "on trust-region boundary")
-    npcurv              && (status = "nonpositive curvature")
     solved              && (status = "solution good enough given atol and rtol")
     user_requested_exit && (status = "user-requested exit")
     overtimed           && (status = "time limit exceeded")
+    npcurv              && (status = "nonpositive curvature")
+    on_boundary         && (status = "on trust-region boundary")
 
     # Update x
     warm_start && kaxpy!(n, one(FC), Δx, x)
diff --git a/test/test_cr.jl b/test/test_cr.jl
index 52da6c013..6ea9de7d2 100644
--- a/test/test_cr.jl
+++ b/test/test_cr.jl
@@ -66,6 +66,14 @@
       A, b = symmetric_indefinite(FC=FC)
       x, stats = cr(A, b, linesearch=true)
       @test stats.status == "nonpositive curvature"
+      # @test real(dot(x, A * x)) ≤ 0
+
+      # Test on trust-region boundary when radius > 0
+      A, b = symmetric_indefinite(FC=FC, shift = 5)
+      x, stats = cr(A, b, radius = one(Float64))
+      @test stats.status == "on trust-region boundary"
+      @test real(dot(x, A * x)) ≤ 0
+      @test norm(x) ≈ 1.0
 
       # Test Linesearch which would stop on the first call since A is negative definite
       A, b = symmetric_indefinite(FC=FC; shift = 5)
@@ -90,7 +98,7 @@
       @test stats.solved == true
 
  
-      # test callback function
+      # Test callback function
       A, b = symmetric_definite(FC=FC)
       solver = CrSolver(A, b)
       tol = 1.0e-1