Scalar indexing is disallowed is Cuda.jl so use allowscalar

Author: Fe-r-oz

src/steadystate.jl

@@ -369,8 +369,7 @@ function _steadystate_fourier(
     n_fourier = 2 * n_max + 1
     n_list = (-n_max):n_max
 
-    # Create arrays with proper GPU/CPU backend
-    weight = one(T)
+    weight = 1
     rows = _dense_similar(L_0_mat, Ns)
     cols = _dense_similar(L_0_mat, Ns)
     vals = _dense_similar(L_0_mat, Ns)
@@ -382,7 +381,6 @@ function _steadystate_fourier(
     Mn = _sparse_similar(L_0_mat, rows, cols, vals, N, N)
     L = L_0_mat + Mn
 
-    # Initialize the big matrix M with proper backend
     M = _sparse_similar(L_0_mat, n_fourier * N, n_fourier * N)
 
     # Add superdiagonal blocks (L_m)
@@ -391,7 +389,7 @@ function _steadystate_fourier(
         cols_block = _dense_similar(L_0_mat, N)
         fill!(rows_block, i)
         fill!(cols_block, i+1)
-        block = _sparse_similar(L_0_mat, rows_block, cols_block, ones(T, N), n_fourier, n_fourier)
+        block = _sparse_similar(L_0_mat, rows_block, cols_block, ones(N), n_fourier, n_fourier)
         M += kron(block, L_m_mat)
     end
 
@@ -401,7 +399,7 @@ function _steadystate_fourier(
         cols_block = _dense_similar(L_0_mat, N)
         fill!(rows_block, i+1)
         fill!(cols_block, i)
-        block = _sparse_similar(L_0_mat, rows_block, cols_block, ones(T, N), n_fourier, n_fourier)
+        block = _sparse_similar(L_0_mat, rows_block, cols_block, ones(N), n_fourier, n_fourier)
         M += kron(block, L_p_mat)
     end
 
@@ -413,16 +411,14 @@ function _steadystate_fourier(
         cols_block = _dense_similar(L_0_mat, N)
         fill!(rows_block, i)
         fill!(cols_block, i)
-        block = _sparse_similar(L_0_mat, rows_block, cols_block, ones(T, N), n_fourier, n_fourier)
+        block = _sparse_similar(L_0_mat, rows_block, cols_block, ones(N), n_fourier, n_fourier)
         M += kron(block, block_diag)
     end
 
-    # Initialize solution vector with proper backend
     v0 = _dense_similar(L_0_mat, n_fourier * N)
-    fill!(v0, zero(T))
-    allowed_setindex!(v0, weight, n_max * N + 1)
+    fill!(v0, 0)
+    @allowscalar v0[n_max * N + 1] = weight
 
-    # Prepare preconditioners if needed
     if !isnothing(solver.Pl)
         kwargs = merge((; kwargs...), (Pl = solver.Pl(M),))
     elseif isa(M, SparseMatrixCSC)
@@ -431,11 +427,9 @@ function _steadystate_fourier(
     !isnothing(solver.Pr) && (kwargs = merge((; kwargs...), (Pr = solver.Pr(M),)))
     kwargs = merge((abstol = tol, reltol = tol), kwargs)
 
-    # Solve the linear system
     prob = LinearProblem(M, v0)
     ρtot = solve(prob, solver.alg; kwargs...).u
 
-    # Extract results
     offset1 = n_max * N
     offset2 = (n_max + 1) * N
     ρ0 = reshape(ρtot[(offset1+1):offset2], Ns, Ns)