From 569d9d9c405a99e74f6f5e92a8ae57c38f20a5e1 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Mon, 29 Aug 2022 09:17:23 -0400
Subject: [PATCH] Improve rNorm estimate in CR

---
 src/cr.jl       | 10 ++++++++--
 test/test_cr.jl |  1 -
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/cr.jl b/src/cr.jl
index 003b5a595..c678c7d29 100644
--- 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
index 358ea3ce4..ffe154b7c 100644
--- 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)