diff --git a/README.md b/README.md
index c7b5dc1..e2fe32b 100644
--- a/README.md
+++ b/README.md
@@ -40,6 +40,14 @@ Under Linux: Make the file executable with `chmod +x tgv_qsm.jl` and run directl
 <folder>/tgv_qsm.jl --help
 ```
 
+### Number of threads
+
+In case Julia uses only 1 CPU thread, you can use the command
+
+```bash
+julia --threads=auto <folder>/tgv_qsm.jl <arguments>
+```
+
 ## Run in Julia
 
 ### Setup
diff --git a/tgv_qsm.jl b/tgv_qsm.jl
index d3f0e5b..d560a71 100644
--- a/tgv_qsm.jl
+++ b/tgv_qsm.jl
@@ -6,7 +6,7 @@
 # On windows use: `julia --threads=auto <path-to-file>/tgv_qsm.jl ARGS`
 
 # Example call:
-# `./tgv_qsm.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz --no-gpu
+# `./tgv_qsm.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz
 
 import Pkg
 
@@ -31,7 +31,7 @@ Comonicon.get_version(::Module) = version
     mask = niread(fn_mask) .!= 0
     res = header(phase).pixdim[2:4]
     println("Resolution from NIfTI header [mm]: $(round.(Float64.(res); digits=2))")
-    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=!no_gpu, step_size, type, nblocks)
+    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, step_size, type, nblocks)
     println("Writing output")
     savenii(chi, output; header=header(phase))
 end
diff --git a/tgv_qsm_amdgpu.jl b/tgv_qsm_amdgpu.jl
new file mode 100644
index 0000000..9d44d05
--- /dev/null
+++ b/tgv_qsm_amdgpu.jl
@@ -0,0 +1,37 @@
+#!/usr/bin/env -S julia --color=yes --startup-file=no --threads=auto
+
+## Usage
+
+# Call with: `<path-to-file>/tgv_qsm_amdgpu.jl ARGS`
+# On windows use: `julia --threads=auto <path-to-file>/tgv_qsm_amdgpu.jl ARGS`
+
+# Example call:
+# `./tgv_qsm_amdgpu.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz
+
+import Pkg
+
+## Uncomment to use a local julia package directory instead of the global one
+# package_dir = joinpath(@__DIR__, ".tgv_cmd_packages")
+# mkpath(package_dir)
+# Pkg.activate(package_dir)
+
+try
+    using AMDGPU, QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+catch
+    Pkg.add(["AMDGPU", "QuantitativeSusceptibilityMappingTGV", "MriResearchTools", "Comonicon"])
+    using AMDGPU, QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+end
+
+version = Comonicon.get_version(QuantitativeSusceptibilityMappingTGV)
+Comonicon.get_version(::Module) = version
+
+@main function tgv_qsm(fn_phase, fn_mask; TE::Float64, output::String="output.nii.gz", fieldstrength::Float64=3.0, regularization::Float64=2.0, erosions::Int=3, B0_dir::Array{Int}=[0,0,1], dedimensionalize::Bool=false, no_laplacian_correction::Bool=false, step_size::Float64=3.0, type::DataType=Float32, nblocks::Int=32)
+    println("Starting calculation...")
+    phase = readphase(fn_phase)
+    mask = niread(fn_mask) .!= 0
+    res = header(phase).pixdim[2:4]
+    println("Resolution from NIfTI header [mm]: $(round.(Float64.(res); digits=2))")
+    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=AMDGPU, step_size, type, nblocks)
+    println("Writing output")
+    savenii(chi, output; header=header(phase))
+end
diff --git a/tgv_qsm_cuda.jl b/tgv_qsm_cuda.jl
index 5745afe..71fd942 100644
--- a/tgv_qsm_cuda.jl
+++ b/tgv_qsm_cuda.jl
@@ -2,11 +2,11 @@
 
 ## Usage
 
-# Call with: `<path-to-file>/tgv_qsm.jl ARGS`
-# On windows use: `julia --threads=auto <path-to-file>/tgv_qsm.jl ARGS`
+# Call with: `<path-to-file>/tgv_qsm_cuda.jl ARGS`
+# On windows use: `julia --threads=auto <path-to-file>/tgv_qsm_cuda.jl ARGS`
 
 # Example call:
-# `./tgv_qsm.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz --no-gpu
+# `./tgv_qsm_cuda.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz
 
 import Pkg
 
@@ -31,7 +31,7 @@ Comonicon.get_version(::Module) = version
     mask = niread(fn_mask) .!= 0
     res = header(phase).pixdim[2:4]
     println("Resolution from NIfTI header [mm]: $(round.(Float64.(res); digits=2))")
-    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=!no_gpu, step_size, type, nblocks, gpu=CUDA)
+    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=CUDA, step_size, type, nblocks)
     println("Writing output")
     savenii(chi, output; header=header(phase))
 end
diff --git a/tgv_qsm_metal.jl b/tgv_qsm_metal.jl
new file mode 100644
index 0000000..38a8c0d
--- /dev/null
+++ b/tgv_qsm_metal.jl
@@ -0,0 +1,37 @@
+#!/usr/bin/env -S julia --color=yes --startup-file=no --threads=auto
+
+## Usage
+
+# Call with: `<path-to-file>/tgv_qsm_metal.jl ARGS`
+# On windows use: `julia --threads=auto <path-to-file>/tgv_qsm_metal.jl ARGS`
+
+# Example call:
+# `./tgv_qsm_metal.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz
+
+import Pkg
+
+## Uncomment to use a local julia package directory instead of the global one
+# package_dir = joinpath(@__DIR__, ".tgv_cmd_packages")
+# mkpath(package_dir)
+# Pkg.activate(package_dir)
+
+try
+    using Metal, QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+catch
+    Pkg.add(["Metal", "QuantitativeSusceptibilityMappingTGV", "MriResearchTools", "Comonicon"])
+    using Metal, QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+end
+
+version = Comonicon.get_version(QuantitativeSusceptibilityMappingTGV)
+Comonicon.get_version(::Module) = version
+
+@main function tgv_qsm(fn_phase, fn_mask; TE::Float64, output::String="output.nii.gz", fieldstrength::Float64=3.0, regularization::Float64=2.0, erosions::Int=3, B0_dir::Array{Int}=[0,0,1], dedimensionalize::Bool=false, no_laplacian_correction::Bool=false, step_size::Float64=3.0, type::DataType=Float32, nblocks::Int=32)
+    println("Starting calculation...")
+    phase = readphase(fn_phase)
+    mask = niread(fn_mask) .!= 0
+    res = header(phase).pixdim[2:4]
+    println("Resolution from NIfTI header [mm]: $(round.(Float64.(res); digits=2))")
+    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=Metal, step_size, type, nblocks)
+    println("Writing output")
+    savenii(chi, output; header=header(phase))
+end
diff --git a/tgv_qsm_oneapi.jl b/tgv_qsm_oneapi.jl
new file mode 100644
index 0000000..a12ed46
--- /dev/null
+++ b/tgv_qsm_oneapi.jl
@@ -0,0 +1,37 @@
+#!/usr/bin/env -S julia --color=yes --startup-file=no --threads=auto
+
+## Usage
+
+# Call with: `<path-to-file>/tgv_qsm_oneapi.jl ARGS`
+# On windows use: `julia --threads=auto <path-to-file>/tgv_qsm_oneapi.jl ARGS`
+
+# Example call:
+# `./tgv_qsm_oneapi.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz
+
+import Pkg
+
+## Uncomment to use a local julia package directory instead of the global one
+# package_dir = joinpath(@__DIR__, ".tgv_cmd_packages")
+# mkpath(package_dir)
+# Pkg.activate(package_dir)
+
+try
+    using oneAPI, QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+catch
+    Pkg.add(["oneAPI", "QuantitativeSusceptibilityMappingTGV", "MriResearchTools", "Comonicon"])
+    using oneAPI, QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+end
+
+version = Comonicon.get_version(QuantitativeSusceptibilityMappingTGV)
+Comonicon.get_version(::Module) = version
+
+@main function tgv_qsm(fn_phase, fn_mask; TE::Float64, output::String="output.nii.gz", fieldstrength::Float64=3.0, regularization::Float64=2.0, erosions::Int=3, B0_dir::Array{Int}=[0,0,1], dedimensionalize::Bool=false, no_laplacian_correction::Bool=false, step_size::Float64=3.0, type::DataType=Float32, nblocks::Int=32)
+    println("Starting calculation...")
+    phase = readphase(fn_phase)
+    mask = niread(fn_mask) .!= 0
+    res = header(phase).pixdim[2:4]
+    println("Resolution from NIfTI header [mm]: $(round.(Float64.(res); digits=2))")
+    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=oneAPI, step_size, type, nblocks)
+    println("Writing output")
+    savenii(chi, output; header=header(phase))
+end