# Solve linear systems inside CUDA kernel function
**URL:** https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202
**Category:** GPU
**Created:** [February 14, 2024, 8:52am UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202 "2024-02-14T08:52:59Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@FR13ndSDP](https://discourse.julialang.org/u/FR13ndSDP)
#### Post date: [February 14, 2024, 8:52am UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/1 "2024-02-14T08:52:59Z")
Hello everyone, I was wondering how to solve a small linear system Ax = b inside a CUDA kernel function, where A, x and b are `MMatrix` and `MVector`. I want to do this with `x = A\b`, but it does not work inside a kernel, is there a solution to this problem?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [February 14, 2024, 11:14am UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/2 "2024-02-14T11:14:52Z")
If you use static arrays it won’t be an issue. This is done in DiffEqGPU.jl
---
### Author:  [@FR13ndSDP](https://discourse.julialang.org/u/FR13ndSDP)
#### Post date: [February 14, 2024, 11:47am UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/3 "2024-02-14T11:47:59Z")
Hi Chris, unfortunately I haven’t been able to solve this problem yet. I hope this code snippet illustrates the problem:
```julia
using CUDA, StaticArrays
const N = 20
function test()
i = (blockIdx().x-1)* blockDim().x + threadIdx().x
# each thread has unique A and b
a = @MMatrix rand(Float64, N, N)
b = @MVector rand(Float64, N)
c = MVector{N, Float64}(undef)
# This works
c = a * b
# But this does not
c = a \ b
return
end
@cuda threads=10 test()
```
Also, is it possible to do batched small linear systems solving with CUBLAS or CUSOLVER?
---
### Author:  [@Zentrik](https://discourse.julialang.org/u/Zentrik)
#### Post date: [February 14, 2024, 12:07pm UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/4 "2024-02-14T12:07:36Z")
You want to use a SMatrix and Svector not a MMatrix, Mvector I assume.
---
### Author:  [@FR13ndSDP](https://discourse.julialang.org/u/FR13ndSDP)
#### Post date: [February 14, 2024, 1:47pm UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/5 "2024-02-14T13:47:20Z")
In my case, A and b will be constructed inside the kernel, so they have to be mutable. Moreover, use `SMatrix` and `SVector` does not help in the snippet above.
---
### Author:  [@trahflow](https://discourse.julialang.org/u/trahflow)
#### Post date: [February 14, 2024, 1:58pm UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/6 "2024-02-14T13:58:08Z")
> [@FR13ndSDP](#):
>
> In my case, AAA and bbb will be constructed inside the kernel, so they have to be mutable.
That’s very unlikely to work. You cannot dynamically allocate memory inside a GPU kernel (see also this recent post: [Modifying a thread-local vector within CUDA Dynamic Parallelism - #2 by vchuravy](https://discourse.julialang.org/t/modifying-a-thread-local-vector-within-cuda-dynamic-parallelism/110110/2)).
What should work though is to allocate all CuArrays outside the kernel, then inside the kernel convert the relevant `view`s into your arrays into `SMatrix`/`SVector`s and do the solve on StaticArrays only. (I don’t have access to a GPU atm to check)
---
### Author:  [@FR13ndSDP](https://discourse.julialang.org/u/FR13ndSDP)
#### Post date: [February 14, 2024, 2:55pm UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/7 "2024-02-14T14:55:12Z")
Allocate memory with `MMatrix` and `MVector` works fine, I think the problem is the `\` operation has some allocations, I’ve tried to implement the Gauss elimination method to solve linear equations, it works well on GPU, but I’m worried about its performance.
---
### Author:  [@utkarsh530](https://discourse.julialang.org/u/utkarsh530)
#### Post date: [February 14, 2024, 4:21pm UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/8 "2024-02-14T16:21:27Z")
Hi, can you try for some `N < 14`? IIUC, there were some allocations here: [StaticArrays.jl/src/solve.jl at master · JuliaArrays/StaticArrays.jl · GitHub](https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/solve.jl#L73)
We can probably try to get that dispatch setup in LinearSolve.jl but not sure as the previous approach may be done for performance reasons.
---
### Author:  [@FR13ndSDP](https://discourse.julialang.org/u/FR13ndSDP)
#### Post date: [February 14, 2024, 4:54pm UTC](https://discourse.julialang.org/t/solve-linear-systems-inside-cuda-kernel-function/110202/9 "2024-02-14T16:54:12Z")
Yes, you are right! For N \leq 14 it works well. But in my case, the typical size is N=[20,200]. I think with my implementation of Gauss elimination will be faster than `A\b` if N \leq 14.