diff --git a/benchmarks/matrix_benchmarks.jl b/benchmarks/matrix_benchmarks.jl
index 9f1a1c5..725b7b1 100644
--- a/benchmarks/matrix_benchmarks.jl
+++ b/benchmarks/matrix_benchmarks.jl
@@ -227,3 +227,68 @@ function benchmark_sparse_dense_add!(
 
     return nothing
 end
+
+"""
+    benchmark_kron!(SUITE, array_constructor, array_type_name; N=100, T=Float64)
+
+Benchmark Kronecker product (kron) for CSC, CSR, and COO formats.
+
+# Arguments
+- `SUITE`: The BenchmarkGroup to add benchmarks to
+- `array_constructor`: Function to construct arrays (e.g., `Array`, `JLArray`)
+- `array_type_name`: String name for the array type (for display)
+
+# Keyword Arguments
+- `N`: Size of the matrices (default: 100)
+- `T`: Element type (default: Float64)
+"""
+function benchmark_kron!(
+    SUITE,
+    array_constructor,
+    array_type_name;
+    N = 100,
+    T = Float64,
+)
+    # Create sparse matrices with 1% density (smaller matrices since kron grows quadratically)
+    sm_a_std = sprand(T, N, N, 0.01)
+    sm_b_std = sprand(T, N, N, 0.01)
+
+    # Convert to different formats
+    sm_a_csc = DeviceSparseMatrixCSC(sm_a_std)
+    sm_b_csc = DeviceSparseMatrixCSC(sm_b_std)
+    
+    sm_a_csr = DeviceSparseMatrixCSR(sm_a_std)
+    sm_b_csr = DeviceSparseMatrixCSR(sm_b_std)
+    
+    sm_a_coo = DeviceSparseMatrixCOO(sm_a_std)
+    sm_b_coo = DeviceSparseMatrixCOO(sm_b_std)
+
+    # Adapt to device
+    dsm_a_csc = adapt(array_constructor, sm_a_csc)
+    dsm_b_csc = adapt(array_constructor, sm_b_csc)
+    
+    dsm_a_csr = adapt(array_constructor, sm_a_csr)
+    dsm_b_csr = adapt(array_constructor, sm_b_csr)
+    
+    dsm_a_coo = adapt(array_constructor, sm_a_coo)
+    dsm_b_coo = adapt(array_constructor, sm_b_coo)
+
+    # Level 3: Format (CSC, CSR, COO - will be plotted together)
+    SUITE["Kronecker Product"][array_type_name]["CSC"] = @benchmarkable begin
+        kron($dsm_a_csc, $dsm_b_csc)
+        _synchronize_backend($dsm_a_csc)
+    end
+
+    SUITE["Kronecker Product"][array_type_name]["CSR"] = @benchmarkable begin
+        kron($dsm_a_csr, $dsm_b_csr)
+        _synchronize_backend($dsm_a_csr)
+    end
+
+    SUITE["Kronecker Product"][array_type_name]["COO"] = @benchmarkable begin
+        kron($dsm_a_coo, $dsm_b_coo)
+        _synchronize_backend($dsm_a_coo)
+    end
+
+    return nothing
+end
+
diff --git a/benchmarks/runbenchmarks.jl b/benchmarks/runbenchmarks.jl
index 7aea965..ca78d25 100644
--- a/benchmarks/runbenchmarks.jl
+++ b/benchmarks/runbenchmarks.jl
@@ -26,17 +26,19 @@ benchmark_matrix_vector_mul!(SUITE, Array, "Array")
 benchmark_matrix_matrix_mul!(SUITE, Array, "Array")
 benchmark_three_arg_dot!(SUITE, Array, "Array")
 benchmark_sparse_dense_add!(SUITE, Array, "Array")
+benchmark_kron!(SUITE, Array, "Array")
 benchmark_conversions!(SUITE, Array, "Array")
 
 # Run benchmarks for JLArrays
 println("Running benchmarks for JLArrays...")
-benchmark_vector_sum!(SUITE, jl, "JLArray")
-benchmark_vector_sparse_dense_dot!(SUITE, jl, "JLArray")
-benchmark_matrix_vector_mul!(SUITE, jl, "JLArray")
-benchmark_matrix_matrix_mul!(SUITE, jl, "JLArray")
-benchmark_three_arg_dot!(SUITE, jl, "JLArray")
-benchmark_sparse_dense_add!(SUITE, jl, "JLArray")
-benchmark_conversions!(SUITE, jl, "JLArray")
+benchmark_vector_sum!(SUITE, JLArray, "JLArray")
+benchmark_vector_sparse_dense_dot!(SUITE, JLArray, "JLArray")
+benchmark_matrix_vector_mul!(SUITE, JLArray, "JLArray")
+benchmark_matrix_matrix_mul!(SUITE, JLArray, "JLArray")
+benchmark_three_arg_dot!(SUITE, JLArray, "JLArray")
+benchmark_sparse_dense_add!(SUITE, JLArray, "JLArray")
+benchmark_kron!(SUITE, JLArray, "JLArray")
+benchmark_conversions!(SUITE, JLArray, "JLArray")
 
 # Tune and run benchmarks
 println("\nTuning benchmarks...")
diff --git a/src/DeviceSparseArrays.jl b/src/DeviceSparseArrays.jl
index 17d96da..2c89380 100644
--- a/src/DeviceSparseArrays.jl
+++ b/src/DeviceSparseArrays.jl
@@ -1,7 +1,7 @@
 module DeviceSparseArrays
 
 using LinearAlgebra
-import LinearAlgebra: wrap, copymutable_oftype, __normalize!
+import LinearAlgebra: wrap, copymutable_oftype, __normalize!, kron
 using SparseArrays
 import SparseArrays: SparseVector, SparseMatrixCSC
 import SparseArrays: getcolptr, getrowval, getnzval, nonzeroinds
diff --git a/src/conversions/conversion_kernels.jl b/src/conversions/conversion_kernels.jl
index 39c21aa..22c96a4 100644
--- a/src/conversions/conversion_kernels.jl
+++ b/src/conversions/conversion_kernels.jl
@@ -39,10 +39,10 @@ end
     keys,
     @Const(rowind),
     @Const(colind),
-    @Const(n),
+    @Const(m),  # Number of rows - needed for proper column-major lexicographic ordering
 )
     i = @index(Global)
-    keys[i] = colind[i] * n + rowind[i]
+    keys[i] = colind[i] * m + rowind[i]
 end
 
 # Kernel for creating sort keys for COO → CSR conversion
@@ -50,10 +50,10 @@ end
     keys,
     @Const(rowind),
     @Const(colind),
-    @Const(m),
+    @Const(n),
 )
     i = @index(Global)
-    keys[i] = rowind[i] * m + colind[i]
+    keys[i] = rowind[i] * n + colind[i]
 end
 
 # Kernel for counting entries per column (for COO → CSC)
diff --git a/src/conversions/conversions.jl b/src/conversions/conversions.jl
index 8b254bf..9aafb61 100644
--- a/src/conversions/conversions.jl
+++ b/src/conversions/conversions.jl
@@ -5,26 +5,6 @@
 # CSC ↔ COO Conversions
 # ============================================================================
 
-"""
-    DeviceSparseMatrixCOO(A::DeviceSparseMatrixCSC)
-
-Convert a Compressed Sparse Column (CSC) matrix to Coordinate (COO) format.
-
-The conversion preserves all matrix data and maintains backend compatibility.
-The result will be on the same backend (CPU/GPU) as the input matrix.
-
-# Examples
-```julia
-using DeviceSparseArrays, SparseArrays
-
-# Create a CSC matrix
-A_sparse = sparse([1, 2, 3], [1, 2, 3], [1.0, 2.0, 3.0], 3, 3)
-A_csc = DeviceSparseMatrixCSC(A_sparse)
-
-# Convert to COO format
-A_coo = DeviceSparseMatrixCOO(A_csc)
-```
-"""
 function DeviceSparseMatrixCOO(A::DeviceSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
     m, n = size(A)
     nnz_count = nnz(A)
@@ -43,27 +23,6 @@ function DeviceSparseMatrixCOO(A::DeviceSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
     return DeviceSparseMatrixCOO(m, n, rowind, colind, nzval)
 end
 
-"""
-    DeviceSparseMatrixCSC(A::DeviceSparseMatrixCOO)
-
-Convert a Coordinate (COO) matrix to Compressed Sparse Column (CSC) format.
-
-The conversion sorts the COO entries by column (then by row within each column)
-and builds the column pointer structure. The result maintains backend compatibility
-with the input matrix.
-
-# Examples
-```julia
-using DeviceSparseArrays, SparseArrays
-
-# Create a COO matrix
-A_sparse = sparse([1, 2, 3], [1, 2, 3], [1.0, 2.0, 3.0], 3, 3)
-A_coo = DeviceSparseMatrixCOO(A_sparse)
-
-# Convert to CSC format
-A_csc = DeviceSparseMatrixCSC(A_coo)
-```
-"""
 function DeviceSparseMatrixCSC(A::DeviceSparseMatrixCOO{Tv,Ti}) where {Tv,Ti}
     m, n = size(A)
     nnz_count = nnz(A)
@@ -75,7 +34,7 @@ function DeviceSparseMatrixCSC(A::DeviceSparseMatrixCOO{Tv,Ti}) where {Tv,Ti}
 
     # Create keys on device
     kernel! = kernel_make_csc_keys!(backend)
-    kernel!(keys, A.rowind, A.colind, n; ndrange = (nnz_count,))
+    kernel!(keys, A.rowind, A.colind, m; ndrange = (nnz_count,))
 
     perm = AcceleratedKernels.sortperm(keys)
 
@@ -103,26 +62,6 @@ end
 # CSR ↔ COO Conversions
 # ============================================================================
 
-"""
-    DeviceSparseMatrixCOO(A::DeviceSparseMatrixCSR)
-
-Convert a Compressed Sparse Row (CSR) matrix to Coordinate (COO) format.
-
-The conversion preserves all matrix data and maintains backend compatibility.
-The result will be on the same backend (CPU/GPU) as the input matrix.
-
-# Examples
-```julia
-using DeviceSparseArrays, SparseArrays
-
-# Create a CSR matrix
-A_sparse = sparse([1, 2, 3], [1, 2, 3], [1.0, 2.0, 3.0], 3, 3)
-A_csr = DeviceSparseMatrixCSR(A_sparse)
-
-# Convert to COO format
-A_coo = DeviceSparseMatrixCOO(A_csr)
-```
-"""
 function DeviceSparseMatrixCOO(A::DeviceSparseMatrixCSR{Tv,Ti}) where {Tv,Ti}
     m, n = size(A)
     nnz_count = nnz(A)
@@ -141,27 +80,6 @@ function DeviceSparseMatrixCOO(A::DeviceSparseMatrixCSR{Tv,Ti}) where {Tv,Ti}
     return DeviceSparseMatrixCOO(m, n, rowind, colind, nzval)
 end
 
-"""
-    DeviceSparseMatrixCSR(A::DeviceSparseMatrixCOO)
-
-Convert a Coordinate (COO) matrix to Compressed Sparse Row (CSR) format.
-
-The conversion sorts the COO entries by row (then by column within each row)
-and builds the row pointer structure. The result maintains backend compatibility
-with the input matrix.
-
-# Examples
-```julia
-using DeviceSparseArrays, SparseArrays
-
-# Create a COO matrix
-A_sparse = sparse([1, 2, 3], [1, 2, 3], [1.0, 2.0, 3.0], 3, 3)
-A_coo = DeviceSparseMatrixCOO(A_sparse)
-
-# Convert to CSR format
-A_csr = DeviceSparseMatrixCSR(A_coo)
-```
-"""
 function DeviceSparseMatrixCSR(A::DeviceSparseMatrixCOO{Tv,Ti}) where {Tv,Ti}
     m, n = size(A)
     nnz_count = nnz(A)
@@ -173,7 +91,7 @@ function DeviceSparseMatrixCSR(A::DeviceSparseMatrixCOO{Tv,Ti}) where {Tv,Ti}
 
     # Create keys on device
     kernel! = kernel_make_csr_keys!(backend)
-    kernel!(keys, A.rowind, A.colind, m; ndrange = (nnz_count,))
+    kernel!(keys, A.rowind, A.colind, n; ndrange = (nnz_count,))
 
     # Sort - use AcceleratedKernels
     perm = AcceleratedKernels.sortperm(keys)
diff --git a/src/matrix_coo/matrix_coo.jl b/src/matrix_coo/matrix_coo.jl
index f90b67b..4d6567c 100644
--- a/src/matrix_coo/matrix_coo.jl
+++ b/src/matrix_coo/matrix_coo.jl
@@ -331,3 +331,72 @@ function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCOO)
 
     return C
 end
+
+"""
+    kron(A::DeviceSparseMatrixCOO, B::DeviceSparseMatrixCOO)
+
+Compute the Kronecker product of two sparse matrices in COO format.
+
+The Kronecker product of two matrices `A` (size m×n) and `B` (size p×q) 
+is an (m*p)×(n*q) matrix formed by multiplying each element of `A` by the 
+entire matrix `B`.
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A = DeviceSparseMatrixCOO(sparse([1, 2], [1, 2], [1.0, 2.0], 2, 2));
+
+julia> B = DeviceSparseMatrixCOO(sparse([1, 2], [1, 2], [3.0, 4.0], 2, 2));
+
+julia> C = kron(A, B);
+
+julia> size(C)
+(4, 4)
+
+julia> nnz(C)
+4
+```
+"""
+function LinearAlgebra.kron(
+    A::DeviceSparseMatrixCOO{Tv1,Ti1},
+    B::DeviceSparseMatrixCOO{Tv2,Ti2},
+) where {Tv1,Ti1,Tv2,Ti2}
+    # Result dimensions
+    m_C = size(A, 1) * size(B, 1)
+    n_C = size(A, 2) * size(B, 2)
+    nnz_C = nnz(A) * nnz(B)
+
+    # Determine result types
+    Tv = promote_type(Tv1, Tv2)
+    Ti = promote_type(Ti1, Ti2)
+
+    # Check backend compatibility
+    backend_A = get_backend(A)
+    backend_B = get_backend(B)
+    backend_A == backend_B || throw(ArgumentError("Both arrays must have the same backend"))
+
+    # Allocate output arrays
+    rowind_C = similar(A.rowind, Ti, nnz_C)
+    colind_C = similar(A.colind, Ti, nnz_C)
+    nzval_C = similar(A.nzval, Tv, nnz_C)
+
+    # Launch kernel
+    kernel! = kernel_kron_coo!(backend_A)
+    kernel!(
+        A.rowind,
+        A.colind,
+        A.nzval,
+        B.rowind,
+        B.colind,
+        B.nzval,
+        rowind_C,
+        colind_C,
+        nzval_C,
+        size(B, 1),
+        size(B, 2);
+        ndrange = nnz_C,
+    )
+
+    return DeviceSparseMatrixCOO(m_C, n_C, rowind_C, colind_C, nzval_C)
+end
diff --git a/src/matrix_coo/matrix_coo_kernels.jl b/src/matrix_coo/matrix_coo_kernels.jl
index 1849c49..2184223 100644
--- a/src/matrix_coo/matrix_coo_kernels.jl
+++ b/src/matrix_coo/matrix_coo_kernels.jl
@@ -137,5 +137,47 @@ end
 )
     i = @index(Global)
 
-    @inbounds C[rowind[i], colind[i]] += nzval[i]
+    C[rowind[i], colind[i]] += nzval[i]
+end
+
+# Kernel for computing Kronecker product in COO format
+@kernel inbounds=true function kernel_kron_coo!(
+    @Const(rowind_A),
+    @Const(colind_A),
+    @Const(nzval_A),
+    @Const(rowind_B),
+    @Const(colind_B),
+    @Const(nzval_B),
+    rowind_C,
+    colind_C,
+    nzval_C,
+    @Const(m_B::Int),
+    @Const(n_B::Int),
+)
+    idx = @index(Global, Linear)
+
+    nnz_A = length(nzval_A)
+    nnz_B = length(nzval_B)
+
+    if idx <= nnz_A * nnz_B
+        # Compute which element from A and B we're combining
+        idx_A = div(idx - 1, nnz_B) + 1
+        idx_B = mod(idx - 1, nnz_B) + 1
+
+        # Get the row and column indices from A and B
+        i_A = rowind_A[idx_A]
+        j_A = colind_A[idx_A]
+        val_A = nzval_A[idx_A]
+
+        i_B = rowind_B[idx_B]
+        j_B = colind_B[idx_B]
+        val_B = nzval_B[idx_B]
+
+        # Compute the row and column indices in C
+        # C[i,j] = A[i_A, j_A] * B[i_B, j_B]
+        # where i = (i_A - 1) * m_B + i_B and j = (j_A - 1) * n_B + j_B
+        rowind_C[idx] = (i_A - 1) * m_B + i_B
+        colind_C[idx] = (j_A - 1) * n_B + j_B
+        nzval_C[idx] = val_A * val_B
+    end
 end
diff --git a/src/matrix_csc/matrix_csc.jl b/src/matrix_csc/matrix_csc.jl
index d9d8351..ad5280c 100644
--- a/src/matrix_csc/matrix_csc.jl
+++ b/src/matrix_csc/matrix_csc.jl
@@ -318,3 +318,36 @@ function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCSC)
 
     return C
 end
+
+"""
+    kron(A::DeviceSparseMatrixCSC, B::DeviceSparseMatrixCSC)
+
+Compute the Kronecker product of two sparse matrices in CSC format.
+
+The Kronecker product is computed by converting to COO format, computing the 
+product, and converting back to CSC format.
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A = DeviceSparseMatrixCSC(sparse([1, 2], [1, 2], [1.0, 2.0], 2, 2));
+
+julia> B = DeviceSparseMatrixCSC(sparse([1, 2], [1, 2], [3.0, 4.0], 2, 2));
+
+julia> C = kron(A, B);
+
+julia> size(C)
+(4, 4)
+
+julia> nnz(C)
+4
+```
+"""
+function LinearAlgebra.kron(A::DeviceSparseMatrixCSC, B::DeviceSparseMatrixCSC)
+    # Convert to COO, compute kron, convert back to CSC
+    A_coo = DeviceSparseMatrixCOO(A)
+    B_coo = DeviceSparseMatrixCOO(B)
+    C_coo = kron(A_coo, B_coo)
+    return DeviceSparseMatrixCSC(C_coo)
+end
diff --git a/src/matrix_csr/matrix_csr.jl b/src/matrix_csr/matrix_csr.jl
index a1fd171..e47fd1a 100644
--- a/src/matrix_csr/matrix_csr.jl
+++ b/src/matrix_csr/matrix_csr.jl
@@ -331,3 +331,40 @@ function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCSR)
 
     return C
 end
+
+"""
+    kron(A::DeviceSparseMatrixCSR, B::DeviceSparseMatrixCSR)
+
+Compute the Kronecker product of two sparse matrices in CSR format.
+
+The Kronecker product is computed by converting to COO format, computing the 
+product, and converting back to CSR format.
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A_coo = DeviceSparseMatrixCOO(sparse([1, 2], [1, 2], [1.0, 2.0], 2, 2));
+
+julia> B_coo = DeviceSparseMatrixCOO(sparse([1, 2], [1, 2], [3.0, 4.0], 2, 2));
+
+julia> A = DeviceSparseMatrixCSR(A_coo);
+
+julia> B = DeviceSparseMatrixCSR(B_coo);
+
+julia> C = kron(A, B);
+
+julia> size(C)
+(4, 4)
+
+julia> nnz(C)
+4
+```
+"""
+function LinearAlgebra.kron(A::DeviceSparseMatrixCSR, B::DeviceSparseMatrixCSR)
+    # Convert to COO, compute kron, convert back to CSR
+    A_coo = DeviceSparseMatrixCOO(A)
+    B_coo = DeviceSparseMatrixCOO(B)
+    C_coo = kron(A_coo, B_coo)
+    return DeviceSparseMatrixCSR(C_coo)
+end
diff --git a/test/shared/conversions.jl b/test/shared/conversions.jl
index 411c152..6ff77c8 100644
--- a/test/shared/conversions.jl
+++ b/test/shared/conversions.jl
@@ -10,16 +10,12 @@ function shared_test_conversions(
         # which is not supported on JLBackend. Therefore, we skip conversion
         # tests for JLArray.
         if array_type != "JLArray"
+            Tv = float_types[end]
+            Ti = int_types[end]
+            A = SparseMatrixCSC{Tv,Ti}(sprand(100, 200, 0.05))
+
             # Test CSC → COO → CSC round-trip
             @testset "CSC ↔ COO" begin
-                A = sparse(
-                    int_types[end][1, 2, 3, 1, 2],
-                    int_types[end][1, 2, 3, 2, 3],
-                    float_types[end][1.0, 2.0, 3.0, 4.0, 5.0],
-                    3,
-                    3,
-                )
-
                 # CSC → COO
                 A_csc = adapt(op, DeviceSparseMatrixCSC(A))
                 A_coo_from_csc = DeviceSparseMatrixCOO(A_csc)
@@ -37,14 +33,6 @@ function shared_test_conversions(
 
             # Test CSR → COO → CSR round-trip
             @testset "CSR ↔ COO" begin
-                A = sparse(
-                    int_types[end][1, 2, 3, 1, 2],
-                    int_types[end][1, 2, 3, 2, 3],
-                    float_types[end][1.0, 2.0, 3.0, 4.0, 5.0],
-                    3,
-                    3,
-                )
-
                 # CSR → COO
                 A_csr = adapt(op, DeviceSparseMatrixCSR(A))
                 A_coo_from_csr = DeviceSparseMatrixCOO(A_csr)
@@ -63,29 +51,11 @@ function shared_test_conversions(
             # Test with empty matrices
             @testset "Edge Cases" begin
                 # Empty matrix
-                A_empty = spzeros(float_types[end], int_types[end], 3, 3)
+                A_empty = spzeros(Tv, Ti, 3, 4)
                 A_csc_empty = adapt(op, DeviceSparseMatrixCSC(A_empty))
                 A_coo_empty = DeviceSparseMatrixCOO(A_csc_empty)
                 @test nnz(A_coo_empty) == 0
-                @test size(A_coo_empty) == (3, 3)
-            end
-
-            # Test large matrix conversion
-            @testset "Large Matrix" begin
-                A_large =
-                    SparseMatrixCSC{float_types[end],int_types[end]}(sprand(100, 100, 0.05))
-
-                # CSC → COO → CSC
-                A_csc_large = adapt(op, DeviceSparseMatrixCSC(A_large))
-                A_coo_large = DeviceSparseMatrixCOO(A_csc_large)
-                A_csc_back = DeviceSparseMatrixCSC(A_coo_large)
-                @test collect(SparseMatrixCSC(A_csc_back)) ≈ collect(A_large)
-
-                # CSR → COO → CSR
-                A_csr_large = adapt(op, DeviceSparseMatrixCSR(A_large))
-                A_coo_large2 = DeviceSparseMatrixCOO(A_csr_large)
-                A_csr_back = DeviceSparseMatrixCSR(A_coo_large2)
-                @test collect(SparseMatrixCSC(A_csr_back)) ≈ collect(A_large)
+                @test size(A_coo_empty) == (3, 4)
             end
         end
     end
diff --git a/test/shared/matrix_coo.jl b/test/shared/matrix_coo.jl
index 5d65848..861a2fc 100644
--- a/test/shared/matrix_coo.jl
+++ b/test/shared/matrix_coo.jl
@@ -277,4 +277,23 @@ function shared_test_linearalgebra_matrix_coo(
             @test_throws DimensionMismatch dA + dB_wrong
         end
     end
+
+    @testset "Kronecker Product" begin
+        for T in (int_types..., float_types..., complex_types...)
+            # Test with rectangular matrices
+            A_sparse = sprand(T, 30, 25, 0.1)
+            B_sparse = sprand(T, 20, 15, 0.1)
+
+            A = adapt(op, DeviceSparseMatrixCOO(A_sparse))
+            B = adapt(op, DeviceSparseMatrixCOO(B_sparse))
+
+            C = kron(A, B)
+            C_expected = kron(A_sparse, B_sparse)
+
+            @test size(C) == size(C_expected)
+            @test nnz(C) == nnz(C_expected)
+            @test Matrix(SparseMatrixCSC(C)) ≈ Matrix(C_expected)
+            @test C isa DeviceSparseMatrixCOO
+        end
+    end
 end
diff --git a/test/shared/matrix_csc.jl b/test/shared/matrix_csc.jl
index f46b702..49490db 100644
--- a/test/shared/matrix_csc.jl
+++ b/test/shared/matrix_csc.jl
@@ -275,4 +275,25 @@ function shared_test_linearalgebra_matrix_csc(
             @test_throws DimensionMismatch dA + dB_wrong
         end
     end
+
+    @testset "Kronecker Product" begin
+        if array_type != "JLArray"
+            for T in (int_types..., float_types..., complex_types...)
+                # Test with rectangular matrices
+                A_sparse = sprand(T, 30, 25, 0.1)
+                B_sparse = sprand(T, 20, 15, 0.1)
+
+                A = adapt(op, DeviceSparseMatrixCSC(A_sparse))
+                B = adapt(op, DeviceSparseMatrixCSC(B_sparse))
+
+                C = kron(A, B)
+                C_expected = kron(A_sparse, B_sparse)
+
+                @test size(C) == size(C_expected)
+                @test nnz(C) == nnz(C_expected)
+                @test Matrix(SparseMatrixCSC(C)) ≈ Matrix(C_expected)
+                @test C isa DeviceSparseMatrixCSC
+            end
+        end
+    end
 end
diff --git a/test/shared/matrix_csr.jl b/test/shared/matrix_csr.jl
index 550d93d..7f93b08 100644
--- a/test/shared/matrix_csr.jl
+++ b/test/shared/matrix_csr.jl
@@ -274,4 +274,25 @@ function shared_test_linearalgebra_matrix_csr(
             @test_throws DimensionMismatch dA + dB_wrong
         end
     end
+
+    @testset "Kronecker Product" begin
+        if array_type != "JLArray"
+            for T in (int_types..., float_types..., complex_types...)
+                # Test with rectangular matrices
+                A_sparse = sprand(T, 30, 25, 0.1)
+                B_sparse = sprand(T, 20, 15, 0.1)
+
+                A = adapt(op, DeviceSparseMatrixCSR(A_sparse))
+                B = adapt(op, DeviceSparseMatrixCSR(B_sparse))
+
+                C = kron(A, B)
+                C_expected = kron(A_sparse, B_sparse)
+
+                @test size(C) == size(C_expected)
+                @test nnz(C) == nnz(C_expected)
+                @test Matrix(SparseMatrixCSC(C)) ≈ Matrix(C_expected)
+                @test C isa DeviceSparseMatrixCSR
+            end
+        end
+    end
 end