From c5b440ee678f37a06235aa54db91cd7fa550a098 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 25 Jan 2021 00:05:06 -0500 Subject: [PATCH 01/15] Add conjugate orthogonal conjugate gradient (COCG) method for complex symmetric linear systems --- src/IterativeSolvers.jl | 1 + src/cocg.jl | 241 ++++++++++++++++++++++++++++++++++++++++ test/cocg.jl | 71 ++++++++++++ test/runtests.jl | 3 + 4 files changed, 316 insertions(+) create mode 100644 src/cocg.jl create mode 100644 test/cocg.jl diff --git a/src/IterativeSolvers.jl b/src/IterativeSolvers.jl index c3f11860..e479c72a 100644 --- a/src/IterativeSolvers.jl +++ b/src/IterativeSolvers.jl @@ -17,6 +17,7 @@ include("hessenberg.jl") include("stationary.jl") include("stationary_sparse.jl") include("cg.jl") +include("cocg.jl") include("minres.jl") include("bicgstabl.jl") include("gmres.jl") diff --git a/src/cocg.jl b/src/cocg.jl new file mode 100644 index 00000000..9cd9bd76 --- /dev/null +++ b/src/cocg.jl @@ -0,0 +1,241 @@ +import Base: iterate +using Printf +export cocg, cocg!, COCGIterable, PCOCGIterable, cocg_iterator!, COCGStateVariables + +mutable struct COCGIterable{matT, solT, vecT, numT <: Real, paramT <: Number} + A::matT + x::solT + r::vecT + c::vecT + u::vecT + tol::numT + residual::numT + rr_prev::paramT + maxiter::Int + mv_products::Int +end + +mutable struct PCOCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number} + Pl::precT + A::matT + x::solT + r::vecT + c::vecT + u::vecT + tol::numT + residual::numT + rc_prev::paramT + maxiter::Int + mv_products::Int +end + +@inline converged(it::Union{COCGIterable, PCOCGIterable}) = it.residual ≤ it.tol + +@inline start(it::Union{COCGIterable, PCOCGIterable}) = 0 + +@inline done(it::Union{COCGIterable, PCOCGIterable}, iteration::Int) = iteration ≥ it.maxiter || converged(it) + + +################# +# Ordinary COCG # +################# + +function iterate(it::COCGIterable, iteration::Int=start(it)) + # Check for termination first + if done(it, iteration) + return nothing + end + + # u := r + βu (almost an axpy) + rr = sum(rₖ^2 for rₖ in it.r) # rᵀ * r + β = rr / it.rr_prev + it.u .= it.r .+ β .* it.u + + # c = A * u + mul!(it.c, it.A, it.u) + uc = sum(uₖ*cₖ for (uₖ,cₖ) in zip(it.u,it.c)) # uᵀ * c + α = rr / uc + + # Improve solution and residual + it.rr_prev = rr + it.x .+= α .* it.u + it.r .-= α .* it.c + + it.residual = norm(it.r) + + # Return the residual at item and iteration number as state + it.residual, iteration + 1 +end + +####################### +# Preconditioned COCG # +####################### + +function iterate(it::PCOCGIterable, iteration::Int=start(it)) + # Check for termination first + if done(it, iteration) + return nothing + end + + # Apply left preconditioner: c = Pl⁻¹ r + ldiv!(it.c, it.Pl, it.r) + + # u := c + βu (almost an axpy) + rc = sum(rₖ*cₖ for (rₖ,cₖ) in zip(it.r,it.c)) # rᵀ * c + β = rc / it.rc_prev + it.u .= it.c .+ β .* it.u + + # c = A * u + mul!(it.c, it.A, it.u) + uc = sum(uₖ*cₖ for (uₖ,cₖ) in zip(it.u,it.c)) # uᵀ * c + α = rc / uc + + # Improve solution and residual + it.rc_prev = rc + it.x .+= α .* it.u + it.r .-= α .* it.c + + it.residual = norm(it.r) + + # Return the residual at item and iteration number as state + it.residual, iteration + 1 +end + +# Utility functions + +""" +Intermediate COCG state variables to be used inside cocg and cocg!. `u`, `r` and `c` should be of the same type as the solution of `cocg` or `cocg!`. +``` +struct COCGStateVariables{T,Tx<:AbstractArray{T}} + u::Tx + r::Tx + c::Tx +end +``` +""" +struct COCGStateVariables{T,Tx<:AbstractArray{T}} + u::Tx + r::Tx + c::Tx +end + +function cocg_iterator!(x, A, b, Pl = Identity(); + abstol::Real = zero(real(eltype(b))), + reltol::Real = sqrt(eps(real(eltype(b)))), + maxiter::Int = size(A, 2), + statevars::COCGStateVariables = COCGStateVariables(zero(x), similar(x), similar(x)), + initially_zero::Bool = false) + u = statevars.u + r = statevars.r + c = statevars.c + u .= zero(eltype(x)) + copyto!(r, b) + + # Compute r with an MV-product or not. + if initially_zero + mv_products = 0 + else + mv_products = 1 + mul!(c, A, x) + r .-= c + end + residual = norm(r) + tolerance = max(reltol * residual, abstol) + + # Return the iterable + if isa(Pl, Identity) + return COCGIterable(A, x, r, c, u, + tolerance, residual, one(eltype(r)), + maxiter, mv_products + ) + else + return PCOCGIterable(Pl, A, x, r, c, u, + tolerance, residual, one(eltype(r)), + maxiter, mv_products + ) + end +end + +""" + cocg(A, b; kwargs...) -> x, [history] + +Same as [`cocg!`](@ref), but allocates a solution vector `x` initialized with zeros. +""" +cocg(A, b; kwargs...) = cocg!(zerox(A, b), A, b; initially_zero = true, kwargs...) + +""" + cocg!(x, A, b; kwargs...) -> x, [history] + +# Arguments + +- `x`: Initial guess, will be updated in-place; +- `A`: linear operator; +- `b`: right-hand side. + +## Keywords + +- `statevars::COCGStateVariables`: Has 3 arrays similar to `x` to hold intermediate results; +- `initially_zero::Bool`: If `true` assumes that `iszero(x)` so that one + matrix-vector product can be saved when computing the initial + residual vector; +- `Pl = Identity()`: left preconditioner of the method. Should be symmetric, + positive-definite like `A`; +- `abstol::Real = zero(real(eltype(b)))`, + `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative + tolerance for the stopping condition + `|r_k| / |r_0| ≤ max(reltol * resnorm, abstol)`, where `r_k = A * x_k - b` + is the residual in the `k`th iteration; +- `maxiter::Int = size(A,2)`: maximum number of iterations; +- `verbose::Bool = false`: print method information; +- `log::Bool = false`: keep track of the residual norm in each iteration. + +# Output + +**if `log` is `false`** + +- `x`: approximated solution. + +**if `log` is `true`** + +- `x`: approximated solution. +- `ch`: convergence history. + +**ConvergenceHistory keys** + +- `:tol` => `::Real`: stopping tolerance. +- `:resnom` => `::Vector`: residual norm at each iteration. +""" +function cocg!(x, A, b; + abstol::Real = zero(real(eltype(b))), + reltol::Real = sqrt(eps(real(eltype(b)))), + maxiter::Int = size(A, 2), + log::Bool = false, + statevars::COCGStateVariables = COCGStateVariables(zero(x), similar(x), similar(x)), + verbose::Bool = false, + Pl = Identity(), + kwargs...) + history = ConvergenceHistory(partial = !log) + history[:abstol] = abstol + history[:reltol] = reltol + log && reserve!(history, :resnorm, maxiter + 1) + + # Actually perform COCG + iterable = cocg_iterator!(x, A, b, Pl; abstol = abstol, reltol = reltol, maxiter = maxiter, + statevars = statevars, kwargs...) + if log + history.mvps = iterable.mv_products + end + for (iteration, item) = enumerate(iterable) + if log + nextiter!(history, mvps = 1) + push!(history, :resnorm, iterable.residual) + end + verbose && @printf("%3d\t%1.2e\n", iteration, iterable.residual) + end + + verbose && println() + log && setconv(history, converged(iterable)) + log && shrink!(history) + + log ? (iterable.x, history) : iterable.x +end diff --git a/test/cocg.jl b/test/cocg.jl new file mode 100644 index 00000000..19b32f6e --- /dev/null +++ b/test/cocg.jl @@ -0,0 +1,71 @@ +module TestCOCG + +using IterativeSolvers +using Test +using Random +using LinearAlgebra + +@testset ("Conjugate Orthogonal Conjugate Gradient") begin + +Random.seed!(1234321) +n = 20 + +@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64) + A = rand(T, n, n) + A = A + transpose(A) + 15I + x = ones(T, n) + b = A * x + + reltol = √eps(real(T)) + + # Solve without preconditioner + x1, his1 = cocg(A, b, reltol = reltol, maxiter = 100, log = true) + @test isa(his1, ConvergenceHistory) + @test norm(A * x1 - b) / norm(b) ≤ reltol + + # With an initial guess + x_guess = rand(T, n) + x2, his2 = cocg!(x_guess, A, b, reltol = reltol, maxiter = 100, log = true) + @test isa(his2, ConvergenceHistory) + @test x2 == x_guess + @test norm(A * x2 - b) / norm(b) ≤ reltol + + # The following tests fails CI on Windows and Ubuntu due to a + # `SingularException(4)` + if T == Float32 && (Sys.iswindows() || Sys.islinux()) + continue + end + # Do an exact LU decomp of a nearby matrix + F = lu(A + rand(T, n, n)) + x3, his3 = cocg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true) + @test norm(A * x3 - b) / norm(b) ≤ reltol +end + +@testset "Termination criterion" begin + for T in (Float32, Float64, ComplexF32, ComplexF64) + A = T[ 2 -1 0 + -1 2 -1 + 0 -1 2] + n = size(A, 2) + b = ones(T, n) + x0 = A \ b + perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n] + + # If the initial residual is small and a small relative tolerance is used, + # many iterations are necessary + x = x0 + perturbation + initial_residual = norm(A * x - b) + x, ch = cocg!(x, A, b, log=true) + @test 2 ≤ niters(ch) ≤ n + + # If the initial residual is small and a large absolute tolerance is used, + # no iterations are necessary + x = x0 + perturbation + initial_residual = norm(A * x - b) + x, ch = cocg!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true) + @test niters(ch) == 0 + end +end +end + +end # module diff --git a/test/runtests.jl b/test/runtests.jl index e612a53a..01e6e582 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,6 +12,9 @@ include("stationary.jl") #Conjugate gradients include("cg.jl") +#Conjugate Orthogonal Conjugate gradient +include("cocg.jl") + #BiCGStab(l) include("bicgstabl.jl") From 91208cd0dfdbac65bdea78799d3500ed2dc831ad Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 25 Jan 2021 16:52:04 -0500 Subject: [PATCH 02/15] Remove cocg.jl; implement COCG as modified CG using unconjugated dot product --- src/IterativeSolvers.jl | 1 - src/cg.jl | 52 +++++---- src/cocg.jl | 241 ---------------------------------------- test/cg.jl | 33 +++++- test/cocg.jl | 71 ------------ test/runtests.jl | 3 - 6 files changed, 63 insertions(+), 338 deletions(-) delete mode 100644 src/cocg.jl delete mode 100644 test/cocg.jl diff --git a/src/IterativeSolvers.jl b/src/IterativeSolvers.jl index e479c72a..c3f11860 100644 --- a/src/IterativeSolvers.jl +++ b/src/IterativeSolvers.jl @@ -17,7 +17,6 @@ include("hessenberg.jl") include("stationary.jl") include("stationary_sparse.jl") include("cg.jl") -include("cocg.jl") include("minres.jl") include("bicgstabl.jl") include("gmres.jl") diff --git a/src/cg.jl b/src/cg.jl index 764ada16..f95d087f 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -2,7 +2,15 @@ import Base: iterate using Printf export cg, cg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables -mutable struct CGIterable{matT, solT, vecT, numT <: Real} +# Conjugated dot product +_dot(x, ::Val{true}) = sum(abs2(xₖ) for xₖ in x) +_dot(x, y, ::Val{true}) = dot(x, y) + +# Unconjugated dot product +_dot(x, ::Val{false}) = sum(xₖ^2 for xₖ in x) +_dot(x, y, ::Val{false}) = sum(prod, zip(x,y)) + +mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, boolT <: Union{Val{true},Val{false}}} A::matT x::solT r::vecT @@ -10,12 +18,13 @@ mutable struct CGIterable{matT, solT, vecT, numT <: Real} u::vecT tol::numT residual::numT - prev_residual::numT + ρ_prev::paramT maxiter::Int mv_products::Int + conjugate_dot::boolT end -mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number} +mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, boolT <: Union{Val{true},Val{false}}} Pl::precT A::matT x::solT @@ -24,9 +33,10 @@ mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Numb u::vecT tol::numT residual::numT - ρ::paramT + ρ_prev::paramT maxiter::Int mv_products::Int + conjugate_dot::boolT end @inline converged(it::Union{CGIterable, PCGIterable}) = it.residual ≤ it.tol @@ -47,18 +57,19 @@ function iterate(it::CGIterable, iteration::Int=start(it)) end # u := r + βu (almost an axpy) - β = it.residual^2 / it.prev_residual^2 + ρ = _dot(it.r, it.conjugate_dot) + β = ρ / it.ρ_prev it.u .= it.r .+ β .* it.u # c = A * u mul!(it.c, it.A, it.u) - α = it.residual^2 / dot(it.u, it.c) + α = ρ / _dot(it.u, it.c, it.conjugate_dot) # Improve solution and residual + it.ρ_prev = ρ it.x .+= α .* it.u it.r .-= α .* it.c - it.prev_residual = it.residual it.residual = norm(it.r) # Return the residual at item and iteration number as state @@ -78,18 +89,17 @@ function iterate(it::PCGIterable, iteration::Int=start(it)) # Apply left preconditioner ldiv!(it.c, it.Pl, it.r) - ρ_prev = it.ρ - it.ρ = dot(it.c, it.r) - # u := c + βu (almost an axpy) - β = it.ρ / ρ_prev + ρ = _dot(it.r, it.c, it.conjugate_dot) + β = ρ / it.ρ_prev it.u .= it.c .+ β .* it.u # c = A * u mul!(it.c, it.A, it.u) - α = it.ρ / dot(it.u, it.c) + α = ρ / _dot(it.u, it.c, it.conjugate_dot) # Improve solution and residual + it.ρ_prev = ρ it.x .+= α .* it.u it.r .-= α .* it.c @@ -122,7 +132,8 @@ function cg_iterator!(x, A, b, Pl = Identity(); reltol::Real = sqrt(eps(real(eltype(b)))), maxiter::Int = size(A, 2), statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), - initially_zero::Bool = false) + initially_zero::Bool = false, + conjugate_dot::Bool = true) u = statevars.u r = statevars.r c = statevars.c @@ -142,15 +153,13 @@ function cg_iterator!(x, A, b, Pl = Identity(); # Return the iterable if isa(Pl, Identity) - return CGIterable(A, x, r, c, u, - tolerance, residual, one(residual), - maxiter, mv_products - ) + return CGIterable(A, x, r, c, u, tolerance, residual, + conjugate_dot ? one(real(eltype(r))) : one(eltype(r)), # for conjugated dot, ρ_prev remains real + maxiter, mv_products, Val(conjugate_dot)) else return PCGIterable(Pl, A, x, r, c, u, - tolerance, residual, one(eltype(x)), - maxiter, mv_products - ) + tolerance, residual, one(eltype(r)), + maxiter, mv_products, Val(conjugate_dot)) end end @@ -211,6 +220,7 @@ function cg!(x, A, b; statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), verbose::Bool = false, Pl = Identity(), + conjugate_dot::Bool = true, kwargs...) history = ConvergenceHistory(partial = !log) history[:abstol] = abstol @@ -219,7 +229,7 @@ function cg!(x, A, b; # Actually perform CG iterable = cg_iterator!(x, A, b, Pl; abstol = abstol, reltol = reltol, maxiter = maxiter, - statevars = statevars, kwargs...) + statevars = statevars, conjugate_dot = conjugate_dot, kwargs...) if log history.mvps = iterable.mv_products end diff --git a/src/cocg.jl b/src/cocg.jl deleted file mode 100644 index 9cd9bd76..00000000 --- a/src/cocg.jl +++ /dev/null @@ -1,241 +0,0 @@ -import Base: iterate -using Printf -export cocg, cocg!, COCGIterable, PCOCGIterable, cocg_iterator!, COCGStateVariables - -mutable struct COCGIterable{matT, solT, vecT, numT <: Real, paramT <: Number} - A::matT - x::solT - r::vecT - c::vecT - u::vecT - tol::numT - residual::numT - rr_prev::paramT - maxiter::Int - mv_products::Int -end - -mutable struct PCOCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number} - Pl::precT - A::matT - x::solT - r::vecT - c::vecT - u::vecT - tol::numT - residual::numT - rc_prev::paramT - maxiter::Int - mv_products::Int -end - -@inline converged(it::Union{COCGIterable, PCOCGIterable}) = it.residual ≤ it.tol - -@inline start(it::Union{COCGIterable, PCOCGIterable}) = 0 - -@inline done(it::Union{COCGIterable, PCOCGIterable}, iteration::Int) = iteration ≥ it.maxiter || converged(it) - - -################# -# Ordinary COCG # -################# - -function iterate(it::COCGIterable, iteration::Int=start(it)) - # Check for termination first - if done(it, iteration) - return nothing - end - - # u := r + βu (almost an axpy) - rr = sum(rₖ^2 for rₖ in it.r) # rᵀ * r - β = rr / it.rr_prev - it.u .= it.r .+ β .* it.u - - # c = A * u - mul!(it.c, it.A, it.u) - uc = sum(uₖ*cₖ for (uₖ,cₖ) in zip(it.u,it.c)) # uᵀ * c - α = rr / uc - - # Improve solution and residual - it.rr_prev = rr - it.x .+= α .* it.u - it.r .-= α .* it.c - - it.residual = norm(it.r) - - # Return the residual at item and iteration number as state - it.residual, iteration + 1 -end - -####################### -# Preconditioned COCG # -####################### - -function iterate(it::PCOCGIterable, iteration::Int=start(it)) - # Check for termination first - if done(it, iteration) - return nothing - end - - # Apply left preconditioner: c = Pl⁻¹ r - ldiv!(it.c, it.Pl, it.r) - - # u := c + βu (almost an axpy) - rc = sum(rₖ*cₖ for (rₖ,cₖ) in zip(it.r,it.c)) # rᵀ * c - β = rc / it.rc_prev - it.u .= it.c .+ β .* it.u - - # c = A * u - mul!(it.c, it.A, it.u) - uc = sum(uₖ*cₖ for (uₖ,cₖ) in zip(it.u,it.c)) # uᵀ * c - α = rc / uc - - # Improve solution and residual - it.rc_prev = rc - it.x .+= α .* it.u - it.r .-= α .* it.c - - it.residual = norm(it.r) - - # Return the residual at item and iteration number as state - it.residual, iteration + 1 -end - -# Utility functions - -""" -Intermediate COCG state variables to be used inside cocg and cocg!. `u`, `r` and `c` should be of the same type as the solution of `cocg` or `cocg!`. -``` -struct COCGStateVariables{T,Tx<:AbstractArray{T}} - u::Tx - r::Tx - c::Tx -end -``` -""" -struct COCGStateVariables{T,Tx<:AbstractArray{T}} - u::Tx - r::Tx - c::Tx -end - -function cocg_iterator!(x, A, b, Pl = Identity(); - abstol::Real = zero(real(eltype(b))), - reltol::Real = sqrt(eps(real(eltype(b)))), - maxiter::Int = size(A, 2), - statevars::COCGStateVariables = COCGStateVariables(zero(x), similar(x), similar(x)), - initially_zero::Bool = false) - u = statevars.u - r = statevars.r - c = statevars.c - u .= zero(eltype(x)) - copyto!(r, b) - - # Compute r with an MV-product or not. - if initially_zero - mv_products = 0 - else - mv_products = 1 - mul!(c, A, x) - r .-= c - end - residual = norm(r) - tolerance = max(reltol * residual, abstol) - - # Return the iterable - if isa(Pl, Identity) - return COCGIterable(A, x, r, c, u, - tolerance, residual, one(eltype(r)), - maxiter, mv_products - ) - else - return PCOCGIterable(Pl, A, x, r, c, u, - tolerance, residual, one(eltype(r)), - maxiter, mv_products - ) - end -end - -""" - cocg(A, b; kwargs...) -> x, [history] - -Same as [`cocg!`](@ref), but allocates a solution vector `x` initialized with zeros. -""" -cocg(A, b; kwargs...) = cocg!(zerox(A, b), A, b; initially_zero = true, kwargs...) - -""" - cocg!(x, A, b; kwargs...) -> x, [history] - -# Arguments - -- `x`: Initial guess, will be updated in-place; -- `A`: linear operator; -- `b`: right-hand side. - -## Keywords - -- `statevars::COCGStateVariables`: Has 3 arrays similar to `x` to hold intermediate results; -- `initially_zero::Bool`: If `true` assumes that `iszero(x)` so that one - matrix-vector product can be saved when computing the initial - residual vector; -- `Pl = Identity()`: left preconditioner of the method. Should be symmetric, - positive-definite like `A`; -- `abstol::Real = zero(real(eltype(b)))`, - `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative - tolerance for the stopping condition - `|r_k| / |r_0| ≤ max(reltol * resnorm, abstol)`, where `r_k = A * x_k - b` - is the residual in the `k`th iteration; -- `maxiter::Int = size(A,2)`: maximum number of iterations; -- `verbose::Bool = false`: print method information; -- `log::Bool = false`: keep track of the residual norm in each iteration. - -# Output - -**if `log` is `false`** - -- `x`: approximated solution. - -**if `log` is `true`** - -- `x`: approximated solution. -- `ch`: convergence history. - -**ConvergenceHistory keys** - -- `:tol` => `::Real`: stopping tolerance. -- `:resnom` => `::Vector`: residual norm at each iteration. -""" -function cocg!(x, A, b; - abstol::Real = zero(real(eltype(b))), - reltol::Real = sqrt(eps(real(eltype(b)))), - maxiter::Int = size(A, 2), - log::Bool = false, - statevars::COCGStateVariables = COCGStateVariables(zero(x), similar(x), similar(x)), - verbose::Bool = false, - Pl = Identity(), - kwargs...) - history = ConvergenceHistory(partial = !log) - history[:abstol] = abstol - history[:reltol] = reltol - log && reserve!(history, :resnorm, maxiter + 1) - - # Actually perform COCG - iterable = cocg_iterator!(x, A, b, Pl; abstol = abstol, reltol = reltol, maxiter = maxiter, - statevars = statevars, kwargs...) - if log - history.mvps = iterable.mv_products - end - for (iteration, item) = enumerate(iterable) - if log - nextiter!(history, mvps = 1) - push!(history, :resnorm, iterable.residual) - end - verbose && @printf("%3d\t%1.2e\n", iteration, iterable.residual) - end - - verbose && println() - log && setconv(history, converged(iterable)) - log && shrink!(history) - - log ? (iterable.x, history) : iterable.x -end diff --git a/test/cg.jl b/test/cg.jl index 9ed84b1b..1cbdfd5e 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -24,7 +24,7 @@ Random.seed!(1234321) @testset "Small full system" begin n = 10 - @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64) + @testset "Matrix{$T}, Conjugated Dot Product" for T in (Float32, Float64, ComplexF32, ComplexF64) A = rand(T, n, n) A = A' * A + I b = rand(T, n) @@ -50,6 +50,37 @@ Random.seed!(1234321) x0 = cg(A, zeros(T, n)) @test x0 == zeros(T, n) end + + @testset "Matrix{$T}, Unconjugated Dot Product" for T in (Float32, Float64, ComplexF32, ComplexF64) + A = rand(T, n, n) + A = A + transpose(A) + 15I + x = ones(T, n) + b = A * x + + reltol = √eps(real(T)) + + # Solve without preconditioner + x1, his1 = cg(A, b, reltol = reltol, maxiter = 100, log = true, conjugate_dot = false) + @test isa(his1, ConvergenceHistory) + @test norm(A * x1 - b) / norm(b) ≤ reltol + + # With an initial guess + x_guess = rand(T, n) + x2, his2 = cg!(x_guess, A, b, reltol = reltol, maxiter = 100, log = true, conjugate_dot = false) + @test isa(his2, ConvergenceHistory) + @test x2 == x_guess + @test norm(A * x2 - b) / norm(b) ≤ reltol + + # The following tests fails CI on Windows and Ubuntu due to a + # `SingularException(4)` + if T == Float32 && (Sys.iswindows() || Sys.islinux()) + continue + end + # Do an exact LU decomp of a nearby matrix + F = lu(A + rand(T, n, n)) + x3, his3 = cg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true, conjugate_dot = false) + @test norm(A * x3 - b) / norm(b) ≤ reltol + end end @testset "Sparse Laplacian" begin diff --git a/test/cocg.jl b/test/cocg.jl deleted file mode 100644 index 19b32f6e..00000000 --- a/test/cocg.jl +++ /dev/null @@ -1,71 +0,0 @@ -module TestCOCG - -using IterativeSolvers -using Test -using Random -using LinearAlgebra - -@testset ("Conjugate Orthogonal Conjugate Gradient") begin - -Random.seed!(1234321) -n = 20 - -@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64) - A = rand(T, n, n) - A = A + transpose(A) + 15I - x = ones(T, n) - b = A * x - - reltol = √eps(real(T)) - - # Solve without preconditioner - x1, his1 = cocg(A, b, reltol = reltol, maxiter = 100, log = true) - @test isa(his1, ConvergenceHistory) - @test norm(A * x1 - b) / norm(b) ≤ reltol - - # With an initial guess - x_guess = rand(T, n) - x2, his2 = cocg!(x_guess, A, b, reltol = reltol, maxiter = 100, log = true) - @test isa(his2, ConvergenceHistory) - @test x2 == x_guess - @test norm(A * x2 - b) / norm(b) ≤ reltol - - # The following tests fails CI on Windows and Ubuntu due to a - # `SingularException(4)` - if T == Float32 && (Sys.iswindows() || Sys.islinux()) - continue - end - # Do an exact LU decomp of a nearby matrix - F = lu(A + rand(T, n, n)) - x3, his3 = cocg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true) - @test norm(A * x3 - b) / norm(b) ≤ reltol -end - -@testset "Termination criterion" begin - for T in (Float32, Float64, ComplexF32, ComplexF64) - A = T[ 2 -1 0 - -1 2 -1 - 0 -1 2] - n = size(A, 2) - b = ones(T, n) - x0 = A \ b - perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n] - - # If the initial residual is small and a small relative tolerance is used, - # many iterations are necessary - x = x0 + perturbation - initial_residual = norm(A * x - b) - x, ch = cocg!(x, A, b, log=true) - @test 2 ≤ niters(ch) ≤ n - - # If the initial residual is small and a large absolute tolerance is used, - # no iterations are necessary - x = x0 + perturbation - initial_residual = norm(A * x - b) - x, ch = cocg!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true) - @test niters(ch) == 0 - end -end -end - -end # module diff --git a/test/runtests.jl b/test/runtests.jl index 01e6e582..e612a53a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,9 +12,6 @@ include("stationary.jl") #Conjugate gradients include("cg.jl") -#Conjugate Orthogonal Conjugate gradient -include("cocg.jl") - #BiCGStab(l) include("bicgstabl.jl") From 33378848ba021569162ed71562939527ca0319b8 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Tue, 26 Jan 2021 08:29:19 -0500 Subject: [PATCH 03/15] Add more unit tests for COCG --- src/cg.jl | 2 +- test/cg.jl | 20 ++++++++++++++++++-- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index f95d087f..8986c663 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -3,7 +3,7 @@ using Printf export cg, cg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables # Conjugated dot product -_dot(x, ::Val{true}) = sum(abs2(xₖ) for xₖ in x) +_dot(x, ::Val{true}) = sum(abs2, x) # for x::Complex, returns Real _dot(x, y, ::Val{true}) = dot(x, y) # Unconjugated dot product diff --git a/test/cg.jl b/test/cg.jl index 1cbdfd5e..04666e2e 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -21,10 +21,26 @@ ldiv!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal Random.seed!(1234321) +@testset "Vector{$T}, conjugated and unconjugated dot products" for T in (ComplexF32, ComplexF64) + n = 100 + x = rand(T, n) + y = rand(T, n) + + # Conjugated dot product + @test IterativeSolvers._dot(x, Val(true)) ≈ x'x + @test IterativeSolvers._dot(x, y, Val(true)) ≈ x'y + @test IterativeSolvers._dot(x, Val(true)) ≈ IterativeSolvers._dot(x, x, Val(true)) + + # Unonjugated dot product + @test IterativeSolvers._dot(x, Val(false)) ≈ transpose(x) * x + @test IterativeSolvers._dot(x, y, Val(false)) ≈ transpose(x) * y + @test IterativeSolvers._dot(x, Val(false)) ≈ IterativeSolvers._dot(x, x, Val(false)) +end + @testset "Small full system" begin n = 10 - @testset "Matrix{$T}, Conjugated Dot Product" for T in (Float32, Float64, ComplexF32, ComplexF64) + @testset "Matrix{$T}, conjugated dot product" for T in (Float32, Float64, ComplexF32, ComplexF64) A = rand(T, n, n) A = A' * A + I b = rand(T, n) @@ -51,7 +67,7 @@ Random.seed!(1234321) @test x0 == zeros(T, n) end - @testset "Matrix{$T}, Unconjugated Dot Product" for T in (Float32, Float64, ComplexF32, ComplexF64) + @testset "Matrix{$T}, unconjugated dot product" for T in (Float32, Float64, ComplexF32, ComplexF64) A = rand(T, n, n) A = A + transpose(A) + 15I x = ones(T, n) From f7df543de8b0de37b8bca9eface4f508b3554279 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 22 Feb 2021 11:01:04 -0500 Subject: [PATCH 04/15] Move definition of unconjugated dot product to common.jl The usage of the unconjugated dot product is not limited to implementing COCG from CG, so it is more appropriate to move its definition outside cg.jl. For example, the unconjugated dot product can be used in implementing COCR from CR. --- src/cg.jl | 8 -------- src/common.jl | 11 +++++++++++ test/cg.jl | 16 ---------------- test/common.jl | 16 ++++++++++++++++ 4 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index 8986c663..d8a0c11b 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -2,14 +2,6 @@ import Base: iterate using Printf export cg, cg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables -# Conjugated dot product -_dot(x, ::Val{true}) = sum(abs2, x) # for x::Complex, returns Real -_dot(x, y, ::Val{true}) = dot(x, y) - -# Unconjugated dot product -_dot(x, ::Val{false}) = sum(xₖ^2 for xₖ in x) -_dot(x, y, ::Val{false}) = sum(prod, zip(x,y)) - mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, boolT <: Union{Val{true},Val{false}}} A::matT x::solT diff --git a/src/common.jl b/src/common.jl index 165c841e..9129ab5a 100644 --- a/src/common.jl +++ b/src/common.jl @@ -30,3 +30,14 @@ struct Identity end \(::Identity, x) = copy(x) ldiv!(::Identity, x) = x ldiv!(y, ::Identity, x) = copyto!(y, x) + +""" +Conjugated and unconjugated dot products +""" +# Conjugated dot product +_dot(x, ::Val{true}) = sum(abs2, x) # for x::Complex, returns Real +_dot(x, y, ::Val{true}) = dot(x, y) + +# Unconjugated dot product +_dot(x, ::Val{false}) = sum(xₖ^2 for xₖ in x) +_dot(x, y, ::Val{false}) = sum(prod, zip(x,y)) diff --git a/test/cg.jl b/test/cg.jl index 04666e2e..c882e4a7 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -21,22 +21,6 @@ ldiv!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal Random.seed!(1234321) -@testset "Vector{$T}, conjugated and unconjugated dot products" for T in (ComplexF32, ComplexF64) - n = 100 - x = rand(T, n) - y = rand(T, n) - - # Conjugated dot product - @test IterativeSolvers._dot(x, Val(true)) ≈ x'x - @test IterativeSolvers._dot(x, y, Val(true)) ≈ x'y - @test IterativeSolvers._dot(x, Val(true)) ≈ IterativeSolvers._dot(x, x, Val(true)) - - # Unonjugated dot product - @test IterativeSolvers._dot(x, Val(false)) ≈ transpose(x) * x - @test IterativeSolvers._dot(x, y, Val(false)) ≈ transpose(x) * y - @test IterativeSolvers._dot(x, Val(false)) ≈ IterativeSolvers._dot(x, x, Val(false)) -end - @testset "Small full system" begin n = 10 diff --git a/test/common.jl b/test/common.jl index 7e08ff69..2f51d9a0 100644 --- a/test/common.jl +++ b/test/common.jl @@ -27,6 +27,22 @@ end @test ldiv!(y, P, copy(x)) == x end +@testset "Vector{$T}, conjugated and unconjugated dot products" for T in (ComplexF32, ComplexF64) + n = 100 + x = rand(T, n) + y = rand(T, n) + + # Conjugated dot product + @test IterativeSolvers._dot(x, Val(true)) ≈ x'x + @test IterativeSolvers._dot(x, y, Val(true)) ≈ x'y + @test IterativeSolvers._dot(x, Val(true)) ≈ IterativeSolvers._dot(x, x, Val(true)) + + # Unonjugated dot product + @test IterativeSolvers._dot(x, Val(false)) ≈ transpose(x) * x + @test IterativeSolvers._dot(x, y, Val(false)) ≈ transpose(x) * y + @test IterativeSolvers._dot(x, Val(false)) ≈ IterativeSolvers._dot(x, x, Val(false)) +end + end DocMeta.setdocmeta!(IterativeSolvers, :DocTestSetup, :(using IterativeSolvers); recursive=true) From d5b31f42d4c9f1df57b9ecff96027dc58d9d8781 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 6 Mar 2021 16:39:42 -0500 Subject: [PATCH 05/15] Define cocg() using DotType in cg() --- src/cg.jl | 48 ++++++++++++++++++++++++++++++++---------------- src/common.jl | 16 +++++++++------- test/cg.jl | 6 +++--- test/common.jl | 10 ++++------ 4 files changed, 48 insertions(+), 32 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index d8a0c11b..36c34a0f 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -1,8 +1,8 @@ import Base: iterate using Printf -export cg, cg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables +export cg, cg!, cocg, cocg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables -mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, boolT <: Union{Val{true},Val{false}}} +mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: DotType} A::matT x::solT r::vecT @@ -13,10 +13,10 @@ mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, bool ρ_prev::paramT maxiter::Int mv_products::Int - conjugate_dot::boolT + dot_type::dotT end -mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, boolT <: Union{Val{true},Val{false}}} +mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: DotType} Pl::precT A::matT x::solT @@ -28,7 +28,7 @@ mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Numb ρ_prev::paramT maxiter::Int mv_products::Int - conjugate_dot::boolT + dot_type::dotT end @inline converged(it::Union{CGIterable, PCGIterable}) = it.residual ≤ it.tol @@ -49,13 +49,14 @@ function iterate(it::CGIterable, iteration::Int=start(it)) end # u := r + βu (almost an axpy) - ρ = _dot(it.r, it.conjugate_dot) + ρ = isa(it.dot_type, ConjugatedDot) ? it.residual^2 : _norm(it.r, it.dot_type)^2 β = ρ / it.ρ_prev + it.u .= it.r .+ β .* it.u # c = A * u mul!(it.c, it.A, it.u) - α = ρ / _dot(it.u, it.c, it.conjugate_dot) + α = ρ / _dot(it.u, it.c, it.dot_type) # Improve solution and residual it.ρ_prev = ρ @@ -82,13 +83,13 @@ function iterate(it::PCGIterable, iteration::Int=start(it)) ldiv!(it.c, it.Pl, it.r) # u := c + βu (almost an axpy) - ρ = _dot(it.r, it.c, it.conjugate_dot) + ρ = _dot(it.r, it.c, it.dot_type) β = ρ / it.ρ_prev it.u .= it.c .+ β .* it.u # c = A * u mul!(it.c, it.A, it.u) - α = ρ / _dot(it.u, it.c, it.conjugate_dot) + α = ρ / _dot(it.u, it.c, it.dot_type) # Improve solution and residual it.ρ_prev = ρ @@ -125,7 +126,7 @@ function cg_iterator!(x, A, b, Pl = Identity(); maxiter::Int = size(A, 2), statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), initially_zero::Bool = false, - conjugate_dot::Bool = true) + dot_type::DotType = ConjugatedDot()) u = statevars.u r = statevars.r c = statevars.c @@ -145,13 +146,13 @@ function cg_iterator!(x, A, b, Pl = Identity(); # Return the iterable if isa(Pl, Identity) - return CGIterable(A, x, r, c, u, tolerance, residual, - conjugate_dot ? one(real(eltype(r))) : one(eltype(r)), # for conjugated dot, ρ_prev remains real - maxiter, mv_products, Val(conjugate_dot)) + return CGIterable(A, x, r, c, u, + tolerance, residual, one(eltype(r)), + maxiter, mv_products, dot_type) else return PCGIterable(Pl, A, x, r, c, u, tolerance, residual, one(eltype(r)), - maxiter, mv_products, Val(conjugate_dot)) + maxiter, mv_products, dot_type) end end @@ -212,7 +213,7 @@ function cg!(x, A, b; statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), verbose::Bool = false, Pl = Identity(), - conjugate_dot::Bool = true, + dot_type::DotType = ConjugatedDot(), kwargs...) history = ConvergenceHistory(partial = !log) history[:abstol] = abstol @@ -221,7 +222,7 @@ function cg!(x, A, b; # Actually perform CG iterable = cg_iterator!(x, A, b, Pl; abstol = abstol, reltol = reltol, maxiter = maxiter, - statevars = statevars, conjugate_dot = conjugate_dot, kwargs...) + statevars = statevars, dot_type = dot_type, kwargs...) if log history.mvps = iterable.mv_products end @@ -239,3 +240,18 @@ function cg!(x, A, b; log ? (iterable.x, history) : iterable.x end + +""" + cocg(A, b; kwargs...) -> x, [history] + +Same as [`cocg!`](@ref), but allocates a solution vector `x` initialized with zeros. +""" +cocg(A, b; kwargs...) = cocg!(zerox(A, b), A, b; initially_zero = true, kwargs...) + +""" + cocg!(x, A, b; kwargs...) -> x, [history] + +Same as [`cg!`](@ref), but uses the unconjugated dot product instead of the usual, +conjugated dot product. +""" +cocg!(x, A, b; kwargs...) = cg!(x, A, b; dot_type = UnconjugatedDot(), kwargs...) diff --git a/src/common.jl b/src/common.jl index 9129ab5a..a7e1d1aa 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,6 +1,6 @@ import LinearAlgebra: ldiv!, \ -export Identity +export Identity, ConjugatedDot, UnconjugatedDot #### Type-handling """ @@ -34,10 +34,12 @@ ldiv!(y, ::Identity, x) = copyto!(y, x) """ Conjugated and unconjugated dot products """ -# Conjugated dot product -_dot(x, ::Val{true}) = sum(abs2, x) # for x::Complex, returns Real -_dot(x, y, ::Val{true}) = dot(x, y) +abstract type DotType end +struct ConjugatedDot <: DotType end +struct UnconjugatedDot <: DotType end -# Unconjugated dot product -_dot(x, ::Val{false}) = sum(xₖ^2 for xₖ in x) -_dot(x, y, ::Val{false}) = sum(prod, zip(x,y)) +_norm(x, ::ConjugatedDot) = norm(x) +_dot(x, y, ::ConjugatedDot) = dot(x, y) + +_norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ^2 for xₖ in x)) +_dot(x, y, ::UnconjugatedDot) = sum(prod, zip(x, y)) diff --git a/test/cg.jl b/test/cg.jl index c882e4a7..acc7e882 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -60,13 +60,13 @@ Random.seed!(1234321) reltol = √eps(real(T)) # Solve without preconditioner - x1, his1 = cg(A, b, reltol = reltol, maxiter = 100, log = true, conjugate_dot = false) + x1, his1 = cocg(A, b, reltol = reltol, maxiter = 100, log = true) @test isa(his1, ConvergenceHistory) @test norm(A * x1 - b) / norm(b) ≤ reltol # With an initial guess x_guess = rand(T, n) - x2, his2 = cg!(x_guess, A, b, reltol = reltol, maxiter = 100, log = true, conjugate_dot = false) + x2, his2 = cocg!(x_guess, A, b, reltol = reltol, maxiter = 100, log = true) @test isa(his2, ConvergenceHistory) @test x2 == x_guess @test norm(A * x2 - b) / norm(b) ≤ reltol @@ -78,7 +78,7 @@ Random.seed!(1234321) end # Do an exact LU decomp of a nearby matrix F = lu(A + rand(T, n, n)) - x3, his3 = cg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true, conjugate_dot = false) + x3, his3 = cocg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true) @test norm(A * x3 - b) / norm(b) ≤ reltol end end diff --git a/test/common.jl b/test/common.jl index 2f51d9a0..36b10860 100644 --- a/test/common.jl +++ b/test/common.jl @@ -33,14 +33,12 @@ end y = rand(T, n) # Conjugated dot product - @test IterativeSolvers._dot(x, Val(true)) ≈ x'x - @test IterativeSolvers._dot(x, y, Val(true)) ≈ x'y - @test IterativeSolvers._dot(x, Val(true)) ≈ IterativeSolvers._dot(x, x, Val(true)) + @test IterativeSolvers._norm(x, ConjugatedDot()) ≈ sqrt(x'x) + @test IterativeSolvers._dot(x, y, ConjugatedDot()) ≈ x'y # Unonjugated dot product - @test IterativeSolvers._dot(x, Val(false)) ≈ transpose(x) * x - @test IterativeSolvers._dot(x, y, Val(false)) ≈ transpose(x) * y - @test IterativeSolvers._dot(x, Val(false)) ≈ IterativeSolvers._dot(x, x, Val(false)) + @test IterativeSolvers._norm(x, UnconjugatedDot()) ≈ sqrt(transpose(x) * x) + @test IterativeSolvers._dot(x, y, UnconjugatedDot()) ≈ transpose(x) * y end end From 97315be26fc501e1e8c6b6ce814a6f0b3702661a Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Wed, 10 Mar 2021 19:12:03 -0500 Subject: [PATCH 06/15] Rename DotType AbstractDot --- src/cg.jl | 8 ++++---- src/common.jl | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index 36c34a0f..9559e98c 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -2,7 +2,7 @@ import Base: iterate using Printf export cg, cg!, cocg, cocg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables -mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: DotType} +mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: AbstractDot} A::matT x::solT r::vecT @@ -16,7 +16,7 @@ mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, dotT dot_type::dotT end -mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: DotType} +mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: AbstractDot} Pl::precT A::matT x::solT @@ -126,7 +126,7 @@ function cg_iterator!(x, A, b, Pl = Identity(); maxiter::Int = size(A, 2), statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), initially_zero::Bool = false, - dot_type::DotType = ConjugatedDot()) + dot_type::AbstractDot = ConjugatedDot()) u = statevars.u r = statevars.r c = statevars.c @@ -213,7 +213,7 @@ function cg!(x, A, b; statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), verbose::Bool = false, Pl = Identity(), - dot_type::DotType = ConjugatedDot(), + dot_type::AbstractDot = ConjugatedDot(), kwargs...) history = ConvergenceHistory(partial = !log) history[:abstol] = abstol diff --git a/src/common.jl b/src/common.jl index a7e1d1aa..62868d7f 100644 --- a/src/common.jl +++ b/src/common.jl @@ -34,9 +34,9 @@ ldiv!(y, ::Identity, x) = copyto!(y, x) """ Conjugated and unconjugated dot products """ -abstract type DotType end -struct ConjugatedDot <: DotType end -struct UnconjugatedDot <: DotType end +abstract type AbstractDot end +struct ConjugatedDot <: AbstractDot end +struct UnconjugatedDot <: AbstractDot end _norm(x, ::ConjugatedDot) = norm(x) _dot(x, y, ::ConjugatedDot) = dot(x, y) From 593e88cbe37d97331d4d379f3a73568072aaa342 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Wed, 10 Mar 2021 22:38:48 -0500 Subject: [PATCH 07/15] Rename dot_type dotproduct --- src/cg.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index 9559e98c..fce607bc 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -13,7 +13,7 @@ mutable struct CGIterable{matT, solT, vecT, numT <: Real, paramT <: Number, dotT ρ_prev::paramT maxiter::Int mv_products::Int - dot_type::dotT + dotproduct::dotT end mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number, dotT <: AbstractDot} @@ -28,7 +28,7 @@ mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Numb ρ_prev::paramT maxiter::Int mv_products::Int - dot_type::dotT + dotproduct::dotT end @inline converged(it::Union{CGIterable, PCGIterable}) = it.residual ≤ it.tol @@ -49,14 +49,14 @@ function iterate(it::CGIterable, iteration::Int=start(it)) end # u := r + βu (almost an axpy) - ρ = isa(it.dot_type, ConjugatedDot) ? it.residual^2 : _norm(it.r, it.dot_type)^2 + ρ = isa(it.dotproduct, ConjugatedDot) ? it.residual^2 : _norm(it.r, it.dotproduct)^2 β = ρ / it.ρ_prev it.u .= it.r .+ β .* it.u # c = A * u mul!(it.c, it.A, it.u) - α = ρ / _dot(it.u, it.c, it.dot_type) + α = ρ / _dot(it.u, it.c, it.dotproduct) # Improve solution and residual it.ρ_prev = ρ @@ -83,13 +83,13 @@ function iterate(it::PCGIterable, iteration::Int=start(it)) ldiv!(it.c, it.Pl, it.r) # u := c + βu (almost an axpy) - ρ = _dot(it.r, it.c, it.dot_type) + ρ = _dot(it.r, it.c, it.dotproduct) β = ρ / it.ρ_prev it.u .= it.c .+ β .* it.u # c = A * u mul!(it.c, it.A, it.u) - α = ρ / _dot(it.u, it.c, it.dot_type) + α = ρ / _dot(it.u, it.c, it.dotproduct) # Improve solution and residual it.ρ_prev = ρ @@ -126,7 +126,7 @@ function cg_iterator!(x, A, b, Pl = Identity(); maxiter::Int = size(A, 2), statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), initially_zero::Bool = false, - dot_type::AbstractDot = ConjugatedDot()) + dotproduct::AbstractDot = ConjugatedDot()) u = statevars.u r = statevars.r c = statevars.c @@ -148,11 +148,11 @@ function cg_iterator!(x, A, b, Pl = Identity(); if isa(Pl, Identity) return CGIterable(A, x, r, c, u, tolerance, residual, one(eltype(r)), - maxiter, mv_products, dot_type) + maxiter, mv_products, dotproduct) else return PCGIterable(Pl, A, x, r, c, u, tolerance, residual, one(eltype(r)), - maxiter, mv_products, dot_type) + maxiter, mv_products, dotproduct) end end @@ -213,7 +213,7 @@ function cg!(x, A, b; statevars::CGStateVariables = CGStateVariables(zero(x), similar(x), similar(x)), verbose::Bool = false, Pl = Identity(), - dot_type::AbstractDot = ConjugatedDot(), + dotproduct::AbstractDot = ConjugatedDot(), kwargs...) history = ConvergenceHistory(partial = !log) history[:abstol] = abstol @@ -222,7 +222,7 @@ function cg!(x, A, b; # Actually perform CG iterable = cg_iterator!(x, A, b, Pl; abstol = abstol, reltol = reltol, maxiter = maxiter, - statevars = statevars, dot_type = dot_type, kwargs...) + statevars = statevars, dotproduct = dotproduct, kwargs...) if log history.mvps = iterable.mv_products end @@ -254,4 +254,4 @@ cocg(A, b; kwargs...) = cocg!(zerox(A, b), A, b; initially_zero = true, kwargs.. Same as [`cg!`](@ref), but uses the unconjugated dot product instead of the usual, conjugated dot product. """ -cocg!(x, A, b; kwargs...) = cg!(x, A, b; dot_type = UnconjugatedDot(), kwargs...) +cocg!(x, A, b; kwargs...) = cg!(x, A, b; dotproduct = UnconjugatedDot(), kwargs...) From 835a894d5cd0c7c4a54d312065c3f66bdeb13c7b Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 14 Mar 2021 15:34:02 -0400 Subject: [PATCH 08/15] Improve efficiency of unconjugated norm and dot --- src/common.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/common.jl b/src/common.jl index 62868d7f..b0acaecb 100644 --- a/src/common.jl +++ b/src/common.jl @@ -41,5 +41,5 @@ struct UnconjugatedDot <: AbstractDot end _norm(x, ::ConjugatedDot) = norm(x) _dot(x, y, ::ConjugatedDot) = dot(x, y) -_norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ^2 for xₖ in x)) -_dot(x, y, ::UnconjugatedDot) = sum(prod, zip(x, y)) +_norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ->xₖ^2, x)) +_dot(x, y, ::UnconjugatedDot) = transpose(@view(x[:])) * @view(y[:]) # allocating, but faster than sum(prod, zip(x,y)) From 9c47e7fef0d17a109d1ffec9031a210ab0df41f6 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Wed, 17 Mar 2021 23:05:16 -0400 Subject: [PATCH 09/15] =?UTF-8?q?Utilize=20=E2=89=88=20in=20test/cg.jl?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/cg.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/cg.jl b/test/cg.jl index acc7e882..640ac35e 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -32,7 +32,7 @@ Random.seed!(1234321) x,ch = cg(A, b; reltol=reltol, maxiter=2n, log=true) @test isa(ch, ConvergenceHistory) - @test norm(A*x - b) / norm(b) ≤ reltol + @test A*x ≈ b rtol=reltol @test ch.isconverged # If you start from the exact solution, you should converge immediately @@ -62,14 +62,14 @@ Random.seed!(1234321) # Solve without preconditioner x1, his1 = cocg(A, b, reltol = reltol, maxiter = 100, log = true) @test isa(his1, ConvergenceHistory) - @test norm(A * x1 - b) / norm(b) ≤ reltol + @test A*x1 ≈ b rtol=reltol # With an initial guess x_guess = rand(T, n) x2, his2 = cocg!(x_guess, A, b, reltol = reltol, maxiter = 100, log = true) @test isa(his2, ConvergenceHistory) @test x2 == x_guess - @test norm(A * x2 - b) / norm(b) ≤ reltol + @test A*x2 ≈ b rtol=reltol # The following tests fails CI on Windows and Ubuntu due to a # `SingularException(4)` @@ -79,7 +79,7 @@ Random.seed!(1234321) # Do an exact LU decomp of a nearby matrix F = lu(A + rand(T, n, n)) x3, his3 = cocg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true) - @test norm(A * x3 - b) / norm(b) ≤ reltol + @test A*x3 ≈ b rtol=reltol end end @@ -95,24 +95,24 @@ end @testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, Float32), Ti in (Int64, Int32) xCG = cg(A, rhs; reltol=reltol, maxiter=100) xJAC = cg(A, rhs; Pl=P, reltol=reltol, maxiter=100) - @test norm(A * xCG - rhs) ≤ reltol - @test norm(A * xJAC - rhs) ≤ reltol + @test A*xCG ≈ rhs rtol=reltol + @test A*xJAC ≈ rhs rtol=reltol end Af = LinearMap(A) @testset "Function" begin xCG = cg(Af, rhs; reltol=reltol, maxiter=100) xJAC = cg(Af, rhs; Pl=P, reltol=reltol, maxiter=100) - @test norm(A * xCG - rhs) ≤ reltol - @test norm(A * xJAC - rhs) ≤ reltol + @test A*xCG ≈ rhs rtol=reltol + @test A*xJAC ≈ rhs rtol=reltol end @testset "Function with specified starting guess" begin x0 = randn(size(rhs)) xCG, hCG = cg!(copy(x0), Af, rhs; abstol=abstol, reltol=0.0, maxiter=100, log=true) xJAC, hJAC = cg!(copy(x0), Af, rhs; Pl=P, abstol=abstol, reltol=0.0, maxiter=100, log=true) - @test norm(A * xCG - rhs) ≤ reltol - @test norm(A * xJAC - rhs) ≤ reltol + @test A*xCG ≈ rhs rtol=reltol + @test A*xJAC ≈ rhs rtol=reltol @test niters(hJAC) == niters(hCG) end end From 3cd796971aa9fc30f00ba19b51cd2674c3adf14f Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Wed, 17 Mar 2021 23:06:50 -0400 Subject: [PATCH 10/15] Remove unnecessary if ...; continue; end block from test/cg.jl --- test/cg.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/test/cg.jl b/test/cg.jl index 640ac35e..34c7121d 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -71,11 +71,6 @@ Random.seed!(1234321) @test x2 == x_guess @test A*x2 ≈ b rtol=reltol - # The following tests fails CI on Windows and Ubuntu due to a - # `SingularException(4)` - if T == Float32 && (Sys.iswindows() || Sys.islinux()) - continue - end # Do an exact LU decomp of a nearby matrix F = lu(A + rand(T, n, n)) x3, his3 = cocg(A, b, Pl = F, maxiter = 100, reltol = reltol, log = true) From d16fecba2030adc5ce370698d63efa7045740d74 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 21 Mar 2021 19:22:18 -0400 Subject: [PATCH 11/15] Update docstring for cocg!() --- src/cg.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cg.jl b/src/cg.jl index fce607bc..28934e5e 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -251,7 +251,8 @@ cocg(A, b; kwargs...) = cocg!(zerox(A, b), A, b; initially_zero = true, kwargs.. """ cocg!(x, A, b; kwargs...) -> x, [history] -Same as [`cg!`](@ref), but uses the unconjugated dot product instead of the usual, -conjugated dot product. +Same as [`cg!`](@ref), but uses the unconjugated dot product (`xᵀy`) instead of the usual, +conjugated dot product (`x'y`) in the algorithm. It is for solving linear systems with +matrices `A` that are complex-symmetric (`Aᵀ == A`) rahter than Hermitian (`A' == A`). """ cocg!(x, A, b; kwargs...) = cg!(x, A, b; dotproduct = UnconjugatedDot(), kwargs...) From 04c3c16eaac1651f0af7ffb3b2b103a1d6253b69 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 12 Apr 2021 12:26:17 -0400 Subject: [PATCH 12/15] Fix typo in comment --- src/cg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cg.jl b/src/cg.jl index 28934e5e..dcea8d69 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -253,6 +253,6 @@ cocg(A, b; kwargs...) = cocg!(zerox(A, b), A, b; initially_zero = true, kwargs.. Same as [`cg!`](@ref), but uses the unconjugated dot product (`xᵀy`) instead of the usual, conjugated dot product (`x'y`) in the algorithm. It is for solving linear systems with -matrices `A` that are complex-symmetric (`Aᵀ == A`) rahter than Hermitian (`A' == A`). +matrices `A` that are complex-symmetric (`Aᵀ == A`) rather than Hermitian (`A' == A`). """ cocg!(x, A, b; kwargs...) = cg!(x, A, b; dotproduct = UnconjugatedDot(), kwargs...) From 48001292f565159ae204e74a222799122816a990 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 12 Apr 2021 12:27:10 -0400 Subject: [PATCH 13/15] Remove unnecessary comments --- src/common.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common.jl b/src/common.jl index b0acaecb..0035b23b 100644 --- a/src/common.jl +++ b/src/common.jl @@ -42,4 +42,4 @@ _norm(x, ::ConjugatedDot) = norm(x) _dot(x, y, ::ConjugatedDot) = dot(x, y) _norm(x, ::UnconjugatedDot) = sqrt(sum(xₖ->xₖ^2, x)) -_dot(x, y, ::UnconjugatedDot) = transpose(@view(x[:])) * @view(y[:]) # allocating, but faster than sum(prod, zip(x,y)) +_dot(x, y, ::UnconjugatedDot) = transpose(@view(x[:])) * @view(y[:]) From f3710c35c8436dcc5838aaf74405268a5509ea23 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Tue, 4 Jan 2022 21:22:42 -0500 Subject: [PATCH 14/15] Use Conjugate Gradient (no plural) as full name of CG --- docs/make.jl | 2 +- docs/src/index.md | 4 ++-- docs/src/linear_systems/cg.md | 4 ++-- docs/src/linear_systems/chebyshev.md | 2 +- test/cg.jl | 2 +- test/runtests.jl | 2 +- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index a8111f0e..4eefccd2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,7 +16,7 @@ makedocs( "Getting started" => "getting_started.md", "Preconditioning" => "preconditioning.md", "Linear systems" => [ - "Conjugate Gradients" => "linear_systems/cg.md", + "Conjugate Gradient" => "linear_systems/cg.md", "Chebyshev iteration" => "linear_systems/chebyshev.md", "MINRES" => "linear_systems/minres.md", "BiCGStab(l)" => "linear_systems/bicgstabl.md", diff --git a/docs/src/index.md b/docs/src/index.md index 2881c612..f3cbfa34 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -12,13 +12,13 @@ When solving linear systems $Ax = b$ for a square matrix $A$ there are quite som | Method | When to use it | |---------------------|--------------------------------------------------------------------------| -| [Conjugate Gradients](@ref CG) | Best choice for **symmetric**, **positive-definite** matrices | +| [Conjugate Gradient](@ref CG) | Best choice for **symmetric**, **positive-definite** matrices | | [MINRES](@ref MINRES) | For **symmetric**, **indefinite** matrices | | [GMRES](@ref GMRES) | For **nonsymmetric** matrices when a good [preconditioner](@ref Preconditioning) is available | | [IDR(s)](@ref IDRs) | For **nonsymmetric**, **strongly indefinite** problems without a good preconditioner | | [BiCGStab(l)](@ref BiCGStabl) | Otherwise for **nonsymmetric** problems | -We also offer [Chebyshev iteration](@ref Chebyshev) as an alternative to Conjugate Gradients when bounds on the spectrum are known. +We also offer [Chebyshev iteration](@ref Chebyshev) as an alternative to Conjugate Gradient when bounds on the spectrum are known. Stationary methods like [Jacobi](@ref), [Gauss-Seidel](@ref), [SOR](@ref) and [SSOR](@ref) can be used as smoothers to reduce high-frequency components in the error in just a few iterations. diff --git a/docs/src/linear_systems/cg.md b/docs/src/linear_systems/cg.md index 7fde0a13..8b688272 100644 --- a/docs/src/linear_systems/cg.md +++ b/docs/src/linear_systems/cg.md @@ -1,6 +1,6 @@ -# [Conjugate Gradients (CG)](@id CG) +# [Conjugate Gradient (CG)](@id CG) -Conjugate Gradients solves $Ax = b$ approximately for $x$ where $A$ is a symmetric, positive-definite linear operator and $b$ the right-hand side vector. The method uses short recurrences and therefore has fixed memory costs and fixed computational costs per iteration. +Conjugate Gradient solves $Ax = b$ approximately for $x$ where $A$ is a symmetric, positive-definite linear operator and $b$ the right-hand side vector. The method uses short recurrences and therefore has fixed memory costs and fixed computational costs per iteration. ## Usage diff --git a/docs/src/linear_systems/chebyshev.md b/docs/src/linear_systems/chebyshev.md index 4b55c918..fdc4d825 100644 --- a/docs/src/linear_systems/chebyshev.md +++ b/docs/src/linear_systems/chebyshev.md @@ -2,7 +2,7 @@ Chebyshev iteration solves the problem $Ax=b$ approximately for $x$ where $A$ is a symmetric, definite linear operator and $b$ the right-hand side vector. The methods assumes the interval $[\lambda_{min}, \lambda_{max}]$ containing all eigenvalues of $A$ is known, so that $x$ can be iteratively constructed via a Chebyshev polynomial with zeros in this interval. This polynomial ultimately acts as a filter that removes components in the direction of the eigenvectors from the initial residual. -The main advantage with respect to Conjugate Gradients is that BLAS1 operations such as inner products are avoided. +The main advantage with respect to Conjugate Gradient is that BLAS1 operations such as inner products are avoided. ## Usage diff --git a/test/cg.jl b/test/cg.jl index 34c7121d..a77857a0 100644 --- a/test/cg.jl +++ b/test/cg.jl @@ -17,7 +17,7 @@ end ldiv!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal -@testset "Conjugate Gradients" begin +@testset "Conjugate Gradient" begin Random.seed!(1234321) diff --git a/test/runtests.jl b/test/runtests.jl index e612a53a..d7830d46 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,7 @@ include("hessenberg.jl") #Stationary solvers include("stationary.jl") -#Conjugate gradients +#Conjugate gradient include("cg.jl") #BiCGStab(l) From b457247dd6098fff21b50ee4f7738e0cb61db2dc Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Tue, 4 Jan 2022 22:14:23 -0500 Subject: [PATCH 15/15] Add COCG in documentation --- docs/src/linear_systems/cg.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/src/linear_systems/cg.md b/docs/src/linear_systems/cg.md index 8b688272..6ce25566 100644 --- a/docs/src/linear_systems/cg.md +++ b/docs/src/linear_systems/cg.md @@ -1,12 +1,16 @@ # [Conjugate Gradient (CG)](@id CG) -Conjugate Gradient solves $Ax = b$ approximately for $x$ where $A$ is a symmetric, positive-definite linear operator and $b$ the right-hand side vector. The method uses short recurrences and therefore has fixed memory costs and fixed computational costs per iteration. +CG solves $Ax = b$ approximately for $x$ where $A$ is a Hermitian, positive-definite linear operator and $b$ the right-hand side vector. The method uses short recurrences and therefore has fixed memory costs and fixed computational costs per iteration. + +A simple variant called Conjugate Orthogonal Conjugate Gradient (COCG) is obtained by replacing the dot product ($x^*y$) in the CG algorithm with the unconjugated dot product ($xᵀy$). For complex symmetric matrices ($A = Aᵀ$), COCG is equivalent to Biconjugate Gradient (BiCG) but takes only one matrix-vector multiplication per iteration rather than two of BiCG. Therefore, COCG converges twice as fast in time as BiCG. Also, like in BiCG, $A$ in COCG does not need to be positive-definite. ## Usage ```@docs cg cg! +cocg +cocg! ``` ## On the GPU