diff --git a/src/mprk.jl b/src/mprk.jl
index 6e54e4b2..5b2dbb6f 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -351,7 +351,9 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
         # as well as to store the system matrix of the linear system
         # Right hand side of linear system is always uprev
         linprob = LinearProblem(P, _vec(uprev))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         MPEConservativeCache(P, σ, tab, linsolve)
@@ -360,7 +362,9 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
         # We use P to store the evaluation of the PDS
         # as well as to store the system matrix of the linear system
         linprob = LinearProblem(P, _vec(linsolve_rhs))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         MPECache(P, similar(u), σ, tab, linsolve_rhs, linsolve)
@@ -658,7 +662,9 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
         # not be altered, since it is needed to compute the adaptive time step
         # size.
         linprob = LinearProblem(P2, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P2),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         MPRK22ConservativeCache(tmp, P, P2, σ,
@@ -666,7 +672,9 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                                 linsolve)
     elseif f isa PDSFunction
         linprob = LinearProblem(P2, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P2),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         MPRK22Cache(tmp, P, P2,
@@ -1242,7 +1250,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt
         # not be altered, since it is needed to compute the adaptive time step
         # size.
         linprob = LinearProblem(P3, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P3),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
         MPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, tab, linsolve)
     elseif f isa PDSFunction
@@ -1251,7 +1261,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt
         D3 = similar(u)
 
         linprob = LinearProblem(P3, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P3),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         MPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, tab, linsolve)
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index c25c79a3..16f0f443 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -240,7 +240,9 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
 
     if f isa ConservativePDSFunction
         linprob = LinearProblem(P2, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P2),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         SSPMPRK22ConservativeCache(tmp, P, P2, σ,
@@ -248,7 +250,9 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                                    linsolve)
     elseif f isa PDSFunction
         linprob = LinearProblem(P2, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P2),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         SSPMPRK22Cache(tmp, P, P2,
@@ -756,7 +760,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
 
     if f isa ConservativePDSFunction
         linprob = LinearProblem(P3, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P3),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
         SSPMPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, ρ, tab, linsolve)
     elseif f isa PDSFunction
@@ -765,7 +771,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
         D3 = similar(u)
 
         linprob = LinearProblem(P3, _vec(tmp))
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve;
+                        alias_A = !issparse(P3),
+                        alias_b = true,
                         assumptions = LinearSolve.OperatorAssumptions(true))
 
         SSPMPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, tab, linsolve)