# GPU Parallel Computing for Global Sensitivity Analysis

**URL:** <https://discourse.julialang.org/t/gpu-parallel-computing-for-global-sensitivity-analysis/101545>\
**Category:** GPU\
**Created:** [July 12, 2023, 4:57pm UTC](https://discourse.julialang.org/t/gpu-parallel-computing-for-global-sensitivity-analysis/101545 "2023-07-12T16:57:19Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [July 12, 2023, 4:57pm UTC](https://discourse.julialang.org/t/gpu-parallel-computing-for-global-sensitivity-analysis/101545/1 "2023-07-12T16:57:19Z")

</div>

Hi everyone. I try to speed up Global Sensitivity Analysis calculations based on the code in this link ([Global sensitivity analysis in a large parameter space - computational cost - #11 by TS-CUBED](https://discourse.julialang.org/t/global-sensitivity-analysis-in-a-large-parameter-space-computational-cost/93598/11)). Essentially, I changed the

> sol = solve(ensemble\_prob,GPUTsit5(), reltol=1e-8, abstol=1e-8, EnsembleThreads();saveat = saveat,trajectories=chunksize)

to

> sol = solve(ensemble\_prob,GPUTsit5(), reltol=1e-8, abstol=1e-8, EnsembleGPUArray(CUDA.CUDABackend());saveat = saveat,trajectories=chunksize)

However, I get the following error:

PosDefException: matrix is not positive definite; Cholesky factorization failed

Does anyone know what does it happen? Also, does the Threads.@threads macro work on GPUs?

Thank you! 🙂

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 13, 2023, 10:51am UTC](https://discourse.julialang.org/t/gpu-parallel-computing-for-global-sensitivity-analysis/101545/2 "2023-07-13T10:51:43Z")

</div>

LOdd, there shouldn’t be a factorization anywhere. Please share an MWE.

> [@Nick1](#):
>
> does the Threads.@threads macro work on GPUs?

No.

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [July 14, 2023, 8:24pm UTC](https://discourse.julialang.org/t/gpu-parallel-computing-for-global-sensitivity-analysis/101545/3 "2023-07-14T20:24:31Z")

</div>

I share a toy model that I use for this purpose. The first is the function to define the ODEs, the second is the main code.

```julia
function pk_model!(du, u, p, t)

    # PROBLEM
    # Model: `pk1`
    # - a simple one compartment pk mode
    
    # PARAMETERS
    CL = p[1]; # Clearance (volume/time)
    V = p[2]; # Central volume (volume)
    KA = p[3]; # Absorption rate constant 1 (1/time)
    
    # CMT
    EV = u[1]; # extravascular compartment (mass)
    CENT = u[2]; # Central compartment (mass)
    
    #define Central compartment concentration
    CP = CENT/V;
    
    # ODEs
    dxdt_EV = -KA *EV;
    dxdt_CENT = KA *EV - CL*CP;
    
    # Function outputs
    du[1] = dxdt_EV;
    du[2] = dxdt_CENT

    # Auxiliary output to calculate AUC, for sensitivity analysis
    dxdt_AUC = CENT;
    du[3] = dxdt_AUC;

end

```

```julia
using DifferentialEquations
using GlobalSensitivity, QuasiMonteCarlo, Statistics, Random
using DiffEqGPU, CUDA
using JLD2
using ComponentArrays
using OrderedCollections

include("pk_model_sens.jl")

Dose_mg = 100.0f0;
u0 = [Dose_mg, 0.0f0, 0.0f0];
param = ComponentArray(CL = 1.0f0, V = 20.0f0, KA = 1.0f0);

#Define the simulation timespan
tspan = (0.f0, 20.f0); #hr
t = collect(range(0.0f0, stop=20.0f0, length=20_0));

# Define an ode problem
# note; the pk_model! is a function defined in 'pk_model.jl'
prob = ODEProblem(pk_model!, u0, tspan, collect(param));

#define parameters to be explored and their ranges
params = OrderedDict(:CL=>[0.1f0, 10.f0], :KA=>[0.01f0,10.0f0]);
param_names = collect(keys(params));
param_values = collect(values(params));

Random.seed!(1235) 
# number of samples
N = 10_000;
# lower and upper bound for model parameters
lb = reduce(hcat,param_values)[1,:];
ub = reduce(hcat,param_values)[2,:];

# Generate design matrices based on QMC
sampler = SobolSample();
A,B = QuasiMonteCarlo.generate_design_matrices(N,lb,ub,sampler);

### GSA using chunks
ensemble_chunked = function (p)
    # global parameterset = p
   
    Np::Int64 = size(p,2)
    chunks::Int64 = 10_0;
# println("Running ", Np, " cases.")
  
    chunksize::Int64 = Np/chunks
  
# println("Using ", chunks, " chunk(s) of size ", chunksize, ".")
  
  
    out = zeros(2,Np)
  
    for k in 1:chunks
        offset = Int((k - 1) * Np/chunks)
        startindx = Int(offset + 1)
        endindx = Int(k * Np/chunks)
  
# println("Starting chunk ", k, ", from ", startindx, " to ", endindx, " with offset ", offset, ".")
  
        pchunk = p[:,startindx:endindx]
  
       # function prob_func(prob,i,repeat)
       # remake(prob; p=pchunk[:,i])
       # end
       function prob_func(prob,i,repeat) 
          p_pop = deepcopy(param);
          p_pop[param_names] =pchunk[:,i];
          remake(prob;p=collect(p_pop))
        end
  
        ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
  
        ensemble_sol = solve(ensemble_prob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend());saveat = t,trajectories=chunksize)
  
         Threads.@threads for i in 1:chunksize
            out[1,i+offset] = ensemble_sol[i][3,end]
            out[2,i+offset] = maximum(ensemble_sol[i][2,:])
        end
    end
    out
end

sobol_chunks = gsa(ensemble_chunked,Sobol(),A,B,batch=true);
save_object("chunks_GPU.jld2", (sobol_chunks=sobol_chunks, param_names=param_names))

sobol_chunks, param_names = load_object("chunks_GPU.jld2");

```