From 0aa0760d73812d9de504068ee48e131f666cc142 Mon Sep 17 00:00:00 2001 From: Harmen Stoppels Date: Mon, 28 Aug 2017 16:33:24 +0200 Subject: [PATCH] Fix element types (#163) * Remove unused methods * Clean up common.jl * Ensure r has the type of x, not b. Remove b from the iterable since it is not used anyway * Make residual type in BiCGStab(l) equal to solution type and remove rhs from iterable * Chebyshev: residual type should equal x's type and remove b from iterable * Similar story for GMRES * Use x element type everywhere * Use different types for solution, temporary and rhs in iterables of stationary methods * zerox over zeros --- src/bicgstabl.jl | 11 +++++----- src/cg.jl | 26 ++++++++++------------- src/chebyshev.jl | 15 ++++++------- src/common.jl | 41 ++--------------------------------- src/gmres.jl | 18 ++++++++-------- src/minres.jl | 19 +++++++++-------- src/stationary.jl | 24 ++++++++++----------- src/stationary_sparse.jl | 46 +++++++++++++++++++++++----------------- 8 files changed, 82 insertions(+), 118 deletions(-) diff --git a/src/bicgstabl.jl b/src/bicgstabl.jl index 013f5bf4..0d4c8964 100644 --- a/src/bicgstabl.jl +++ b/src/bicgstabl.jl @@ -2,12 +2,11 @@ export bicgstabl, bicgstabl!, bicgstabl_iterator, bicgstabl_iterator!, BiCGStabI import Base: start, next, done -mutable struct BiCGStabIterable{precT, matT, vecT <: AbstractVector, smallMatT <: AbstractMatrix, realT <: Real, scalarT <: Number} +mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, smallMatT <: AbstractMatrix, realT <: Real, scalarT <: Number} A::matT - b::vecT l::Int - x::vecT + x::solT r_shadow::vecT rs::smallMatT us::smallMatT @@ -33,7 +32,7 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2; initial_zero = false, tol = sqrt(eps(real(eltype(b)))) ) - T = eltype(b) + T = eltype(x) n = size(A, 1) mv_products = 0 @@ -69,7 +68,7 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2; # Stopping condition based on relative tolerance. reltol = nrm * tol - BiCGStabIterable(A, b, l, x, r_shadow, rs, us, + BiCGStabIterable(A, l, x, r_shadow, rs, us, max_mv_products, mv_products, reltol, nrm, Pl, γ, ω, σ, M @@ -81,7 +80,7 @@ end @inline done(it::BiCGStabIterable, iteration::Int) = it.mv_products ≥ it.max_mv_products || converged(it) function next(it::BiCGStabIterable, iteration::Int) - T = eltype(it.b) + T = eltype(it.x) L = 2 : it.l + 1 it.σ = -it.ω * it.σ diff --git a/src/cg.jl b/src/cg.jl index e8cfdfe3..19f02f3f 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -2,10 +2,9 @@ import Base: start, next, done export cg, cg!, CGIterable, PCGIterable, cg_iterator, cg_iterator! -mutable struct CGIterable{matT, vecT <: AbstractVector, numT <: Real} +mutable struct CGIterable{matT, solT, vecT, numT <: Real} A::matT - x::vecT - b::vecT + x::solT r::vecT c::vecT u::vecT @@ -16,11 +15,10 @@ mutable struct CGIterable{matT, vecT <: AbstractVector, numT <: Real} mv_products::Int end -mutable struct PCGIterable{precT, matT, vecT <: AbstractVector, numT <: Real, paramT <: Number} +mutable struct PCGIterable{precT, matT, solT, vecT, numT <: Real, paramT <: Number} Pl::precT A::matT - x::vecT - b::vecT + x::solT r::vecT c::vecT u::vecT @@ -46,7 +44,7 @@ function next(it::CGIterable, iteration::Int) # u := r + βu (almost an axpy) β = it.residual^2 / it.prev_residual^2 @blas! it.u *= β - @blas! it.u += one(eltype(it.b)) * it.r + @blas! it.u += one(eltype(it.u)) * it.r # c = A * u A_mul_B!(it.c, it.A, it.u) @@ -76,7 +74,7 @@ function next(it::PCGIterable, iteration::Int) # u := c + βu (almost an axpy) β = it.ρ / ρ_prev @blas! it.u *= β - @blas! it.u += one(eltype(it.b)) * it.c + @blas! it.u += one(eltype(it.u)) * it.c # c = A * u A_mul_B!(it.c, it.A, it.u) @@ -102,7 +100,8 @@ function cg_iterator!(x, A, b, Pl = Identity(); initially_zero::Bool = false ) u = zeros(x) - r = copy(b) + r = similar(x) + copy!(r, b) # Compute r with an MV-product or not. if initially_zero @@ -120,16 +119,13 @@ function cg_iterator!(x, A, b, Pl = Identity(); # Return the iterable if isa(Pl, Identity) - return CGIterable(A, x, b, - r, c, u, + return CGIterable(A, x, r, c, u, reltol, residual, one(residual), maxiter, mv_products ) else - ρ = one(eltype(r)) - return PCGIterable(Pl, A, x, b, - r, c, u, - reltol, residual, ρ, + return PCGIterable(Pl, A, x, r, c, u, + reltol, residual, one(eltype(x)), maxiter, mv_products ) end diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 07425ae0..83b8cb57 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -2,12 +2,11 @@ import Base: next, start, done export chebyshev, chebyshev! -mutable struct ChebyshevIterable{precT, matT, vecT, realT <: Real} +mutable struct ChebyshevIterable{precT, matT, solT, vecT, realT <: Real} Pl::precT A::matT - b::vecT - x::vecT + x::solT r::vecT u::vecT c::vecT @@ -28,7 +27,7 @@ start(::ChebyshevIterable) = 0 done(c::ChebyshevIterable, iteration::Int) = iteration ≥ c.maxiter || converged(c) function next(cheb::ChebyshevIterable, iteration::Int) - T = eltype(cheb.u) + T = eltype(cheb.x) solve!(cheb.c, cheb.Pl, cheb.r) @@ -64,8 +63,9 @@ function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real; λ_avg = (λmax + λmin) / 2 λ_diff = (λmax - λmin) / 2 - T = eltype(b) - r = copy(b) + T = eltype(x) + r = similar(x) + copy!(r, b) u = zeros(x) c = similar(x) @@ -82,8 +82,7 @@ function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real; mv_products = 1 end - ChebyshevIterable(Pl, A, b, - x, r, u, c, + ChebyshevIterable(Pl, A, x, r, u, c, zero(real(T)), λ_avg, λ_diff, resnorm, reltol, maxiter, mv_products diff --git a/src/common.jl b/src/common.jl index 50d78e9b..2a41e308 100644 --- a/src/common.jl +++ b/src/common.jl @@ -12,60 +12,23 @@ Determine type of the division of an element of `b` against an element of `A`: """ Adivtype(A, b) = typeof(one(eltype(b))/one(eltype(A))) -""" - Amultype(A, x) -Determine type of the multiplication of an element of `b` with an element of `A`: -`typeof(one(eltype(A))*one(eltype(x)))` -""" -Amultype(A, x) = typeof(one(eltype(A))*one(eltype(x))) - -""" - randx(A, b) -Build a random unitary vector `Vector{T}`, where `T` is `Adivtype(A,b)`. -""" -function randx(A, b) - T = Adivtype(A, b) - x = initrand!(Array(T, size(A, 2))) -end - """ zerox(A, b) Build a zeros vector `Vector{T}`, where `T` is `Adivtype(A,b)`. """ -function zerox(A, b) - T = Adivtype(A, b) - x = zeros(T, size(A, 2)) -end +zerox(A, b) = zeros(Adivtype(A, b), size(A, 2)) #### Numerics """ solve(A,b) -Solve `A\b` with a direct solver. When `A` is a function `A(b)` is dispatched instead. +Solve `A\\b` with a direct solver. When `A` is a function `A(b)` is dispatched instead. """ solve(A::Function,b) = A(b) - solve(A,b) = A\b - solve!(out::AbstractArray{T},A::Int,b::AbstractArray{T}) where {T} = scale!(out,b, 1/A) - solve!(out::AbstractArray{T},A,b::AbstractArray{T}) where {T} = A_ldiv_B!(out,A,b) solve!(out::AbstractArray{T},A::Function,b::AbstractArray{T}) where {T} = copy!(out,A(b)) -""" - initrand!(v) -Overwrite `v` with a random unitary vector of the same length. -""" -function initrand!(v::Vector) - _randn!(v) - nv = norm(v) - for i = 1:length(v) - v[i] /= nv - end - v -end -_randn!(v::Array{Float64}) = randn!(v) -_randn!(v) = copy!(v, randn(length(v))) - # Identity preconditioner struct Identity end diff --git a/src/gmres.jl b/src/gmres.jl index 22c866d2..bf9d9ff6 100644 --- a/src/gmres.jl +++ b/src/gmres.jl @@ -28,11 +28,11 @@ Residual(order, T::Type) = Residual{T, real(T)}( one(real(T)) ) -mutable struct GMRESIterable{preclT, precrT, vecT <: AbstractVector, arnoldiT <: ArnoldiDecomp, residualT <: Residual, resT <: Real} +mutable struct GMRESIterable{preclT, precrT, solT, rhsT, vecT, arnoldiT <: ArnoldiDecomp, residualT <: Residual, resT <: Real} Pl::preclT Pr::precrT - x::vecT - b::vecT + x::solT + b::rhsT Ax::vecT # Some room to work in. arnoldi::arnoldiT @@ -98,7 +98,7 @@ function next(g::GMRESIterable, iteration::Int) g.residual.current, iteration + 1 end -gmres_iterable(A, b; kwargs...) = gmres_iterable!(zeros(b), A, b; initially_zero = true, kwargs...) +gmres_iterable(A, b; kwargs...) = gmres_iterable!(zerox(A, b), A, b; initially_zero = true, kwargs...) function gmres_iterable!(x, A, b; Pl = Identity(), @@ -106,17 +106,17 @@ function gmres_iterable!(x, A, b; tol = sqrt(eps(real(eltype(b)))), restart::Int = min(20, length(b)), maxiter::Int = restart, - initially_zero = false + initially_zero::Bool = false ) - T = eltype(b) + T = eltype(x) # Approximate solution arnoldi = ArnoldiDecomp(A, restart, T) residual = Residual(restart, T) - mv_products = initially_zero == true ? 1 : 0 + mv_products = initially_zero ? 1 : 0 # Workspace vector to reduce the # allocs. - Ax = similar(b) + Ax = similar(x) residual.current = init!(arnoldi, x, b, Pl, Ax, initially_zero = initially_zero) init_residual!(residual, residual.current) @@ -133,7 +133,7 @@ end Same as [`gmres!`](@ref), but allocates a solution vector `x` initialized with zeros. """ -gmres(A, b; kwargs...) = gmres!(zeros(b), A, b; initially_zero = true, kwargs...) +gmres(A, b; kwargs...) = gmres!(zerox(A, b), A, b; initially_zero = true, kwargs...) """ gmres!(x, A, b; kwargs...) -> x, [history] diff --git a/src/minres.jl b/src/minres.jl index d9bf99af..0ff891ef 100644 --- a/src/minres.jl +++ b/src/minres.jl @@ -3,10 +3,10 @@ export minres_iterable, minres, minres! import Base.LinAlg: BLAS.axpy!, givensAlgorithm import Base: start, next, done -mutable struct MINRESIterable{matT, vecT <: DenseVector, smallVecT <: DenseVector, rotT <: Number, realT <: Real} +mutable struct MINRESIterable{matT, solT, vecT <: DenseVector, smallVecT <: DenseVector, rotT <: Number, realT <: Real} A::matT skew_hermitian::Bool - x::vecT + x::solT # Krylov basis vectors v_prev::vecT @@ -44,15 +44,16 @@ function minres_iterable!(x, A, b; tol = sqrt(eps(real(eltype(b)))), maxiter = size(A, 1) ) - T = eltype(b) + T = eltype(x) HessenbergT = skew_hermitian ? T : real(T) - v_prev = similar(b) - v_curr = copy(b) - v_next = similar(b) - w_prev = similar(b) - w_curr = similar(b) - w_next = similar(b) + v_prev = similar(x) + v_curr = similar(x) + copy!(v_curr, b) + v_next = similar(x) + w_prev = similar(x) + w_curr = similar(x) + w_next = similar(x) mv_products = 0 diff --git a/src/stationary.jl b/src/stationary.jl index 684d8699..ced0a6ba 100644 --- a/src/stationary.jl +++ b/src/stationary.jl @@ -35,11 +35,11 @@ function jacobi!(x, A::AbstractMatrix, b; maxiter::Int=10) x end -mutable struct DenseJacobiIterable{matT,vecT} +mutable struct DenseJacobiIterable{matT,vecT,solT,rhsT} A::matT - x::vecT + x::solT next::vecT - b::vecT + b::rhsT maxiter::Int end @@ -93,10 +93,10 @@ function gauss_seidel!(x, A::AbstractMatrix, b; maxiter::Int=10) x end -mutable struct DenseGaussSeidelIterable{matT,vecT} +mutable struct DenseGaussSeidelIterable{matT,solT,rhsT} A::matT - x::vecT - b::vecT + x::solT + b::rhsT maxiter::Int end @@ -149,11 +149,11 @@ function sor!(x, A::AbstractMatrix, b, ω::Real; maxiter::Int=10) x end -mutable struct DenseSORIterable{matT,vecT,numT} +mutable struct DenseSORIterable{matT,solT,vecT,rhsT,numT} A::matT - x::vecT + x::solT tmp::vecT - b::vecT + b::rhsT ω::numT maxiter::Int end @@ -207,11 +207,11 @@ function ssor!(x, A::AbstractMatrix, b, ω::Real; maxiter::Int=10) x end -mutable struct DenseSSORIterable{matT,vecT,numT} +mutable struct DenseSSORIterable{matT,solT,vecT,rhsT,numT} A::matT - x::vecT + x::solT tmp::vecT - b::vecT + b::rhsT ω::numT maxiter::Int end diff --git a/src/stationary_sparse.jl b/src/stationary_sparse.jl index 19c1bc25..ca7e2cef 100644 --- a/src/stationary_sparse.jl +++ b/src/stationary_sparse.jl @@ -208,12 +208,12 @@ end ## Jacobi ## -mutable struct JacobiIterable{T <: Number, vecT <: AbstractVector} +mutable struct JacobiIterable{T <: Number, solT, vecT, rhsT} O::OffDiagonal - x::vecT + x::solT next::vecT - b::vecT + b::rhsT maxiter::Int end @@ -229,8 +229,9 @@ function next(j::JacobiIterable{T}, iteration::Int) where {T} nothing, iteration + 1 end -jacobi_iterable(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector; maxiter::Int = 10) = - JacobiIterable{eltype(x), typeof(x)}(OffDiagonal(A, DiagonalIndices(A)), x, similar(x), b, maxiter) +jacobi_iterable(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector; maxiter::Int = 10) = + JacobiIterable{eltype(x), typeof(x), typeof(x), typeof(b)}( + OffDiagonal(A, DiagonalIndices(A)), x, similar(x), b, maxiter) """ @@ -253,12 +254,12 @@ end ## Gauss-Seidel ## -mutable struct GaussSeidelIterable{vecT <: AbstractVector} +mutable struct GaussSeidelIterable{solT, rhsT} U::StrictlyUpperTriangular L::FastLowerTriangular - x::vecT - b::vecT + x::solT + b::rhsT maxiter::Int end @@ -299,14 +300,14 @@ end ## SOR ## -mutable struct SORIterable{T, vecT <: AbstractVector, numT <: Real} +mutable struct SORIterable{T, solT, vecT, rhsT, numT <: Real} U::StrictlyUpperTriangular L::FastLowerTriangular ω::numT - x::vecT + x::solT next::vecT - b::vecT + b::rhsT maxiter::Int end @@ -329,7 +330,10 @@ end function sor_iterable(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector, ω::Real; maxiter::Int = 10) D = DiagonalIndices(A) T = eltype(x) - SORIterable{T,typeof(x),eltype(ω)}(StrictlyUpperTriangular(A, D), FastLowerTriangular(A, D), ω, x, similar(x), b, maxiter) + SORIterable{T,typeof(x),typeof(x),typeof(b),eltype(ω)}( + StrictlyUpperTriangular(A, D), FastLowerTriangular(A, D), ω, + x, similar(x), b, maxiter + ) end """ @@ -352,25 +356,27 @@ end ## SSOR ## -mutable struct SSORIterable{T, vecT, numT <: Real} +mutable struct SSORIterable{T, solT, vecT, rhsT, numT <: Real} sL::StrictlyLowerTriangular sU::StrictlyUpperTriangular L::FastLowerTriangular U::FastUpperTriangular ω::numT - x::vecT + x::solT tmp::vecT - b::vecT + b::rhsT maxiter::Int end function ssor_iterable(x::AbstractVector{T}, A::SparseMatrixCSC, b::AbstractVector, ω::Real; maxiter::Int = 10) where {T} D = DiagonalIndices(A) - sL = StrictlyLowerTriangular(A, D) - sU = StrictlyUpperTriangular(A, D) - L = FastLowerTriangular(A, D) - U = FastUpperTriangular(A, D) - SSORIterable{T,typeof(x),typeof(ω)}(sL, sU, L, U, ω, x, similar(x), b, maxiter) + SSORIterable{T,typeof(x),typeof(x),typeof(b),typeof(ω)}( + StrictlyLowerTriangular(A, D), + StrictlyUpperTriangular(A, D), + FastLowerTriangular(A, D), + FastUpperTriangular(A, D), + ω, x, similar(x), b, maxiter + ) end start(s::SSORIterable) = 1