From 974ff9e96316476484e9ef1ca1b419502a83f9e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 24 Jan 2023 19:47:10 +0100 Subject: [PATCH 1/2] steps towards iterative solvers + tests for general numbers --- src/iterative_wrappers.jl | 8 ++-- test/basictests.jl | 94 ++++++++++++++++++++++++--------------- 2 files changed, 64 insertions(+), 38 deletions(-) diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index 7386e37b7..f35d873b1 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -153,11 +153,13 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...) M = (M === Identity()) ? I : InvPreconditioner(M) N = (N === Identity()) ? I : InvPreconditioner(N) - atol = float(cache.abstol) - rtol = float(cache.reltol) + Ta = eltype(cache.A) + + atol = Ta(float(cache.abstol)) + rtol = Ta(float(cache.reltol)) itmax = cache.maxiters verbose = cache.verbose ? 1 : 0 - + args = (cache.cacheval, cache.A, cache.b) kwargs = (atol = atol, rtol = rtol, itmax = itmax, verbose = verbose, history = true, alg.kwargs...) diff --git a/test/basictests.jl b/test/basictests.jl index dfe6cb707..aed86cac9 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -3,6 +3,8 @@ using Test import Random const Dual64 = ForwardDiff.Dual{Nothing, Float64, 1} +Base.:^(x::MultiFloat{T, N}, y::Int) where {T,N} = MultiFloat{T, N}(BigFloat(x)^y) +Base.:^(x::MultiFloat{T, N}, y::Float64) where {T,N} = MultiFloat{T, N}(BigFloat(x)^y) n = 8 A = Matrix(I, n, n) @@ -19,18 +21,21 @@ prob2 = LinearProblem(A2, b2; u0 = x2) cache_kwargs = (; verbose = true, abstol = 1e-8, reltol = 1e-8, maxiter = 30) -function test_interface(alg, prob1, prob2) - A1 = prob1.A - b1 = prob1.b - x1 = prob1.u0 - A2 = prob2.A - b2 = prob2.b - x2 = prob2.u0 - - y = solve(prob1, alg; cache_kwargs...) +function test_interface(alg, prob1, prob2; T=Float64) + A1 = prob1.A .|> T + b1 = prob1.b .|> T + x1 = prob1.u0 .|> T + A2 = prob2.A .|> T + b2 = prob2.b .|> T + x2 = prob2.u0 .|> T + + myprob1 = LinearProblem(A1, b1; u0 = x1) + myprob2 = LinearProblem(A2, b2; u0 = x2) + + y = solve(myprob1, alg; cache_kwargs...) @test A1 * y ≈ b1 - cache = SciMLBase.init(prob1, alg; cache_kwargs...) # initialize cache + cache = SciMLBase.init(myprob1, alg; cache_kwargs...) # initialize cache y = solve(cache) @test A1 * y ≈ b1 @@ -140,44 +145,44 @@ end end @testset "Sparspak Factorization (Float64x1)" begin - A1 = sparse(A / 1) .|> Float64x1 - b1 = rand(n) .|> Float64x1 - x1 = zero(b) .|> Float64x1 - A2 = sparse(A / 2) .|> Float64x1 - b2 = rand(n) .|> Float64x1 - x2 = zero(b) .|> Float64x1 + A1 = sparse(A / 1) + b1 = rand(n) + x1 = zero(b) + A2 = sparse(A / 2) + b2 = rand(n) + x2 = zero(b) prob1 = LinearProblem(A1, b1; u0 = x1) prob2 = LinearProblem(A2, b2; u0 = x2) - test_interface(SparspakFactorization(), prob1, prob2) + test_interface(SparspakFactorization(), prob1, prob2; T=Float64x1) end @testset "Sparspak Factorization (Float64x2)" begin - A1 = sparse(A / 1) .|> Float64x2 - b1 = rand(n) .|> Float64x2 - x1 = zero(b) .|> Float64x2 - A2 = sparse(A / 2) .|> Float64x2 - b2 = rand(n) .|> Float64x2 - x2 = zero(b) .|> Float64x2 + A1 = sparse(A / 1) + b1 = rand(n) + x1 = zero(b) + A2 = sparse(A / 2) + b2 = rand(n) + x2 = zero(b) prob1 = LinearProblem(A1, b1; u0 = x1) prob2 = LinearProblem(A2, b2; u0 = x2) - test_interface(SparspakFactorization(), prob1, prob2) + test_interface(SparspakFactorization(), prob1, prob2; T=Float64x2) end @testset "Sparspak Factorization (Dual64)" begin - A1 = sparse(A / 1) .|> Dual64 - b1 = rand(n) .|> Dual64 - x1 = zero(b) .|> Dual64 - A2 = sparse(A / 2) .|> Dual64 - b2 = rand(n) .|> Dual64 - x2 = zero(b) .|> Dual64 + A1 = sparse(A / 1) + b1 = rand(n) + x1 = zero(b) + A2 = sparse(A / 2) + b2 = rand(n) + x2 = zero(b) prob1 = LinearProblem(A1, b1; u0 = x1) prob2 = LinearProblem(A2, b2; u0 = x2) - test_interface(SparspakFactorization(), prob1, prob2) + test_interface(SparspakFactorization(), prob1, prob2; T=Dual64) end - + @testset "FastLAPACK Factorizations" begin A1 = A / 1 b1 = rand(n) @@ -225,7 +230,14 @@ end ("GMRES", KrylovJL_GMRES(kwargs...)), # ("BICGSTAB",KrylovJL_BICGSTAB(kwargs...)), ("MINRES", KrylovJL_MINRES(kwargs...))) - @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) end + @testset "$(alg[1])" begin + test_interface(alg[2], prob1, prob2) + test_interface(alg[2], prob1, prob2; T=Float64x1) + test_interface(alg[2], prob1, prob2; T=Float64x2) + # test_interface(alg[2], prob1, prob2; T=Dual64) + # https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/646 + # ForwardDiff.Dual is a Real, not an AbstractFloat + end end end @@ -237,7 +249,14 @@ end # ("BICGSTAB",IterativeSolversJL_BICGSTAB(kwargs...)), # ("MINRES",IterativeSolversJL_MINRES(kwargs...)), ) - @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) end + @testset "$(alg[1])" begin + test_interface(alg[2], prob1, prob2) + test_interface(alg[2], prob1, prob2; T=Float64x1) + test_interface(alg[2], prob1, prob2; T=Float64x2) + # test_interface(alg[2], prob1, prob2; T=Dual64) + # https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/givens.jl#L77 + # ForwardDiff.Dual is a Real, not an AbstractFloat + end end end @@ -246,7 +265,12 @@ end for alg in (("Default", KrylovKitJL(kwargs...)), ("CG", KrylovKitJL_CG(kwargs...)), ("GMRES", KrylovKitJL_GMRES(kwargs...))) - @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) end + @testset "$(alg[1])" begin + test_interface(alg[2], prob1, prob2) + test_interface(alg[2], prob1, prob2; T=Float64x1) + test_interface(alg[2], prob1, prob2; T=Float64x2) + test_interface(alg[2], prob1, prob2; T=Dual64) + end @test alg[2] isa KrylovKitJL end end From fc4ae656cdbff46ca7e40efddcc5bd599001aa55 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 24 Jan 2023 21:04:20 +0100 Subject: [PATCH 2/2] format; let failing tests for duals ring --- src/iterative_wrappers.jl | 4 ++-- test/basictests.jl | 36 ++++++++++++++++++------------------ 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index f35d873b1..86f01ecee 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -154,12 +154,12 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...) N = (N === Identity()) ? I : InvPreconditioner(N) Ta = eltype(cache.A) - + atol = Ta(float(cache.abstol)) rtol = Ta(float(cache.reltol)) itmax = cache.maxiters verbose = cache.verbose ? 1 : 0 - + args = (cache.cacheval, cache.A, cache.b) kwargs = (atol = atol, rtol = rtol, itmax = itmax, verbose = verbose, history = true, alg.kwargs...) diff --git a/test/basictests.jl b/test/basictests.jl index aed86cac9..cf3998f8a 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -3,8 +3,8 @@ using Test import Random const Dual64 = ForwardDiff.Dual{Nothing, Float64, 1} -Base.:^(x::MultiFloat{T, N}, y::Int) where {T,N} = MultiFloat{T, N}(BigFloat(x)^y) -Base.:^(x::MultiFloat{T, N}, y::Float64) where {T,N} = MultiFloat{T, N}(BigFloat(x)^y) +Base.:^(x::MultiFloat{T, N}, y::Int) where {T, N} = MultiFloat{T, N}(BigFloat(x)^y) +Base.:^(x::MultiFloat{T, N}, y::Float64) where {T, N} = MultiFloat{T, N}(BigFloat(x)^y) n = 8 A = Matrix(I, n, n) @@ -21,7 +21,7 @@ prob2 = LinearProblem(A2, b2; u0 = x2) cache_kwargs = (; verbose = true, abstol = 1e-8, reltol = 1e-8, maxiter = 30) -function test_interface(alg, prob1, prob2; T=Float64) +function test_interface(alg, prob1, prob2; T = Float64) A1 = prob1.A .|> T b1 = prob1.b .|> T x1 = prob1.u0 .|> T @@ -31,7 +31,7 @@ function test_interface(alg, prob1, prob2; T=Float64) myprob1 = LinearProblem(A1, b1; u0 = x1) myprob2 = LinearProblem(A2, b2; u0 = x2) - + y = solve(myprob1, alg; cache_kwargs...) @test A1 * y ≈ b1 @@ -154,7 +154,7 @@ end prob1 = LinearProblem(A1, b1; u0 = x1) prob2 = LinearProblem(A2, b2; u0 = x2) - test_interface(SparspakFactorization(), prob1, prob2; T=Float64x1) + test_interface(SparspakFactorization(), prob1, prob2; T = Float64x1) end @testset "Sparspak Factorization (Float64x2)" begin @@ -167,7 +167,7 @@ end prob1 = LinearProblem(A1, b1; u0 = x1) prob2 = LinearProblem(A2, b2; u0 = x2) - test_interface(SparspakFactorization(), prob1, prob2; T=Float64x2) + test_interface(SparspakFactorization(), prob1, prob2; T = Float64x2) end @testset "Sparspak Factorization (Dual64)" begin @@ -180,9 +180,9 @@ end prob1 = LinearProblem(A1, b1; u0 = x1) prob2 = LinearProblem(A2, b2; u0 = x2) - test_interface(SparspakFactorization(), prob1, prob2; T=Dual64) + test_interface(SparspakFactorization(), prob1, prob2; T = Dual64) end - + @testset "FastLAPACK Factorizations" begin A1 = A / 1 b1 = rand(n) @@ -232,9 +232,9 @@ end ("MINRES", KrylovJL_MINRES(kwargs...))) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) - test_interface(alg[2], prob1, prob2; T=Float64x1) - test_interface(alg[2], prob1, prob2; T=Float64x2) - # test_interface(alg[2], prob1, prob2; T=Dual64) + test_interface(alg[2], prob1, prob2; T = Float64x1) + test_interface(alg[2], prob1, prob2; T = Float64x2) + test_interface(alg[2], prob1, prob2; T = Dual64) # https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/646 # ForwardDiff.Dual is a Real, not an AbstractFloat end @@ -251,10 +251,10 @@ end ) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) - test_interface(alg[2], prob1, prob2; T=Float64x1) - test_interface(alg[2], prob1, prob2; T=Float64x2) - # test_interface(alg[2], prob1, prob2; T=Dual64) - # https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/givens.jl#L77 + test_interface(alg[2], prob1, prob2; T = Float64x1) + test_interface(alg[2], prob1, prob2; T = Float64x2) + test_interface(alg[2], prob1, prob2; T = Dual64) + # https://github.com/JuliaLang/julia/issues/41753 # ForwardDiff.Dual is a Real, not an AbstractFloat end end @@ -267,9 +267,9 @@ end ("GMRES", KrylovKitJL_GMRES(kwargs...))) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) - test_interface(alg[2], prob1, prob2; T=Float64x1) - test_interface(alg[2], prob1, prob2; T=Float64x2) - test_interface(alg[2], prob1, prob2; T=Dual64) + test_interface(alg[2], prob1, prob2; T = Float64x1) + test_interface(alg[2], prob1, prob2; T = Float64x2) + test_interface(alg[2], prob1, prob2; T = Dual64) end @test alg[2] isa KrylovKitJL end