diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp
index f251d632117..42b497a9923 100644
--- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp
+++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash_impl.hpp
@@ -358,7 +358,8 @@ void NegativeTwoPhaseFlash::computeDerivatives(
   }
 
   // Solve linear system
-  solveLinearSystem( A.toSlice(), X.toSlice() );
+  bool const solveStatus = solveLinearSystem( A.toSlice(), X.toSlice() );
+  GEOS_ERROR_IF( !solveStatus, "Failed to solve for the derivatives." );
 
   // Fill in the derivatives
   for( integer i = 0; i < numComps; ++i )
diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp
index 3e017c005ea..2b9f3e21f38 100644
--- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp
+++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/StabilityTest.hpp
@@ -325,7 +325,8 @@ struct StabilityTest
     }
 
     // Solve linear system
-    solveLinearSystem( A.toSlice(), X.toSlice() );
+    bool const solveStatus = solveLinearSystem( A.toSlice(), X.toSlice() );
+    GEOS_ERROR_IF( !solveStatus, "Failed to solve for the derivatives." );
 
     for( integer idof = 0; idof < numDofs; ++idof )
     {
diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp
index 666c57a116e..a9916863563 100644
--- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp
+++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/Utilities.hpp
@@ -21,7 +21,7 @@
 #define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_FUNCTIONS_UTILITIES_HPP_
 
 #include "constitutive/fluid/multifluid/MultiFluidConstants.hpp"
-#include "denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp"
+#include "denseLinearAlgebra/denseLASolvers.hpp"
 
 namespace geos
 {
@@ -166,14 +166,21 @@ GEOS_HOST_DEVICE
 static bool solveLinearSystem( arraySlice2d< real64, USD > const & A,
                                arraySlice2d< real64, USD > const & X )
 {
-#if defined(GEOS_DEVICE_COMPILE)
-  GEOS_UNUSED_VAR( A );
-  GEOS_UNUSED_VAR( X );
-  return false;
-#else
-  BlasLapackLA::solveLinearSystem( A, X );
-  return true;
-#endif
+  // The flash and stability-test systems are A: N x N and X: N x (N+1).
+  // denseLinearAlgebra::solve is GEOS_HOST_DEVICE but needs both dimensions at compile time.
+  switch( A.size( 0 ) )
+  {
+    case 2:  return denseLinearAlgebra::solve< 2, 3 >( A, X );
+    case 3:  return denseLinearAlgebra::solve< 3, 4 >( A, X );
+    case 4:  return denseLinearAlgebra::solve< 4, 5 >( A, X );
+    case 5:  return denseLinearAlgebra::solve< 5, 6 >( A, X );
+    case 6:  return denseLinearAlgebra::solve< 6, 7 >( A, X );
+    case 7:  return denseLinearAlgebra::solve< 7, 8 >( A, X );
+    case 8:  return denseLinearAlgebra::solve< 8, 9 >( A, X );
+    case 9:  return denseLinearAlgebra::solve< 9, 10 >( A, X );
+    case 10:  return denseLinearAlgebra::solve< 10, 11 >( A, X );
+    default: return false;
+  }
 }
 
 } // namespace compositional
diff --git a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp
index f0e5e6fbed2..afcbca96e1c 100644
--- a/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp
+++ b/src/coreComponents/denseLinearAlgebra/denseLASolvers.hpp
@@ -519,7 +519,7 @@ inline
 bool solve( MATRIX_TYPE & A, RHS_TYPE & b, SOL_TYPE && x )
 {
   static_assert( N > 0, "N must be greater than 0." );
-  static_assert( N < 10, "N cannot be larger than 9" );
+  static_assert( N < 11, "N cannot be larger than 10" );
   LvArray::tensorOps::internal::checkSizes< N, N >( A );
   LvArray::tensorOps::internal::checkSizes< N >( b );
   LvArray::tensorOps::internal::checkSizes< N >( x );
@@ -587,9 +587,9 @@ inline
 bool solve( MATRIX_TYPE & A, SOL_TYPE && X )
 {
   static_assert( N > 0, "N must be greater than 0." );
-  static_assert( N < 10, "N cannot be larger than 9" );
+  static_assert( N < 11, "N cannot be larger than 10" );
   static_assert( M > 0, "M must be greater than 0." );
-  static_assert( M < 10, "M cannot be larger than 9" );
+  static_assert( M < 12, "M cannot be larger than 11" );
   LvArray::tensorOps::internal::checkSizes< N, N >( A );
   LvArray::tensorOps::internal::checkSizes< N, M >( X );
 
@@ -656,9 +656,9 @@ inline
 bool solve( MATRIX_TYPE & A, RHS_TYPE & B, SOL_TYPE && X )
 {
   static_assert( N > 0, "N must be greater than 0." );
-  static_assert( N < 10, "N cannot be larger than 9" );
+  static_assert( N < 11, "N cannot be larger than 10" );
   static_assert( M > 0, "M must be greater than 0." );
-  static_assert( M < 10, "M cannot be larger than 9" );
+  static_assert( M < 12, "M cannot be larger than 11" );
   LvArray::tensorOps::internal::checkSizes< N, N >( A );
   LvArray::tensorOps::internal::checkSizes< N, M >( B );
   LvArray::tensorOps::internal::checkSizes< N, M >( X );