diff --git a/src/ParameterTypes/passive_ren.jl b/src/ParameterTypes/passive_ren.jl
index 577eaf4..51503ef 100644
--- a/src/ParameterTypes/passive_ren.jl
+++ b/src/ParameterTypes/passive_ren.jl
@@ -25,7 +25,7 @@ Construct direct parameterisation of a passive REN.
 - `ν::Number=0`: Passivity index. Use `ν > 0` for an incrementally strictly input passive model. Set both `ν = 0` and `ρ = 0` for incrementally passive model.
 - `ρ::Number=0`: Passivity index. Use `ρ > 0` for an incrementally strictly output passive model. 
 
-Note that setting both `ν,ρ > 0` or both `ν,ρ < 0` is not currently supported and will throw an error.
+Note that the product of passivity indices ρν has to be less than 1/4 for passive REN.
 
 # Keyword arguments
 
@@ -55,8 +55,8 @@ function PassiveRENParams{T}(
     end
 
     # Check ρ and ν
-    if ρ*ν > 0
-        error("If ρ and ν are both positive, passiveREN could produce incorrect results. Please set at least one of them as zero. ")               
+    if ρ*ν >= 1/4
+        error("ρ and ν can not be arbitrarily large for passiveREN models. Please make sure ρν < 1/4. ")               
     end
 
     if ρ < 0 || ν < 0
@@ -114,7 +114,8 @@ function direct_to_explicit(ps::PassiveRENParams{T}, return_h=false) where T
         H = x_to_h(X, ϵ, polar_param, ρ_polar) + Γ2
     else    
         # For ρ!=0 case, ISOP model
-        D22 = _D22_pass(M, ρ) 
+        D22 = _D22_pass(M, ρ)
+
         C2_imp = _C2_pass(D22, C2, ρ)
         D21_imp = _D21_pass(D22, D21, D12_imp, ρ)
 
@@ -132,8 +133,6 @@ function direct_to_explicit(ps::PassiveRENParams{T}, return_h=false) where T
 
 end
 
-_D22_pass(M, ρ) = ((I+M) \ I) / ρ   
-
 _C2_pass(D22, C2, ρ) = (D22'*(-2ρ*I) + I)*C2
 
 _D21_pass(D22, D21, D12_imp, ρ) = (D22'*(-2ρ*I) + I)*D21 - D12_imp'
@@ -142,6 +141,11 @@ _M_pass(X3, Y3, ϵ) = X3'*X3 + Y3 - Y3' + ϵ*I
 
 _R_pass(D22, ν, ρ) = -2ν*I + D22 + D22' + D22'*(-2ρ*I)*D22
 
+function _D22_pass(M, ρ)
+    Im = _I(M) # Prevents scalar indexing on backwards pass of () / (I + M) on GPU
+    return ((Im + M) \ Im) / ρ  
+end
+
 function _Γ1_pass(nx, ny, C2, D21, ρ, T) 
     [C2'; D21'; zeros(T, nx, ny)] * (-2ρ*I) * [C2 D21 zeros(T, ny, nx)]
 end