From 2b1c402677556019709737da448d58c6e7ce1d99 Mon Sep 17 00:00:00 2001 From: Moritz Carmesin Date: Thu, 21 Nov 2024 10:47:19 +0100 Subject: [PATCH 1/9] Add tests for LinearExponential with GPU --- test/gpu/linear_exp.jl | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 test/gpu/linear_exp.jl diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl new file mode 100644 index 0000000000..634c49678d --- /dev/null +++ b/test/gpu/linear_exp.jl @@ -0,0 +1,35 @@ +using LinearAlgebra +using SparseArrays +using CUDA +using CUDA.CUSPARSE +using OrdinaryDiffEq + +# Linear exponential solvers +A = MatrixOperator([2.0 -1.0; -1.0 2.0]) +u0 = ones(2) + +A_gpu = MatrixOperator(cu([2.0 -1.0; -1.0 2.0])) +u0_gpu = cu(ones(2)) +prob_gpu = ODEProblem(A_gpu, u0_gpu, (0.0, 1.0)) + +sol_analytic = exp(1.0 * Matrix(A)) * u0 + +sol1_gpu = solve(prob_gpu, LinearExponential(krylov = :off))(1.0) |> Vector +sol2_gpu = solve(prob_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector +sol3_gpu = solve(prob_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector + +@test isapprox(sol1_gpu, sol_analytic, rtol = 1e-6) +@test isapprox(sol2_gpu, sol_analytic, rtol = 1e-6) +@test isapprox(sol3_gpu, sol_analytic, rtol = 1e-6) + +A2_gpu = MatrixOperator(cu(sparse([2.0 -1.0; -1.0 2.0]))) +prob2_gpu = ODEProblem(A2_gpu, u0_gpu, (0.0, 1.0)) + +sol2_1_gpu = solve(prob2_gpu, LinearExponential(krylov = :off))(1.0) |> Vector +sol2_2_gpu = solve(prob2_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector +sol2_3_gpu = solve(prob2_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector + +@test isapprox(sol2_1_gpu, sol_analytic, rtol = 1e-6) +@test isapprox(sol2_2_gpu, sol_analytic, rtol = 1e-6) +@test isapprox(sol2_3_gpu, sol_analytic, rtol = 1e-6) +@test isapprox(sol2_4_gpu, sol_analytic, rtol = 1e-4) \ No newline at end of file From 1e6685fa2043818537e4f5c41915bd76cef7bcfe Mon Sep 17 00:00:00 2001 From: Moritz Carmesin Date: Thu, 21 Nov 2024 10:52:28 +0100 Subject: [PATCH 2/9] Enable GPU tests for LinearExponential --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 3928ddad86..55b58486ab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -172,6 +172,7 @@ end end @time @safetestset "Autoswitch GPU" include("gpu/autoswitch.jl") @time @safetestset "Linear LSRK GPU" include("gpu/linear_lsrk.jl") + @time @safetestset "Linear Exponential GPU" include("gpu/linear_exp.jl") @time @safetestset "Reaction-Diffusion Stiff Solver GPU" include("gpu/reaction_diffusion_stiff.jl") @time @safetestset "Scalar indexing bug bypass" include("gpu/hermite_test.jl") end From a542c65642ad63c30cd673289aa493d2d7476e1d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 02:33:12 -0800 Subject: [PATCH 3/9] Update linear_caches.jl --- lib/OrdinaryDiffEqLinear/src/linear_caches.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl index 556444c0ca..a89f3288bb 100644 --- a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl +++ b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl @@ -573,7 +573,7 @@ function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int) u = zero(u_prototype) # stores the current state W = Matrix{T}(undef, n, p + 1) # stores the w vectors P = Matrix{T}(undef, n, p + 2) # stores output from phiv! - Ks = KrylovSubspace{T}(n, maxiter) # stores output from arnoldi! + Ks = KrylovSubspace{T,T,typeof(similar(u_prototype,size(u_prototype,1),2))}(n, maxiter) # stores output from arnoldi! phiv_cache = PhivCache(u_prototype, maxiter, p + 1) # cache used by phiv! (need +1 for error estimation) return u, W, P, Ks, phiv_cache end From cf2fec9671eebc6f0033f412076f455a1cc012dd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 02:46:29 -0800 Subject: [PATCH 4/9] Update lib/OrdinaryDiffEqLinear/src/linear_caches.jl --- lib/OrdinaryDiffEqLinear/src/linear_caches.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl index a89f3288bb..14aaeb0692 100644 --- a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl +++ b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl @@ -571,8 +571,8 @@ function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int) n = length(u_prototype) T = eltype(u_prototype) u = zero(u_prototype) # stores the current state - W = Matrix{T}(undef, n, p + 1) # stores the w vectors - P = Matrix{T}(undef, n, p + 2) # stores output from phiv! + W = similar(u_prototype, n, p+1) # stores the w vectors + P = similar(u_prototype, n, p+2) # stores output from phiv! Ks = KrylovSubspace{T,T,typeof(similar(u_prototype,size(u_prototype,1),2))}(n, maxiter) # stores output from arnoldi! phiv_cache = PhivCache(u_prototype, maxiter, p + 1) # cache used by phiv! (need +1 for error estimation) return u, W, P, Ks, phiv_cache From 236e987b958fc47f64d6e7d3a0fd3cd2b329a0e9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 14:14:57 -0800 Subject: [PATCH 5/9] Update test/gpu/linear_exp.jl --- test/gpu/linear_exp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl index 634c49678d..6c75557877 100644 --- a/test/gpu/linear_exp.jl +++ b/test/gpu/linear_exp.jl @@ -14,7 +14,7 @@ prob_gpu = ODEProblem(A_gpu, u0_gpu, (0.0, 1.0)) sol_analytic = exp(1.0 * Matrix(A)) * u0 -sol1_gpu = solve(prob_gpu, LinearExponential(krylov = :off))(1.0) |> Vector +@test_broken sol1_gpu = solve(prob_gpu, LinearExponential(krylov = :off))(1.0) |> Vector sol2_gpu = solve(prob_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector sol3_gpu = solve(prob_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector From 9243aabcbab9494a923482c1ae7ae30ed26b88cd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 14:15:02 -0800 Subject: [PATCH 6/9] Update test/gpu/linear_exp.jl --- test/gpu/linear_exp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl index 6c75557877..00b4f1b15f 100644 --- a/test/gpu/linear_exp.jl +++ b/test/gpu/linear_exp.jl @@ -18,7 +18,7 @@ sol_analytic = exp(1.0 * Matrix(A)) * u0 sol2_gpu = solve(prob_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector sol3_gpu = solve(prob_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector -@test isapprox(sol1_gpu, sol_analytic, rtol = 1e-6) +@test_broken isapprox(sol1_gpu, sol_analytic, rtol = 1e-6) @test isapprox(sol2_gpu, sol_analytic, rtol = 1e-6) @test isapprox(sol3_gpu, sol_analytic, rtol = 1e-6) From 1fe937d176ef187ffe3ba805de87e456a5a55a87 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 14:15:09 -0800 Subject: [PATCH 7/9] Update test/gpu/linear_exp.jl --- test/gpu/linear_exp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl index 00b4f1b15f..0dac2d3a7c 100644 --- a/test/gpu/linear_exp.jl +++ b/test/gpu/linear_exp.jl @@ -25,7 +25,7 @@ sol3_gpu = solve(prob_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector A2_gpu = MatrixOperator(cu(sparse([2.0 -1.0; -1.0 2.0]))) prob2_gpu = ODEProblem(A2_gpu, u0_gpu, (0.0, 1.0)) -sol2_1_gpu = solve(prob2_gpu, LinearExponential(krylov = :off))(1.0) |> Vector +@test_broken sol2_1_gpu = solve(prob2_gpu, LinearExponential(krylov = :off))(1.0) |> Vector sol2_2_gpu = solve(prob2_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector sol2_3_gpu = solve(prob2_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector From f15a48a4ce47cd922353d2c9e56137d7e28eab8c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 14:15:15 -0800 Subject: [PATCH 8/9] Update test/gpu/linear_exp.jl --- test/gpu/linear_exp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl index 0dac2d3a7c..94f2fb14d5 100644 --- a/test/gpu/linear_exp.jl +++ b/test/gpu/linear_exp.jl @@ -29,7 +29,7 @@ prob2_gpu = ODEProblem(A2_gpu, u0_gpu, (0.0, 1.0)) sol2_2_gpu = solve(prob2_gpu, LinearExponential(krylov = :simple))(1.0) |> Vector sol2_3_gpu = solve(prob2_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vector -@test isapprox(sol2_1_gpu, sol_analytic, rtol = 1e-6) +@test_broken isapprox(sol2_1_gpu, sol_analytic, rtol = 1e-6) @test isapprox(sol2_2_gpu, sol_analytic, rtol = 1e-6) @test isapprox(sol2_3_gpu, sol_analytic, rtol = 1e-6) @test isapprox(sol2_4_gpu, sol_analytic, rtol = 1e-4) \ No newline at end of file From ee8ea06cfa16380f339b952cf7f5737c9dd41b4b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Nov 2024 14:35:13 -0800 Subject: [PATCH 9/9] Update test/gpu/linear_exp.jl --- test/gpu/linear_exp.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/gpu/linear_exp.jl b/test/gpu/linear_exp.jl index 94f2fb14d5..dbbde584cc 100644 --- a/test/gpu/linear_exp.jl +++ b/test/gpu/linear_exp.jl @@ -31,5 +31,4 @@ sol2_3_gpu = solve(prob2_gpu, LinearExponential(krylov = :adaptive))(1.0) |> Vec @test_broken isapprox(sol2_1_gpu, sol_analytic, rtol = 1e-6) @test isapprox(sol2_2_gpu, sol_analytic, rtol = 1e-6) -@test isapprox(sol2_3_gpu, sol_analytic, rtol = 1e-6) -@test isapprox(sol2_4_gpu, sol_analytic, rtol = 1e-4) \ No newline at end of file +@test isapprox(sol2_3_gpu, sol_analytic, rtol = 1e-6) \ No newline at end of file