Improve rNorm estimate in CR

JuliaSmoothOptimizers · Aug 29, 2022 · 569d9d9 · 569d9d9
1 parent a08e934
commit 569d9d9
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 3 deletions.
diff --git a/src/cr.jl b/src/cr.jl
@@ -281,8 +281,14 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
     xNorm = @knrm2(n, x)
     xNorm ≈ radius && (on_boundary = true)
     @kaxpy!(n, -α, Mq, r) # residual
-    rNorm² = abs(rNorm² - α * ρ)
-    rNorm = sqrt(rNorm²)
+    if MisI
+      rNorm² = @kdotr(n, r, r)
+      rNorm = sqrt(rNorm²)
+    else
+      ω = sqrt(α) * sqrt(ρ)
+      rNorm = sqrt(abs(rNorm + ω)) * sqrt(abs(rNorm - ω))
+      rNorm² = rNorm * rNorm  # rNorm² = rNorm² - α * ρ
+    end
     history && push!(rNorms, rNorm)
     mul!(Ar, A, r)
     ArNorm = @knrm2(n, Ar)

diff --git a/test/test_cr.jl b/test/test_cr.jl
@@ -63,7 +63,6 @@
       @test(resid ≤ cr_tol)
       @test(stats.solved)
 
-
       # test callback function
       A, b = symmetric_definite(FC=FC)
       solver = CrSolver(A, b)