diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml
new file mode 100644
index 00000000..3165c49c
--- /dev/null
+++ b/.github/workflows/Test-GPU.yml
@@ -0,0 +1,46 @@
+name: Test-GPU
+
+on:
+  push:
+    branches:
+      - main
+  pull_request:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
+# needed to allow julia-actions/cache to delete old caches that it has created
+permissions:
+  actions: write
+  contents: read
+
+jobs:
+  test:
+    runs-on: self-hosted
+    env:
+      CUDA_VISIBLE_DEVICES: 1
+      JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_smc
+      JULIA_SMC_TEST_GROUP: "GPU"
+    strategy:
+      matrix:
+        julia-version: ['1.10', '1']
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.julia-version }}
+          arch: x64
+      - uses: julia-actions/julia-downgrade-compat@v1
+        if: ${{ matrix.version == '1.10' }}
+        with:
+          skip: LinearAlgebra, Random, SparseArrays
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v5
+        with:
+          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false
diff --git a/Project.toml b/Project.toml
index 33e2f6d9..6656e125 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SparseMatrixColorings"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
 authors = ["Guillaume Dalle", "Alexis Montoison"]
-version = "0.4.20"
+version = "0.4.21"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -12,15 +12,18 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [weakdeps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 
 [extensions]
+SparseMatrixColoringsCUDAExt = "CUDA"
 SparseMatrixColoringsCliqueTreesExt = "CliqueTrees"
 SparseMatrixColoringsColorsExt = "Colors"
 
 [compat]
 ADTypes = "1.2.1"
+CUDA = "5.8.2"
 CliqueTrees = "1"
 Colors = "0.12.11, 0.13"
 DocStringExtensions = "0.8,0.9"
diff --git a/ext/SparseMatrixColoringsCUDAExt.jl b/ext/SparseMatrixColoringsCUDAExt.jl
new file mode 100644
index 00000000..2750964d
--- /dev/null
+++ b/ext/SparseMatrixColoringsCUDAExt.jl
@@ -0,0 +1,155 @@
+module SparseMatrixColoringsCUDAExt
+
+import SparseMatrixColorings as SMC
+using SparseArrays: SparseMatrixCSC, rowvals, nnz, nzrange
+using CUDA: CuVector, CuMatrix
+using CUDA.CUSPARSE: AbstractCuSparseMatrix, CuSparseMatrixCSC, CuSparseMatrixCSR
+
+SMC.matrix_versions(A::AbstractCuSparseMatrix) = (A,)
+
+## Compression (slow, through CPU)
+
+function SMC.compress(
+    A::AbstractCuSparseMatrix, result::SMC.AbstractColoringResult{structure,:column}
+) where {structure}
+    return CuMatrix(SMC.compress(SparseMatrixCSC(A), result))
+end
+
+function SMC.compress(
+    A::AbstractCuSparseMatrix, result::SMC.AbstractColoringResult{structure,:row}
+) where {structure}
+    return CuMatrix(SMC.compress(SparseMatrixCSC(A), result))
+end
+
+## CSC Result
+
+function SMC.ColumnColoringResult(
+    A::CuSparseMatrixCSC, bg::SMC.BipartiteGraph{T}, color::Vector{<:Integer}
+) where {T<:Integer}
+    group = SMC.group_by_color(T, color)
+    compressed_indices = SMC.column_csc_indices(bg, color)
+    additional_info = (; compressed_indices_gpu_csc=CuVector(compressed_indices))
+    return SMC.ColumnColoringResult(
+        A, bg, color, group, compressed_indices, additional_info
+    )
+end
+
+function SMC.RowColoringResult(
+    A::CuSparseMatrixCSC, bg::SMC.BipartiteGraph{T}, color::Vector{<:Integer}
+) where {T<:Integer}
+    group = SMC.group_by_color(T, color)
+    compressed_indices = SMC.row_csc_indices(bg, color)
+    additional_info = (; compressed_indices_gpu_csc=CuVector(compressed_indices))
+    return SMC.RowColoringResult(A, bg, color, group, compressed_indices, additional_info)
+end
+
+function SMC.StarSetColoringResult(
+    A::CuSparseMatrixCSC,
+    ag::SMC.AdjacencyGraph{T},
+    color::Vector{<:Integer},
+    star_set::SMC.StarSet{<:Integer},
+) where {T<:Integer}
+    group = SMC.group_by_color(T, color)
+    compressed_indices = SMC.star_csc_indices(ag, color, star_set)
+    additional_info = (; compressed_indices_gpu_csc=CuVector(compressed_indices))
+    return SMC.StarSetColoringResult(
+        A, ag, color, group, compressed_indices, additional_info
+    )
+end
+
+## CSR Result
+
+function SMC.ColumnColoringResult(
+    A::CuSparseMatrixCSR, bg::SMC.BipartiteGraph{T}, color::Vector{<:Integer}
+) where {T<:Integer}
+    group = SMC.group_by_color(T, color)
+    compressed_indices = SMC.column_csc_indices(bg, color)
+    compressed_indices_csr = SMC.column_csr_indices(bg, color)
+    additional_info = (; compressed_indices_gpu_csr=CuVector(compressed_indices_csr))
+    return SMC.ColumnColoringResult(
+        A, bg, color, group, compressed_indices, additional_info
+    )
+end
+
+function SMC.RowColoringResult(
+    A::CuSparseMatrixCSR, bg::SMC.BipartiteGraph{T}, color::Vector{<:Integer}
+) where {T<:Integer}
+    group = SMC.group_by_color(T, color)
+    compressed_indices = SMC.row_csc_indices(bg, color)
+    compressed_indices_csr = SMC.row_csr_indices(bg, color)
+    additional_info = (; compressed_indices_gpu_csr=CuVector(compressed_indices_csr))
+    return SMC.RowColoringResult(A, bg, color, group, compressed_indices, additional_info)
+end
+
+function SMC.StarSetColoringResult(
+    A::CuSparseMatrixCSR,
+    ag::SMC.AdjacencyGraph{T},
+    color::Vector{<:Integer},
+    star_set::SMC.StarSet{<:Integer},
+) where {T<:Integer}
+    group = SMC.group_by_color(T, color)
+    compressed_indices = SMC.star_csc_indices(ag, color, star_set)
+    additional_info = (; compressed_indices_gpu_csr=CuVector(compressed_indices))
+    return SMC.StarSetColoringResult(
+        A, ag, color, group, compressed_indices, additional_info
+    )
+end
+
+## Decompression
+
+for R in (:ColumnColoringResult, :RowColoringResult)
+    # loop to avoid method ambiguity
+    @eval function SMC.decompress!(
+        A::CuSparseMatrixCSC, B::CuMatrix, result::SMC.$R{<:CuSparseMatrixCSC}
+    )
+        compressed_indices = result.additional_info.compressed_indices_gpu_csc
+        copyto!(A.nzVal, view(B, compressed_indices))
+        return A
+    end
+
+    @eval function SMC.decompress!(
+        A::CuSparseMatrixCSR, B::CuMatrix, result::SMC.$R{<:CuSparseMatrixCSR}
+    )
+        compressed_indices = result.additional_info.compressed_indices_gpu_csr
+        copyto!(A.nzVal, view(B, compressed_indices))
+        return A
+    end
+end
+
+function SMC.decompress!(
+    A::CuSparseMatrixCSC,
+    B::CuMatrix,
+    result::SMC.StarSetColoringResult{<:CuSparseMatrixCSC},
+    uplo::Symbol=:F,
+)
+    if uplo != :F
+        throw(
+            SMC.UnsupportedDecompressionError(
+                "Single-triangle decompression is not supported on GPU matrices"
+            ),
+        )
+    end
+    compressed_indices = result.additional_info.compressed_indices_gpu_csc
+    copyto!(A.nzVal, view(B, compressed_indices))
+    return A
+end
+
+function SMC.decompress!(
+    A::CuSparseMatrixCSR,
+    B::CuMatrix,
+    result::SMC.StarSetColoringResult{<:CuSparseMatrixCSR},
+    uplo::Symbol=:F,
+)
+    if uplo != :F
+        throw(
+            SMC.UnsupportedDecompressionError(
+                "Single-triangle decompression is not supported on GPU matrices"
+            ),
+        )
+    end
+    compressed_indices = result.additional_info.compressed_indices_gpu_csr
+    copyto!(A.nzVal, view(B, compressed_indices))
+    return A
+end
+
+end
diff --git a/src/decompression.jl b/src/decompression.jl
index bc04394c..a3dd471a 100644
--- a/src/decompression.jl
+++ b/src/decompression.jl
@@ -106,6 +106,14 @@ function compress(
     return Br, Bc
 end
 
+struct UnsupportedDecompressionError
+    msg::String
+end
+
+function Base.showerror(io::IO, err::UnsupportedDecompressionError)
+    return print(io, "UnsupportedDecompressionError: $(err.msg)")
+end
+
 """
     decompress(B::AbstractMatrix, result::AbstractColoringResult{_,:column/:row})
     decompress(Br::AbstractMatrix, Bc::AbstractMatrix, result::AbstractColoringResult{_,:bidirectional})
diff --git a/src/graph.jl b/src/graph.jl
index 23331789..9b53a36d 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -100,7 +100,7 @@ end
 Return a [`SparsityPatternCSC`](@ref) corresponding to the matrix `[0 Aᵀ; A 0]`, with a minimum of allocations.
 """
 function bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool)
-    bidirectional_pattern(SparsityPatternCSC(SparseMatrixCSC(A)); symmetric_pattern)
+    return bidirectional_pattern(SparsityPatternCSC(SparseMatrixCSC(A)); symmetric_pattern)
 end
 
 function bidirectional_pattern(S::SparsityPatternCSC{T}; symmetric_pattern::Bool) where {T}
diff --git a/src/matrices.jl b/src/matrices.jl
index 795983c8..eaf79882 100644
--- a/src/matrices.jl
+++ b/src/matrices.jl
@@ -10,7 +10,7 @@ Return various versions of the same matrix:
 
 Used for internal testing.
 """
-function matrix_versions(A)
+function matrix_versions(A::AbstractMatrix)
     A_dense = Matrix(A)
     A_sparse = sparse(A)
     versions = [
diff --git a/src/result.jl b/src/result.jl
index ad2ba348..b5b9cdd2 100644
--- a/src/result.jl
+++ b/src/result.jl
@@ -146,25 +146,38 @@ $TYPEDFIELDS
 - [`AbstractColoringResult`](@ref)
 """
 struct ColumnColoringResult{
-    M<:AbstractMatrix,T<:Integer,G<:BipartiteGraph{T},GT<:AbstractGroups{T}
+    M<:AbstractMatrix,
+    T<:Integer,
+    G<:BipartiteGraph{T},
+    CT<:AbstractVector{T},
+    GT<:AbstractGroups{T},
+    VT<:AbstractVector{T},
+    A,
 } <: AbstractColoringResult{:nonsymmetric,:column,:direct}
     "matrix that was colored"
     A::M
     "bipartite graph that was used for coloring"
     bg::G
     "one integer color for each column or row (depending on `partition`)"
-    color::Vector{T}
+    color::CT
     "color groups for columns or rows (depending on `partition`)"
     group::GT
     "flattened indices mapping the compressed matrix `B` to the uncompressed matrix `A` when `A isa SparseMatrixCSC`. They satisfy `nonzeros(A)[k] = vec(B)[compressed_indices[k]]`"
-    compressed_indices::Vector{T}
+    compressed_indices::VT
+    "optional data used for decompressing into specific matrix types"
+    additional_info::A
 end
 
 function ColumnColoringResult(
     A::AbstractMatrix, bg::BipartiteGraph{T}, color::Vector{<:Integer}
 ) where {T<:Integer}
-    S = bg.S2
     group = group_by_color(T, color)
+    compressed_indices = column_csc_indices(bg, color)
+    return ColumnColoringResult(A, bg, color, group, compressed_indices, nothing)
+end
+
+function column_csc_indices(bg::BipartiteGraph{T}, color::Vector{<:Integer}) where {T}
+    S = bg.S2
     n = size(S, 1)
     rv = rowvals(S)
     compressed_indices = zeros(T, nnz(S))
@@ -176,7 +189,23 @@ function ColumnColoringResult(
             compressed_indices[k] = (c - 1) * n + i
         end
     end
-    return ColumnColoringResult(A, bg, color, group, compressed_indices)
+    return compressed_indices
+end
+
+function column_csr_indices(bg::BipartiteGraph{T}, color::Vector{<:Integer}) where {T}
+    Sᵀ = bg.S1  # CSC storage of transpose(A)
+    n = size(Sᵀ, 2)
+    rv = rowvals(Sᵀ)
+    compressed_indices = zeros(T, nnz(Sᵀ))
+    for i in axes(Sᵀ, 2)
+        for k in nzrange(Sᵀ, i)
+            j = rv[k]
+            c = color[j]
+            # A[i, j] = B[i, c]
+            compressed_indices[k] = (c - 1) * n + i
+        end
+    end
+    return compressed_indices
 end
 
 """
@@ -195,21 +224,33 @@ $TYPEDFIELDS
 - [`AbstractColoringResult`](@ref)
 """
 struct RowColoringResult{
-    M<:AbstractMatrix,T<:Integer,G<:BipartiteGraph{T},GT<:AbstractGroups{T}
+    M<:AbstractMatrix,
+    T<:Integer,
+    G<:BipartiteGraph{T},
+    CT<:AbstractVector{T},
+    GT<:AbstractGroups{T},
+    VT<:AbstractVector{T},
+    A,
 } <: AbstractColoringResult{:nonsymmetric,:row,:direct}
     A::M
     bg::G
-    color::Vector{T}
+    color::CT
     group::GT
-    compressed_indices::Vector{T}
+    compressed_indices::VT
+    additional_info::A
 end
 
 function RowColoringResult(
     A::AbstractMatrix, bg::BipartiteGraph{T}, color::Vector{<:Integer}
 ) where {T<:Integer}
-    S = bg.S2
     group = group_by_color(T, color)
-    C = length(group)  # ncolors
+    compressed_indices = row_csc_indices(bg, color)
+    return RowColoringResult(A, bg, color, group, compressed_indices, nothing)
+end
+
+function row_csc_indices(bg::BipartiteGraph{T}, color::Vector{<:Integer}) where {T}
+    S = bg.S2
+    C = maximum(color)  # ncolors
     rv = rowvals(S)
     compressed_indices = zeros(T, nnz(S))
     for j in axes(S, 2)
@@ -220,7 +261,23 @@ function RowColoringResult(
             compressed_indices[k] = (j - 1) * C + c
         end
     end
-    return RowColoringResult(A, bg, color, group, compressed_indices)
+    return compressed_indices
+end
+
+function row_csr_indices(bg::BipartiteGraph{T}, color::Vector{<:Integer}) where {T}
+    Sᵀ = bg.S1  # CSC storage of transpose(A)
+    C = maximum(color)  # ncolors
+    rv = rowvals(Sᵀ)
+    compressed_indices = zeros(T, nnz(Sᵀ))
+    for i in axes(Sᵀ, 2)
+        for k in nzrange(Sᵀ, i)
+            j = rv[k]
+            c = color[i]
+            # A[i, j] = B[c, j]
+            compressed_indices[k] = (j - 1) * C + c
+        end
+    end
+    return compressed_indices
 end
 
 """
@@ -239,13 +296,20 @@ $TYPEDFIELDS
 - [`AbstractColoringResult`](@ref)
 """
 struct StarSetColoringResult{
-    M<:AbstractMatrix,T<:Integer,G<:AdjacencyGraph{T},GT<:AbstractGroups{T}
+    M<:AbstractMatrix,
+    T<:Integer,
+    G<:AdjacencyGraph{T},
+    CT<:AbstractVector{T},
+    GT<:AbstractGroups{T},
+    VT<:AbstractVector{T},
+    A,
 } <: AbstractColoringResult{:symmetric,:column,:direct}
     A::M
     ag::G
-    color::Vector{T}
+    color::CT
     group::GT
-    compressed_indices::Vector{T}
+    compressed_indices::VT
+    additional_info::A
 end
 
 function StarSetColoringResult(
@@ -254,11 +318,18 @@ function StarSetColoringResult(
     color::Vector{<:Integer},
     star_set::StarSet{<:Integer},
 ) where {T<:Integer}
+    group = group_by_color(T, color)
+    compressed_indices = star_csc_indices(ag, color, star_set)
+    return StarSetColoringResult(A, ag, color, group, compressed_indices, nothing)
+end
+
+function star_csc_indices(
+    ag::AdjacencyGraph{T}, color::Vector{<:Integer}, star_set
+) where {T}
     (; star, hub) = star_set
     S = pattern(ag)
     edge_to_index = edge_indices(ag)
     n = S.n
-    group = group_by_color(T, color)
     rvS = rowvals(S)
     compressed_indices = zeros(T, nnz(S))  # needs to be independent from the storage in the graph, in case the graph gets reused
     for j in axes(S, 2)
@@ -287,8 +358,7 @@ function StarSetColoringResult(
             end
         end
     end
-
-    return StarSetColoringResult(A, ag, color, group, compressed_indices)
+    return compressed_indices
 end
 
 """
diff --git a/test/Project.toml b/test/Project.toml
index 00ee7348..ab595aba 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,6 +6,7 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
 CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
diff --git a/test/cuda.jl b/test/cuda.jl
new file mode 100644
index 00000000..6549b168
--- /dev/null
+++ b/test/cuda.jl
@@ -0,0 +1,62 @@
+using CUDA.CUSPARSE: CuSparseMatrixCSC, CuSparseMatrixCSR
+using LinearAlgebra
+using SparseArrays
+using SparseMatrixColorings
+import SparseMatrixColorings as SMC
+using StableRNGs
+using Test
+
+include("utils.jl")
+
+rng = StableRNG(63)
+
+asymmetric_params = vcat(
+    [(10, 20, p) for p in (0.0:0.2:0.5)],
+    [(20, 10, p) for p in (0.0:0.2:0.5)],
+    [(100, 200, p) for p in (0.01:0.02:0.05)],
+    [(200, 100, p) for p in (0.01:0.02:0.05)],
+)
+
+symmetric_params = vcat(
+    [(10, p) for p in (0.0:0.2:0.5)], #
+    [(100, p) for p in (0.01:0.02:0.05)],
+)
+
+@testset verbose = true "Column coloring & decompression" begin
+    problem = ColoringProblem(; structure=:nonsymmetric, partition=:column)
+    algo = GreedyColoringAlgorithm(; decompression=:direct)
+    @testset for T in (CuSparseMatrixCSC, CuSparseMatrixCSR)
+        @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params
+            A0 = T(sprand(rng, m, n, p))
+            test_coloring_decompression(A0, problem, algo; gpu=true)
+        end
+    end
+end;
+
+@testset verbose = true "Row coloring & decompression" begin
+    problem = ColoringProblem(; structure=:nonsymmetric, partition=:row)
+    algo = GreedyColoringAlgorithm(; decompression=:direct)
+    @testset for T in (CuSparseMatrixCSC, CuSparseMatrixCSR)
+        @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params
+            A0 = T(sprand(rng, m, n, p))
+            test_coloring_decompression(A0, problem, algo; gpu=true)
+        end
+    end
+end;
+
+@testset verbose = true "Symmetric coloring & direct decompression" begin
+    problem = ColoringProblem(; structure=:symmetric, partition=:column)
+    algo = GreedyColoringAlgorithm(; postprocessing=false, decompression=:direct)
+    @testset for T in (CuSparseMatrixCSC, CuSparseMatrixCSR)
+        @testset "$((; n, p))" for (n, p) in symmetric_params
+            A0 = T(sparse(Symmetric(sprand(rng, n, n, p))))
+            test_coloring_decompression(A0, problem, algo; gpu=true)
+        end
+        A0 = T(sparse(Diagonal(ones(10))))
+        result = coloring(A0, problem, algo)
+        B = compress(A0, result)
+        @test_throws SMC.UnsupportedDecompressionError decompress!(
+            similar(A0), B, result, :U
+        )
+    end
+end;
diff --git a/test/result.jl b/test/result.jl
index 661297f0..0bea3353 100644
--- a/test/result.jl
+++ b/test/result.jl
@@ -1,4 +1,4 @@
-using SparseMatrixColorings: group_by_color
+using SparseMatrixColorings: group_by_color, UnsupportedDecompressionError
 using Test
 
 @testset "Group by color" begin
@@ -31,3 +31,8 @@ end
     B = compress(A, coloring(A, problem, algo))
     @test size(B, 1) == 0
 end
+
+@testset "Errors" begin
+    e = SparseMatrixColorings.UnsupportedDecompressionError("hello")
+    @test sprint(showerror, e) == "UnsupportedDecompressionError: hello"
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1ddb8413..81c4ecc9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,77 +11,84 @@ using Colors: Colors
 include("utils.jl")
 
 @testset verbose = true "SparseMatrixColorings" begin
-    @testset verbose = true "Code quality" begin
-        @testset "Aqua" begin
-            Aqua.test_all(SparseMatrixColorings; stale_deps=(; ignore=[:Requires],))
-        end
-        @testset "JET" begin
-            JET.test_package(SparseMatrixColorings; target_defined_modules=true)
-        end
-        @testset "JuliaFormatter" begin
-            @test JuliaFormatter.format(
-                SparseMatrixColorings; verbose=false, overwrite=false
-            )
-        end
-        @testset "Doctests" begin
-            Documenter.doctest(SparseMatrixColorings)
-        end
-    end
-    @testset verbose = true "Internals" begin
-        @testset "Graph" begin
-            include("graph.jl")
-        end
-        @testset "Forest" begin
-            include("forest.jl")
-        end
-        @testset "Order" begin
-            include("order.jl")
-        end
-        @testset "Check" begin
-            include("check.jl")
-        end
-        @testset "Matrices" begin
-            include("matrices.jl")
-        end
-        @testset "Constructors" begin
-            include("constructors.jl")
-        end
-        @testset "Result" begin
-            include("result.jl")
-        end
-        @testset "Constant coloring" begin
-            include("constant.jl")
-        end
-        @testset "ADTypes coloring algorithms" begin
-            include("adtypes.jl")
-        end
-        @testset "Visualization" begin
-            include("show_colors.jl")
-        end
-    end
-    @testset verbose = true "Correctness" begin
-        @testset "Small instances" begin
-            include("small.jl")
-        end
-        @testset "Random instances" begin
-            include("random.jl")
-        end
-        @testset "Structured matrices" begin
-            include("structured.jl")
-        end
-        @testset "Instances with known colorings" begin
-            include("theory.jl")
-        end
-        @testset "SuiteSparse" begin
-            include("suitesparse.jl")
-        end
-    end
-    @testset verbose = true "Performance" begin
-        @testset "Type stability" begin
-            include("type_stability.jl")
-        end
-        @testset "Allocations" begin
-            include("allocations.jl")
+    if get(ENV, "JULIA_SMC_TEST_GROUP", nothing) == "GPU"
+        @testset "CUDA" begin
+            using CUDA
+            include("cuda.jl")
+        end
+    else
+        @testset verbose = true "Code quality" begin
+            @testset "Aqua" begin
+                Aqua.test_all(SparseMatrixColorings; stale_deps=(; ignore=[:Requires],))
+            end
+            @testset "JET" begin
+                JET.test_package(SparseMatrixColorings; target_defined_modules=true)
+            end
+            @testset "JuliaFormatter" begin
+                @test JuliaFormatter.format(
+                    SparseMatrixColorings; verbose=false, overwrite=false
+                )
+            end
+            @testset "Doctests" begin
+                Documenter.doctest(SparseMatrixColorings)
+            end
+        end
+        @testset verbose = true "Internals" begin
+            @testset "Graph" begin
+                include("graph.jl")
+            end
+            @testset "Forest" begin
+                include("forest.jl")
+            end
+            @testset "Order" begin
+                include("order.jl")
+            end
+            @testset "Check" begin
+                include("check.jl")
+            end
+            @testset "Matrices" begin
+                include("matrices.jl")
+            end
+            @testset "Constructors" begin
+                include("constructors.jl")
+            end
+            @testset "Result" begin
+                include("result.jl")
+            end
+            @testset "Constant coloring" begin
+                include("constant.jl")
+            end
+            @testset "ADTypes coloring algorithms" begin
+                include("adtypes.jl")
+            end
+            @testset "Visualization" begin
+                include("show_colors.jl")
+            end
+        end
+        @testset verbose = true "Correctness" begin
+            @testset "Small instances" begin
+                include("small.jl")
+            end
+            @testset "Random instances" begin
+                include("random.jl")
+            end
+            @testset "Structured matrices" begin
+                include("structured.jl")
+            end
+            @testset "Instances with known colorings" begin
+                include("theory.jl")
+            end
+            @testset "SuiteSparse" begin
+                include("suitesparse.jl")
+            end
+        end
+        @testset verbose = true "Performance" begin
+            @testset "Type stability" begin
+                include("type_stability.jl")
+            end
+            @testset "Allocations" begin
+                include("allocations.jl")
+            end
         end
     end
 end
diff --git a/test/utils.jl b/test/utils.jl
index fb3d42c5..bb80f95f 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -22,6 +22,7 @@ function test_coloring_decompression(
     B0=nothing,
     color0=nothing,
     test_fast=false,
+    gpu=false,
 ) where {structure,partition,decompression}
     color_vec = Vector{Int}[]
     @testset "$(typeof(A))" for A in matrix_versions(A0)
@@ -60,6 +61,17 @@ function test_coloring_decompression(
             !isnothing(B0) && @test B == B0
         end
 
+        @testset "Full decompression" begin
+            @test decompress(B, result) ≈ A0
+            @test decompress(B, result) ≈ A0  # check result wasn't modified
+            @test decompress!(respectful_similar(A, eltype(B)), B, result) ≈ A0
+            @test decompress!(respectful_similar(A, eltype(B)), B, result) ≈ A0
+        end
+
+        if gpu
+            continue
+        end
+
         @testset "Recoverability" begin
             # TODO: find tests for recoverability for substitution decompression
             if decompression == :direct
@@ -81,13 +93,6 @@ function test_coloring_decompression(
             end
         end
 
-        @testset "Full decompression" begin
-            @test decompress(B, result) ≈ A0
-            @test decompress(B, result) ≈ A0  # check result wasn't modified
-            @test decompress!(respectful_similar(A, eltype(B)), B, result) ≈ A0
-            @test decompress!(respectful_similar(A, eltype(B)), B, result) ≈ A0
-        end
-
         @testset "Single-color decompression" begin
             if decompression == :direct  # TODO: implement for :substitution too
                 A2 = respectful_similar(A, eltype(B))
@@ -194,11 +199,6 @@ function test_bicoloring_decompression(
             end
         end
 
-        if decompression == :direct
-            @testset "Recoverability" begin
-                @test structurally_biorthogonal(A0, row_color, column_color)
-            end
-        end
         @testset "Full decompression" begin
             @test decompress(Br, Bc, result) ≈ A0
             @test decompress(Br, Bc, result) ≈ A0  # check result wasn't modified
@@ -209,6 +209,12 @@ function test_bicoloring_decompression(
                 respectful_similar(A, promote_eltype(Br, Bc)), Br, Bc, result
             ) ≈ A0
         end
+
+        if decompression == :direct
+            @testset "Recoverability" begin
+                @test structurally_biorthogonal(A0, row_color, column_color)
+            end
+        end
     end
 end