### Install DiffEqGPU.jl

Source: https://docs.sciml.ai/DiffEqGPU/stable

Use the Julia package manager to add the DiffEqGPU package to your environment.

```julia
using Pkg
Pkg.add("DiffEqGPU")
```

--------------------------------

### Batched Reduction Example with EnsembleGPUArray

Source: https://docs.sciml.ai/DiffEqGPU/stable/examples/reductions

This example demonstrates running 20 ODEs at a time, calculating the final value, and reducing the results by summing them. It uses EnsembleGPUArray to process 100 trajectories with a batch size of 20.

```julia
using OrdinaryDiffEq, DiffEqGPU, CUDA

seed = 100
using Random;
Random.seed!(seed);
ra = rand(100)

function f!(du, u, p, t)
    du[1] = 1.01 * u[1]
end

prob = ODEProblem(f!, [0.5], (0.0, 1.0))

function output_func(sol, i)
    last(sol), false
end

function prob_func(prob, i, repeat)
    remake(prob, u0 = ra[i] * prob.u0)
end

function reduction(u, batch, I)
    u .+ sum(batch), false
end

prob2 = EnsembleProblem(prob, prob_func = prob_func, output_func = output_func,
    reduction = reduction, u_init = Vector{eltype(prob.u0)}([0.0]))
sim4 = solve(prob2, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()), trajectories = 100,
    batch_size = 20)
```

```text
EnsembleSolution Solution of length 1 with uType:
Float64
```

--------------------------------

### Setup processes with different CUDA devices

Source: https://docs.sciml.ai/DiffEqGPU/stable/tutorials/multigpu

Initializes distributed processes and assigns each worker to a specific CUDA device.

```julia
# Setup processes with different CUDA devices
using Distributed
numgpus = 1
addprocs(numgpus)
import CUDA

let gpuworkers = asyncmap(collect(zip(workers(), CUDA.devices()))) do (p, d)
        remotecall_wait(CUDA.device!, p, d)
        p
    end
end
```

--------------------------------

### Solve an ODE on the GPU

Source: https://docs.sciml.ai/DiffEqGPU/stable/getting_started

This example shows how to define an ODEProblem and solve it using GPU-accelerated arrays. Note that Float32 is recommended for better performance on GPUs.

```julia
using OrdinaryDiffEq, CUDA, LinearAlgebra
u0 = cu(rand(1000))
A = cu(randn(1000, 1000))
f(du, u, p, t) = mul!(du, A, u)
prob = ODEProblem(f, u0, (0.0f0, 1.0f0)) # Float32 is better on GPUs!
sol = solve(prob, Tsit5())
```

--------------------------------

### Multi-GPU system script

Source: https://docs.sciml.ai/DiffEqGPU/stable/tutorials/multigpu

A complete example for a 2-GPU system where devices are explicitly assigned to worker processes.

```julia
using DiffEqGPU, CUDA, OrdinaryDiffEq, Test
CUDA.device!(0)

using Distributed
addprocs(2)
@everywhere using DiffEqGPU, CUDA, OrdinaryDiffEq, Test, Random

@everywhere begin
    function lorenz_distributed(du, u, p, t)
        du[1] = p[1] * (u[2] - u[1])
        du[2] = u[1] * (p[2] - u[3]) - u[2]
        du[3] = u[1] * u[2] - p[3] * u[3]
    end
    CUDA.allowscalar(false)
    u0 = Float32[1.0; 0.0; 0.0]
    tspan = (0.0f0, 100.0f0)
    p = [10.0f0, 28.0f0, 8 / 3.0f0]
    Random.seed!(1)
    pre_p_distributed = [rand(Float32, 3) for i in 1:100_000]
    function prob_func_distributed(prob, i, repeat)
        remake(prob, p = pre_p_distributed[i] .* p)
    end
end

@sync begin
    @spawnat 2 begin
        CUDA.allowscalar(false)
        CUDA.device!(0)
    end
    @spawnat 3 begin
        CUDA.allowscalar(false)
        CUDA.device!(1)
    end
end

CUDA.allowscalar(false)
prob = ODEProblem(lorenz_distributed, u0, tspan, p)
monteprob = EnsembleProblem(prob, prob_func = prob_func_distributed)

@time sol = solve(monteprob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()),
    trajectories = 100_000,
    batch_size = 50_000, saveat = 1.0f0)
```

--------------------------------

### GPU-accelerated reaction-diffusion ODE setup with CuArray

Source: https://docs.sciml.ai/DiffEqGPU/stable/examples/reaction_diffusion

Sets up the reaction-diffusion ODE problem on the GPU by converting matrices and initial condition to CuArray. The function gf uses views and in-place matrix multiplication (LA.mul!) with broadcasted updates. CUDA.allowscalar(false) disables slow scalar fallbacks.

```julia
import CUDA
const gMx = CUDA.CuArray(Float32.(Mx))
const gMy = CUDA.CuArray(Float32.(My))
const gα₁ = CUDA.CuArray(Float32.(α₁))
gu0 = CUDA.CuArray(Float32.(u0))

const gMyA = CUDA.CuArray(zeros(Float32, N, N))
const AgMx = CUDA.CuArray(zeros(Float32, N, N))
const gDA = CUDA.CuArray(zeros(Float32, N, N))
function gf(du, u, p, t)
    A = @view u[:, :, 1]
    B = @view u[:, :, 2]
    C = @view u[:, :, 3]
    dA = @view du[:, :, 1]
    dB = @view du[:, :, 2]
    dC = @view du[:, :, 3]
    LA.mul!(gMyA, gMy, A)
    LA.mul!(AgMx, A, gMx)
    @. gDA = D * (gMyA + AgMx)
    @. dA = gDA + gα₁ - β₁ * A - r₁ * A * B + r₂ * C
    @. dB = α₂ - β₂ * B - r₁ * A * B + r₂ * C
    @. dC = α₃ - β₃ * C + r₁ * A * B - r₂ * C
end

prob2 = DE.ODEProblem(gf, gu0, (0.0, 100.0))
CUDA.allowscalar(false) # makes sure none of the slow fallbacks are used
@time sol = DE.solve(prob2, DE.ROCK2(), progress = true, dt = 0.003, save_everystep = false,
    save_start = false);
```

--------------------------------

### Forward-mode AD solution output

Source: https://docs.sciml.ai/DiffEqGPU/stable/examples/ad

The resulting EnsembleSolution object after running the forward-mode AD example.

```text
EnsembleSolution Solution of length 10000 with uType:
SciMLBase.ODESolution{ForwardDiff.Dual{Nothing, Float64, 3}, 2, uType, Nothing, Nothing, Vector{Float32}, rateType, SciMLBase.ODEProblem{Vector{ForwardDiff.Dual{Nothing, Float64, 3}}, Tuple{Float32, Float32}, true, Vector{Float32}, SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.lorenz), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, IType, SciMLBase.DEStats, Nothing} where {uType, rateType, IType}
```

--------------------------------

### EnsembleGPUArray Lorenz example

Source: https://docs.sciml.ai/DiffEqGPU/stable/manual/ensemblegpuarray

Example of using EnsembleGPUArray with a Lorenz system on a CUDA GPU. The prob_func randomizes parameters for each trajectory, and the solve call runs 10,000 trajectories with Tsit5 and saveat=1.0f0.

```julia
using DiffEqGPU, CUDA, OrdinaryDiffEq
function lorenz(du, u, p, t)
    du[1] = p[1] * (u[2] - u[1])
    du[2] = u[1] * (p[2] - u[3]) - u[2]
    du[3] = u[1] * u[2] - p[3] * u[3]
end

u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = [10.0f0,28.0f0,8/3f0]
prob = ODEProblem(lorenz,u0,tspan,p)
prob_func = (prob,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(CUDADevice()),trajectories=10_000,saveat=1.0f0)
```

--------------------------------

### Solve Reaction-Diffusion PDE

Source: https://docs.sciml.ai/DiffEqGPU/stable/examples/reaction_diffusion

Setup and solve the ODEProblem using the defined reaction-diffusion function.

```julia
import DifferentialEquations as DE

prob = DE.ODEProblem(f, u0, (0.0, 100.0))
@time sol = DE.solve(prob, DE.ROCK2());
```

```julia
retcode: Success
Interpolation: 3rd order Hermite
t: 125-element Vector{Float64}:
   0.0
   2.796123741859731e-5
   6.43785906507212e-5
   0.00013556965396987797
   0.00023376199344619583
   0.00037760221289476507
   0.0005650578854793943
   0.0008071655769861446
   0.001105353327393115
   0.0014666858344513482
   ⋮
   7.9506232854699155
   9.220778312424493
  11.055370462175839
  13.943634271456988
  18.834203937890415
  27.13897752197475
  43.69184525085535
  74.9721648444126
 100.0
u: 125-element Vector{Array{Float64, 3}}:
 [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
 [3.906344697357476e-10 3.909112473032509e-10 … 2.796119589268521e-5 2.7883165317332937e-5; 3.909112473032509e-10 3.911880248707541e-10 … 2.8039337179064475e-5 2.796119589268521e-5; … ; 3.909112473032509e-10 3.911880248707541e-10 … 2.8039337179064475e-5 2.796119589268521e-5; 3.906344697357476e-10 3.909112473032509e-10 … 2.796119589268521e-5 2.7883165317332937e-5;;; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.7961237409320334e-5 2.79612374093249e-5; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.796123740931577e-5 2.7961237409320334e-5; … ; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.796123740931577e-5 2.7961237409320334e-5; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.7961237409320334e-5 2.79612374093249e-5;;; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.796045559984408e-5 2.7960455599839515e-5; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.7960455599848643e-5 2.796045559984408e-5; … ; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.7960455599848643e-5 2.796045559984408e-5; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.796045559984408e-5 2.7960455599839515e-5]
 [2.064777884980437e-9 2.072178013484094e-9 … 6.437748417993337e-5 6.396710011488427e-5; 2.072178013484094e-9 2.0796095725767602e-9 … 6.479081264774206e-5 6.437748417993337e-5; … ; 2.072178013484094e-9 2.0796095725767602e-9 … 6.479081264774206e-5 6.437748417993337e-5; 2.064777884980437e-9 2.072178013484094e-9 … 6.437748417993337e-5 6.396710011488427e-5;;; 6.437859057628412e-5 6.437859057628411e-5 … 6.437859048984323e-5 6.437859049019145e-5; 6.437859057628411e-5 6.437859057628411e-5 … 6.437859048949346e-5 6.437859048984323e-5; … ; 6.437859057628411e-5 6.437859057628411e-5 … 6.437859048949346e-5 6.437859048984323e-5; 6.437859057628412e-5 6.437859057628411e-5 … 6.437859048984323e-5 6.437859049019145e-5;;; 6.4374446196661e-5 6.437444619666102e-5 … 6.437444628310191e-5 6.437444628275368e-5; 6.437444619666102e-5 6.437444619666102e-5 … 6.437444628345167e-5 6.437444628310191e-5; … ; 6.437444619666102e-5 6.437444619666102e-5 … 6.437444628345167e-5 6.437444628310191e-5; 6.4374446196661e-5 6.437444619666102e-5 … 6.437444628310191e-5 6.437444628275368e-5]
```

--------------------------------

### Solve Lorenz ODE ensemble on GPU with DiffEqGPU

Source: https://docs.sciml.ai/DiffEqGPU/stable/getting_started

Sets up a Lorenz ODE with static arrays and solves 10,000 parameter variations in parallel on a GPU using GPUTsit5 and EnsembleGPUKernel with CUDA.CUDABackend. Requires DiffEqGPU, OrdinaryDiffEq, StaticArrays, and CUDA packages.

```julia
using DiffEqGPU, OrdinaryDiffEq, StaticArrays, CUDA

function lorenz(u, p, t)
    σ = p[1]
    ρ = p[2]
    β = p[3]
    du1 = σ * (u[2] - u[1])
    du2 = u[1] * (ρ - u[3]) - u[2]
    du3 = u[1] * u[2] - β * u[3]
    return SVector{3}(du1, du2, du3)
end

u0 = @SVector [1.0f0; 0.0f0; 0.0f0]
tspan = (0.0f0, 10.0f0)
p = @SVector [10.0f0, 28.0f0, 8 / 3.0f0]
prob = ODEProblem{false}(lorenz, u0, tspan, p)
prob_func = (prob, i, repeat) -> remake(prob, p = (@SVector rand(Float32, 3)) .* p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false)

sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()),
    trajectories = 10_000)
```

```text
EnsembleSolution Solution of length 10000 with uType:
SciMLBase.ODESolution{Float32, 2, SubArray{StaticArraysCore.SVector{3, Float32}, 1, Matrix{StaticArraysCore.SVector{3, Float32}}, Tuple{UnitRange{Int64}, Int64}, true}, Nothing, Nothing, SubArray{Float32, 1, Matrix{Float32}, Tuple{UnitRange{Int64}, Int64}, true}, Nothing, DiffEqGPU.ImmutableODEProblem{StaticArraysCore.SVector{3, Float32}, Tuple{Float32, Float32}, false, StaticArraysCore.SVector{3, Float32}, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.lorenz), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, GPUTsit5, SciMLBase.LinearInterpolation{SubArray{Float32, 1, Matrix{Float32}, Tuple{UnitRange{Int64}, Int64}, true}, SubArray{StaticArraysCore.SVector{3, Float32}, 1, Matrix{StaticArraysCore.SVector{3, Float32}}, Tuple{UnitRange{Int64}, Int64}, true}}, Nothing, Nothing}
```

--------------------------------

### Lower-level GPU ensemble solve with DiffEqGPU.vectorized_solve and vectorized_asolve

Source: https://docs.sciml.ai/DiffEqGPU/stable/tutorials/lower_level_api

Use the lower-level API of EnsembleGPUKernel to reduce overhead by manually converting problems to GPU arrays and calling DiffEqGPU.vectorized_solve (fixed time-stepping) or DiffEqGPU.vectorized_asolve (adaptive time-stepping). The first call for each solver is for compilation; the second call measures the actual solve time. The output shows the returned solution structure with CuArray columns.

```julia
using DiffEqGPU, OrdinaryDiffEq, StaticArrays, CUDA, DiffEqBase

trajectories = 10_000

function lorenz(u, p, t)
    σ = p[1]
    ρ = p[2]
    β = p[3]
    du1 = σ * (u[2] - u[1])
    du2 = u[1] * (ρ - u[3]) - u[2]
    du3 = u[1] * u[2] - β * u[3]
    return SVector{3}(du1, du2, du3)
end

u0 = @SVector [1.0f0; 0.0f0; 0.0f0]
tspan = (0.0f0, 10.0f0)
p = @SVector [10.0f0, 28.0f0, 8 / 3.0f0]
prob = ODEProblem{false}(lorenz, u0, tspan, p)

## Building different problems for different parameters
probs = map(1:trajectories) do i
    DiffEqGPU.make_prob_compatible(remake(prob, p = (@SVector rand(Float32, 3)) .* p))
end

## Move the arrays to the GPU
probs = cu(probs)

## Finally use the lower API for faster solves! (Fixed time-stepping)

# Run once for compilation
@time CUDA.@sync ts, us = DiffEqGPU.vectorized_solve(probs, prob, GPUTsit5();
    save_everystep = false, dt = 0.1f0)

@time CUDA.@sync ts, us = DiffEqGPU.vectorized_solve(probs, prob, GPUTsit5();
    save_everystep = false, dt = 0.1f0)

## Adaptive time-stepping (Notice the boolean argument)
# Run once for compilation
@time CUDA.@sync ts, us = DiffEqGPU.vectorized_asolve(probs, prob, GPUTsit5();
    save_everystep = false, dt = 0.1f0)

@time CUDA.@sync ts, us = DiffEqGPU.vectorized_asolve(probs, prob, GPUTsit5();
    save_everystep = false, dt = 0.1f0)
```

```julia
(Float32[0.0 0.0 … 0.0 0.0; 10.0 10.0 … 10.0 10.0], StaticArraysCore.SVector{3, Float32}[[1.0, 0.0, 0.0] [1.0, 0.0, 0.0] … [1.0, 0.0, 0.0] [1.0, 0.0, 0.0]; [5.2474143f-6, 4.531192f-6, 8.4393784f-11] [2.752024, 2.7527635, 3.5590749] … [-5.1432833, -5.6739984, 14.661653] [-2.471736, -2.490478, 20.2666]])
```

--------------------------------

### Reverse-mode AD with Flux training loop

Source: https://docs.sciml.ai/DiffEqGPU/stable/examples/ad

Example of using EnsembleGPUArray within a Flux training loop. Note that CUDA.allowscalar(false) is used to ensure GPU compatibility.

```julia
using OrdinaryDiffEq, SciMLSensitivity, Flux, DiffEqGPU, CUDA, Test
CUDA.allowscalar(false)

function modelf(du, u, p, t)
    du[1] = 1.01 * u[1] * p[1] * p[2]
end

function model()
    prob = ODEProblem(modelf, u0, (0.0, 1.0), pa)

    function prob_func(prob, i, repeat)
        remake(prob, u0 = 0.5 .+ i / 100 .* prob.u0)
    end

    ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
    solve(ensemble_prob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()), saveat = 0.1,
        trajectories = 10)
end

# loss function
loss() = sum(abs2, 1.0 .- Array(model()))

data = Iterators.repeated((), 10)

cb = function () # callback function to observe training
    @show loss()
end

pa = [1.0, 2.0]
u0 = [3.0]
opt = ADAM(0.1)
println("Starting to train")

l1 = loss()

Flux.train!(loss, Flux.params([pa]), data, opt; cb = cb)
```

--------------------------------

### Project.toml Dependencies

Source: https://docs.sciml.ai/DiffEqGPU/stable

List of direct dependencies used to build the package documentation.

```toml
Status `/var/lib/buildkite-agent/builds/gpuci-17/julialang/diffeqgpu-dot-jl/docs/Project.toml`
  [79e6a3ab] Adapt v3.7.1
  [6e4b80f9] BenchmarkTools v1.3.2
  [052768ef] CUDA v5.1.0
  [2b5f629d] DiffEqBase v6.138.1
  [071ae1c0] DiffEqGPU v3.3.0 `/var/lib/buildkite-agent/builds/gpuci-17/julialang/diffeqgpu-dot-jl`
  [e30172f5] Documenter v1.1.2
  [587475ba] Flux v0.14.6
  [f6369f11] ForwardDiff v0.10.36
  [1dea7af3] OrdinaryDiffEq v6.59.0
  [91a5bcdd] Plots v1.39.0
  [1bc83da4] SafeTestsets v0.1.0
  [1ed8b502] SciMLSensitivity v7.47.0
  [90137ffa] StaticArrays v1.6.5
  [789caeaf] StochasticDiffEq v6.63.2
  [8ba89e20] Distributed
  [37e2e46d] LinearAlgebra
  [9a3f8284] Random
  [10745b16] Statistics v1.9.0
```

--------------------------------

### Solve Lorenz ensemble with EnsembleGPUKernel and GPUTsit5

Source: https://docs.sciml.ai/DiffEqGPU/stable/tutorials/gpu_ensemble_basic

Use EnsembleGPUKernel with GPUTsit5 to solve an ensemble of Lorenz ODEs on a GPU. The problem must be written in out-of-place form with StaticArrays and use special solvers. The output shows the resulting EnsembleSolution type with 10000 trajectories.

```julia
using DiffEqGPU, OrdinaryDiffEq, StaticArrays, CUDA

function lorenz2(u, p, t)
    σ = p[1]
    ρ = p[2]
    β = p[3]
    du1 = σ * (u[2] - u[1])
    du2 = u[1] * (ρ - u[3]) - u[2]
    du3 = u[1] * u[2] - β * u[3]
    return SVector{3}(du1, du2, du3)
end

u0 = @SVector [1.0f0; 0.0f0; 0.0f0]
tspan = (0.0f0, 10.0f0)
p = @SVector [10.0f0, 28.0f0, 8 / 3.0f0]
prob = ODEProblem{false}(lorenz2, u0, tspan, p)
prob_func = (prob, i, repeat) -> remake(prob, p = (@SVector rand(Float32, 3)) .* p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false)
sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()),
    trajectories = 10_000,
    saveat = 1.0f0)
```

```text
EnsembleSolution Solution of length 10000 with uType:
SciMLBase.ODESolution{Float32, 2, SubArray{StaticArraysCore.SVector{3, Float32}, 1, Matrix{StaticArraysCore.SVector{3, Float32}}, Tuple{UnitRange{Int64}, Int64}, true}, Nothing, SubArray{Float32, 1, Matrix{Float32}, Tuple{UnitRange{Int64}, Int64}, true}, Nothing, DiffEqGPU.ImmutableODEProblem{StaticArraysCore.SVector{3, Float32}, Tuple{Float32, Float32}, false, StaticArraysCore.SVector{3, Float32}, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.lorenz2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, GPUTsit5, SciMLBase.LinearInterpolation{SubArray{Float32, 1, Matrix{Float32}, Tuple{UnitRange{Int64}, Int64}, true}, SubArray{StaticArraysCore.SVector{3, Float32}, 1, Matrix{StaticArraysCore.SVector{3, Float32}}, Tuple{UnitRange{Int64}, Int64}, true}}, Nothing, Nothing}
```

--------------------------------

### Within-Method GPU Parallelism with CuArray and Tsit5

Source: https://docs.sciml.ai/DiffEqGPU/stable/getting_started

Solve a large ODE with within-method GPU parallelism using CUDA.jl's CuArray as u0. The example defines a 1000-dimensional ODE with a matrix multiplication and solves it with Tsit5. The output shows the solution's retcode, interpolation type, time vector, and the solution values.

```julia
using OrdinaryDiffEq, LinearAlgebra
u0 = rand(1000)
A = randn(1000, 1000)
f(du, u, p, t) = mul!(du, A, u)
prob = ODEProblem(f, u0, (0.0, 1.0))
sol = solve(prob, Tsit5())
```

```text
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 50-element Vector{Float64}:
 0.0
 0.0015478761855679305
 0.00648580863398825
 0.013842416823521709
 0.022391643841718323
 0.03325242939710622
 0.04561520602822331
 0.05901456341425282
 0.07299580669091357
 0.087920728422522
 ⋮
 0.8051958359073149
 0.8303464269733638
 0.855220989124135
 0.8807960944818966
 0.906252317657025
 0.9322591011728141
 0.9588072465756187
 0.9814855255365785
 1.0
u: 50-element Vector{Vector{Float64}}:
 [0.6077993951026277, 0.7492427325037554, 0.24030948236125937, 0.7876203548496596, 0.9790896488924092, 0.39466961581632154, 0.8624005981887551, 0.44259694848939535, 0.8924620923795857, 0.02291634363312045  …  0.9575521392611473, 0.8803300650478935, 0.5201892514317958, 0.9654442939425254, 0.36768574375921526, 0.6396059294669751, 0.901286174640774, 0.685994963053012, 0.27133784540119166, 0.3416276233225435]
 [0.6388217639528583, 0.7550178758102967, 0.23851886858256294, 0.798157641408362, 0.9991230598196218, 0.3979320485800568, 0.8199383091064789, 0.41758245795126836, 0.8625227443916994, 0.025110186525423902  …  0.9928820975413687, 0.9188273924595998, 0.5424605821203632, 0.9989534800028482, 0.34863444594913573, 0.6129502249468205, 0.9449488603002396, 0.6859185020749737, 0.2372386250598567, 0.35672114399690125]
 [0.7604156908506836, 0.7674550530644948, 0.23847720394552546, 0.8294452834399094, 1.0804114475300428, 0.4060703085606619, 0.6888640957065337, 0.3343223619420136, 0.7636075351231503, 0.03724209063173313  …  1.1052247509058286, 1.036939414566552, 0.590016004446043, 1.0983914945208484, 0.2796647248494942, 0.5103533632849603, 1.0769710826307772, 0.6990725769164804, 0.11904315092961476, 0.41837714984868773]
 [1.0053064029594083, 0.7754094885668386, 0.25381901276200136, 0.8727080609515534, 1.2506728414728154, 0.41008832071754375, 0.5034713184150569, 0.19682926309607676, 0.6106487534432536, 0.0715403130063657  …  1.2729780841888372, 1.1972501755009486, 0.592849788111663, 1.2227273732202086, 0.15679756465134848, 0.31151327408613535, 1.2523068815326168, 0.7558685674918, -0.08699326886042917, 0.5490906533963484]
 [1.383919386374002, 0.7837838905457066, 0.2923610032618746, 0.9236774476427589, 1.523181116666089, 0.3997205500175946, 0.2955086213223752, 0.009936882166138456, 0.4342715986146158, 0.1385042704272604  …  1.472789062410583, 1.35732431035077, 0.4905464950458207, 1.3244911127993384, -0.007538677655676027, 0.02053802821052491, 1.4211369646984968, 0.8788149364361183, -0.37901171481451906, 0.7632432858301461]
 [1.9996969576951213, 0.8237528265358053, 0.36398372741091023, 0.9962401787667565, 1.989912694960366, 0.36144105137963456, 0.022407587889207218, -0.27997753933805464, 0.2327477796056976, 0.26713486066052694  …  1.7455467459548604, 1.5225296611969998, 0.19108356494503692, 1.3716391061004376, -0.22617082700556285, -0.4215580154598785, 1.5721740596708869, 1.1317783194275162, -0.8466567366911893, 1.1446223361289465]
```

--------------------------------

### Solver Output and Interpolation

Source: https://docs.sciml.ai/DiffEqGPU/stable/examples/reaction_diffusion

Displays the solver return code, interpolation method, time steps, and state vectors.

```text
retcode: Success
Interpolation: 3rd order Hermite
t: 125-element Vector{Float64}:
   0.0
   2.796123741859731e-5
   6.43785906507212e-5
   0.00013556965396987797
   0.00023376199344619583
   0.00037760221289476507
   0.0005650578854793943
   0.0008071655769861446
   0.001105353327393115
   0.0014666858344513482
   ⋮
   7.9506232854699155
   9.220778312424493
  11.055370462175839
  13.943634271456988
  18.834203937890415
  27.13897752197475
  43.69184525085535
  74.9721648444126
 100.0
u: 125-element Vector{Array{Float64, 3}}:
 [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
 [3.906344697357476e-10 3.909112473032509e-10 … 2.796119589268521e-5 2.7883165317332937e-5; 3.909112473032509e-10 3.911880248707541e-10 … 2.8039337179064475e-5 2.796119589268521e-5; … ; 3.909112473032509e-10 3.911880248707541e-10 … 2.8039337179064475e-5 2.796119589268521e-5; 3.906344697357476e-10 3.909112473032509e-10 … 2.796119589268521e-5 2.7883165317332937e-5;;; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.7961237409320334e-5 2.79612374093249e-5; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.796123740931577e-5 2.7961237409320334e-5; … ; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.796123740931577e-5 2.7961237409320334e-5; 2.7961237415829508e-5 2.7961237415829508e-5 … 2.7961237409320334e-5 2.79612374093249e-5;;; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.796045559984408e-5 2.7960455599839515e-5; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.7960455599848643e-5 2.796045559984408e-5; … ; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.7960455599848643e-5 2.796045559984408e-5; 2.7960455593334905e-5 2.7960455593334905e-5 … 2.796045559984408e-5 2.7960455599839515e-5]
```