Skip to content

Commit 52c1619

Browse files
farhadrclassdpo
authored andcommitted
fixed the issue and moved rNorm
1 parent a6cb197 commit 52c1619

File tree

3 files changed

+12
-27
lines changed

3 files changed

+12
-27
lines changed

src/cr.jl

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -150,27 +150,28 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema
150150
end
151151
MisI && kcopy!(n, r, p) # r ← p
152152
MisI || mulorldiv!(r, M, p, ldiv)
153-
mul!(Ar, A, r)
154-
ρ = kdotr(n, r, Ar)
155-
153+
156154
rNorm = knorm_elliptic(n, r, p) # ‖r‖
157155
history && push!(rNorms, rNorm) # Values of ‖r‖
158156

159157
if rNorm == 0
160158
stats.niter = 0
161159
stats.solved, stats.inconsistent = true, false
162160
stats.timer = start_time |> ktimer
163-
stats.status = "x = 0 is a zero-residual solution"
161+
stats.status = "x is a zero-residual solution"
164162
history && push!(ArNorms, zero(T))
165163
solver.warm_start = false
166164
return solver
167165
end
166+
167+
mul!(Ar, A, r)
168+
ρ = kdotr(n, r, Ar)
168169

169170
if ρ == 0
170171
stats.niter = 0
171172
stats.solved, stats.inconsistent = true, false
172173
stats.timer = start_time |> ktimer
173-
stats.status = "0 is a zero-curvature direction"
174+
stats.status = "b is a zero-curvature direction"
174175
history && push!(ArNorms, zero(T))
175176
solver.warm_start = false
176177
linesearch && kcopy!(n, x, b) # x ← b

test/test_cr.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@
5252
A, b = zero_rhs(FC=FC)
5353
(x, stats) = cr(A, b)
5454
@test norm(x) == 0
55-
@test stats.status == "x = 0 is a zero-residual solution"
55+
@test stats.status == "x is a zero-residual solution"
5656

5757
# Test with Jacobi (or diagonal) preconditioner
5858
A, b, M = square_preconditioned(FC=FC)
@@ -68,7 +68,7 @@
6868
@test stats.status == "nonpositive curvature"
6969

7070
# Test Linesearch which would stop on the first call since A is negative definite
71-
A, b = indefinite_system(FC=FC)
71+
A, b = symmetric_indefinite(FC=FC; shift = 5)
7272
x, stats = cr(A, b, linesearch=true)
7373
@test stats.status == "nonpositive curvature"
7474
@test stats.niter == 0
@@ -78,14 +78,14 @@
7878
# Test when b^TAb=0 and linesearch is true
7979
A, b = system_zero_quad(FC=FC)
8080
x, stats = cr(A, b, linesearch=true)
81-
@test stats.status == "0 is a zero-curvature direction"
81+
@test stats.status == "b is a zero-curvature direction"
8282
@test all(x .== b)
8383
@test stats.solved == true
8484

8585
# Test when b^TAb=0 and linesearch is false
8686
A, b = system_zero_quad(FC=FC)
8787
x, stats = cr(A,b, linesearch=false)
88-
@test stats.status == "0 is a zero-curvature direction"
88+
@test stats.status == "b is a zero-curvature direction"
8989
@test norm(x) == zero(FC)
9090
@test stats.solved == true
9191

test/test_utils.jl

Lines changed: 2 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -14,29 +14,13 @@ function symmetric_definite(n :: Int=10; FC=Float64)
1414
end
1515

1616
# Symmetric and indefinite systems.
17-
function symmetric_indefinite(n :: Int=10; FC=Float64)
17+
function symmetric_indefinite(n :: Int=10; FC=Float64, shift=0)
1818
α = FC <: Complex ? FC(im) : one(FC)
19-
A = spdiagm(-1 => α * ones(FC, n-1), 0 => ones(FC, n), 1 => conj(α) * ones(FC, n-1))
19+
A = spdiagm(-1 => α * ones(FC, n-1), 0 => ones(FC, n), 1 => conj(α) * ones(FC, n-1))- shift * eye(n)
2020
b = A * FC[1:n;]
2121
return A, b
2222
end
2323

24-
25-
# Symmetric and indefinite system with negative curvature.
26-
# need this to test CR with Linesearch and nonpositive curvature at the first iteration
27-
function indefinite_system(n::Int=10; FC=Float64, shift=5)
28-
α = FC <: Complex ? FC(im) : one(FC)
29-
# Build the tridiagonal matrix with 1's on the diagonal and α on the off-diagonals, and force negative curvature.
30-
A = spdiagm(
31-
-1 => α * ones(FC, n-1),
32-
0 => ones(FC, n),
33-
1 => conj(α) * ones(FC, n-1)
34-
) - shift * eye(n)
35-
36-
b = A * FC[1:n;]
37-
return A, b
38-
end
39-
4024
"""
4125
system_zero_quad(n::Int=2; FC=Float64)
4226

0 commit comments

Comments
 (0)