diff --git a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.cpp b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.cpp
index d05bbf8aef2..541a802861b 100644
--- a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.cpp
+++ b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.cpp
@@ -379,6 +379,151 @@ void matrixLeastSquaresSolutionSolve( arraySlice2d< real64, USD > const & A,
   GEOS_ERROR_IF( INFO != 0, "The algorithm computing matrix linear system failed to converge." );
 }
 
+template< int USD1, int USD2 >
+GEOS_FORCE_INLINE
+void matrixCopy( int const N,
+                 int const M,
+                 arraySlice2d< real64, USD1 > const & A,
+                 arraySlice2d< real64, USD2 > const & B )
+{
+  for( int i = 0; i < N; i++ )
+  {
+    for( int j = 0; j < M; j++ )
+    {
+      B( i, j ) = A( i, j );
+    }
+  }
+}
+
+template< int USD >
+GEOS_FORCE_INLINE
+void matrixTranspose( int const N,
+                      arraySlice2d< real64, USD > const & A )
+{
+  for( int i = 0; i < N; i++ )
+  {
+    for( int j = i+1; j < N; j++ )
+    {
+      std::swap( A( i, j ), A( j, i ) );
+    }
+  }
+}
+
+template< typename T, int USD >
+void solveLinearSystem( arraySlice2d< T, USD > const & A,
+                        arraySlice2d< real64 const, USD > const & B,
+                        arraySlice2d< real64, USD > const & X )
+{
+  // --- Check that source matrix is square
+  int const N = LvArray::integerConversion< int >( A.size( 0 ) );
+  GEOS_ASSERT_MSG( N > 0 &&
+                   N == A.size( 1 ),
+                   "Matrix must be square" );
+
+  // --- Check that rhs B has appropriate dimensions
+  GEOS_ASSERT_MSG( B.size( 0 ) == N,
+                   "right-hand-side matrix has wrong dimensions" );
+  int const M = LvArray::integerConversion< int >( B.size( 1 ) );
+
+  // --- Check that solution X has appropriate dimensions
+  GEOS_ASSERT_MSG( X.size( 0 ) == N &&
+                   X.size( 1 ) == M,
+                   "solution matrix has wrong dimensions" );
+
+  // --- Check that everything is contiguous
+  GEOS_ASSERT_MSG( A.isContiguous(), "Matrix is not contiguous" );
+  GEOS_ASSERT_MSG( B.isContiguous(), "right-hand-side matrix is not contiguous" );
+  GEOS_ASSERT_MSG( X.isContiguous(), "solution matrix is not contiguous" );
+
+  real64 * matrixData = nullptr;
+  array2d< real64 > LU;   // Space for LU-factors
+  if constexpr ( !std::is_const< T >::value )
+  {
+    matrixData = A.dataIfContiguous();
+  }
+  else
+  {
+    LU.resize( N, N );
+    matrixData = LU.data();
+    // Direct copy here ignoring permutation
+    int const INCX = 1;
+    int const INCY = 1;
+    int const K = LvArray::integerConversion< int >( A.size( ) );
+    GEOS_dcopy( &K, A.dataIfContiguous(), &INCX, matrixData, &INCY );
+  }
+
+  array1d< int > IPIV( N );
+  int INFO;
+  char const TRANS = (USD == MatrixLayout::ROW_MAJOR) ? 'T' : 'N';
+
+  GEOS_dgetrf( &N, &N, matrixData, &N, IPIV.data(), &INFO );
+
+  GEOS_ASSERT_MSG( INFO == 0, "LAPACK dgetrf error code: " << INFO );
+
+  if constexpr ( std::is_const< T >::value )
+  {
+    int const INCX = 1;
+    int const INCY = 1;
+    int const K = LvArray::integerConversion< int >( B.size( ) );
+    GEOS_dcopy( &K, B.dataIfContiguous(), &INCX, X.dataIfContiguous(), &INCY );
+  }
+
+  // For row-major form, we need to reorder into col-major form
+  // This might require an extra allocation
+  real64 * solutionData = X.dataIfContiguous();
+  array2d< real64, MatrixLayout::COL_MAJOR_PERM > X0;
+  if constexpr ( USD == MatrixLayout::ROW_MAJOR )
+  {
+    if( 1 < M && M == N )
+    {
+      // Square case: swap in place
+      matrixTranspose( N, X );
+    }
+    else if( 1 < M )
+    {
+      X0.resize( N, M );
+      matrixCopy( N, M, X, X0.toSlice() );
+      solutionData = X0.data();
+    }
+  }
+
+  GEOS_dgetrs( &TRANS, &N, &M, matrixData, &N, IPIV.data(), solutionData, &N, &INFO );
+
+  GEOS_ASSERT_MSG( INFO == 0, "LAPACK dgetrs error code: " << INFO );
+
+  if constexpr ( USD == MatrixLayout::ROW_MAJOR )
+  {
+    if( 1 < M && M == N )
+    {
+      // Square case: swap in place
+      matrixTranspose( N, X );
+    }
+    else if( 1 < M )
+    {
+      matrixCopy( N, M, X0.toSlice(), X );
+    }
+  }
+}
+
+template< typename T, int USD >
+void solveLinearSystem( arraySlice2d< T, USD > const & A,
+                        arraySlice1d< real64 const > const & b,
+                        arraySlice1d< real64 > const & x )
+{
+  // --- Check that b and x have the same size
+  int const N = LvArray::integerConversion< int >( b.size( 0 ) );
+  GEOS_ASSERT_MSG( 0 < N && x.size() == N,
+                   "right-hand-side and/or solution has wrong dimensions" );
+
+  // Create 2d slices
+  int const dims[2] = {N, 1};
+  int const strides[2] = {1, 1};
+  arraySlice2d< real64 const, USD > B( b.dataIfContiguous(), dims, strides );
+  arraySlice2d< real64, USD > X( x.dataIfContiguous(), dims, strides );
+
+  solveLinearSystem( A, B, X );
+}
+
 } // namespace detail
 
 real64 BlasLapackLA::determinant( arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & A )
@@ -885,60 +1030,56 @@ void BlasLapackLA::matrixEigenvalues( MatRowMajor< real64 const > const & A,
   matrixEigenvalues( AT.toSliceConst(), lambda );
 }
 
-void BlasLapackLA::solveLinearSystem( MatColMajor< real64 const > const & A,
-                                      arraySlice1d< real64 const > const & rhs,
-                                      arraySlice1d< real64 > const & solution )
+void BlasLapackLA::solveLinearSystem( MatRowMajor< real64 const > const & A,
+                                      Vec< real64 const > const & rhs,
+                                      Vec< real64 > const & solution )
 {
-  // --- Check that source matrix is square
-  int const NN = LvArray::integerConversion< int >( A.size( 0 ));
-  GEOS_ASSERT_MSG( NN > 0 &&
-                   NN == A.size( 1 ),
-                   "Matrix must be square" );
-
-  // --- Check that rhs and solution have appropriate dimension
-  GEOS_ASSERT_MSG( rhs.size( 0 ) == NN,
-                   "right-hand-side vector has wrong dimensions" );
-
-  GEOS_ASSERT_MSG( solution.size( 0 ) == NN,
-                   "solution vector has wrong dimensions" );
-
-  array1d< int > IPIV;
-  IPIV.resize( NN );
-  int const NRHS = 1; // we only allow for 1 rhs vector.
-  int INFO;
-
-  // make a copy of A, since dgeev destroys contents
-  array2d< real64, MatrixLayout::COL_MAJOR_PERM > ACOPY( A.size( 0 ), A.size( 1 ) );
-  BlasLapackLA::matrixCopy( A, ACOPY );
-
-  // copy the rhs in the solution vector
-  BlasLapackLA::vectorCopy( rhs, solution );
-
-  GEOS_dgetrf( &NN, &NN, ACOPY.data(), &NN, IPIV.data(), &INFO );
+  detail::solveLinearSystem( A, rhs, solution );
+}
 
-  GEOS_ASSERT_MSG( INFO == 0, "LAPACK dgetrf error code: " << INFO );
+void BlasLapackLA::solveLinearSystem( MatColMajor< real64 const > const & A,
+                                      Vec< real64 const > const & rhs,
+                                      Vec< real64 > const & solution )
+{
+  detail::solveLinearSystem( A, rhs, solution );
+}
 
-  GEOS_dgetrs( "N", &NN, &NRHS, ACOPY.data(), &NN, IPIV.data(), solution.dataIfContiguous(), &NN, &INFO );
+void BlasLapackLA::solveLinearSystem( MatRowMajor< real64 > const & A,
+                                      Vec< real64 > const & rhs )
+{
+  detail::solveLinearSystem( A, rhs.toSliceConst(), rhs );
+}
 
-  GEOS_ASSERT_MSG( INFO == 0, "LAPACK dgetrs error code: " << INFO );
+void BlasLapackLA::solveLinearSystem( MatColMajor< real64 > const & A,
+                                      Vec< real64 > const & rhs )
+{
+  detail::solveLinearSystem( A, rhs.toSliceConst(), rhs );
 }
 
 void BlasLapackLA::solveLinearSystem( MatRowMajor< real64 const > const & A,
-                                      arraySlice1d< real64 const > const & rhs,
-                                      arraySlice1d< real64 > const & solution )
+                                      MatRowMajor< real64 const > const & rhs,
+                                      MatRowMajor< real64 > const & solution )
 {
-  array2d< real64, MatrixLayout::COL_MAJOR_PERM > AT( A.size( 0 ), A.size( 1 ) );
+  detail::solveLinearSystem( A, rhs, solution );
+}
 
-  // convert A to a column major format
-  for( int i = 0; i < A.size( 0 ); ++i )
-  {
-    for( int j = 0; j < A.size( 1 ); ++j )
-    {
-      AT( i, j ) = A( i, j );
-    }
-  }
+void BlasLapackLA::solveLinearSystem( MatColMajor< real64 const > const & A,
+                                      MatColMajor< real64 const > const & rhs,
+                                      MatColMajor< real64 > const & solution )
+{
+  detail::solveLinearSystem( A, rhs, solution );
+}
+
+void BlasLapackLA::solveLinearSystem( MatRowMajor< real64 > const & A,
+                                      MatRowMajor< real64 > const & rhs )
+{
+  detail::solveLinearSystem( A, rhs.toSliceConst(), rhs );
+}
 
-  solveLinearSystem( AT.toSliceConst(), rhs, solution );
+void BlasLapackLA::solveLinearSystem( MatColMajor< real64 > const & A,
+                                      MatColMajor< real64 > const & rhs )
+{
+  detail::solveLinearSystem( A, rhs.toSliceConst(), rhs );
 }
 
 void BlasLapackLA::matrixLeastSquaresSolutionSolve( arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & A,
diff --git a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp
index adc72007e49..488ed59a9b9 100644
--- a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp
+++ b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp
@@ -30,6 +30,7 @@ namespace geos
  * \class BlasLapackLA
  * \brief This class contains a collection of BLAS and LAPACK linear
  *        algebra operations for GEOSX array1d and array2d
+ * \warning These methods are currently not supported on GPUs
  */
 struct BlasLapackLA
 {
@@ -454,18 +455,88 @@ struct BlasLapackLA
    * @param [in]  rhs GEOSX array1d.
    * @param [out] solution GEOSX array1d.
    */
+  static void solveLinearSystem( MatColMajor< real64 const > const & A,
+                                 Vec< real64 const > const & rhs,
+                                 Vec< real64 > const & solution );
+
+  /**
+   * @copydoc solveLinearSystem( MatColMajor<real64 const> const &, Vec< real64 const > const &, Vec< real64 const > const & )
+   */
   static void solveLinearSystem( MatRowMajor< real64 const > const & A,
                                  Vec< real64 const > const & rhs,
                                  Vec< real64 > const & solution );
 
   /**
-   * @copydoc solveLinearSystem( MatRowMajor<real64 const> const &, Vec< real64 const > const &, Vec< real64 const > const & )
+   * @brief Solves the linear system ;
+   * \p A  \p solution = \p rhs.
+   *
+   * @details The method is intended for the solution of a small dense linear system.
+   * This solves the system in-place without allocating extra memory for the matrix or the solution. This means
+   * that at on exit the matrix is modified replaced by the LU factors and the right hand side vector is
+   * replaced by the solution.
+   * It employs lapack method dgetr.
+   *
+   * @param [in/out]  A GEOSX array2d. The matrix. On exit this will be replaced by the factorisation of A
+   * @param [in/out]  rhs GEOSX array1d. The right hand side. On exit this will be the solution
+   */
+  static void solveLinearSystem( MatColMajor< real64 > const & A,
+                                 Vec< real64 > const & rhs );
+
+  /**
+   * @copydoc solveLinearSystem( MatColMajor< real64 > const &, Vec< real64 > const & )
+   */
+  static void solveLinearSystem( MatRowMajor< real64 > const & A,
+                                 Vec< real64 > const & rhs );
+
+  /**
+   * @brief Solves the linear system ;
+   * \p A  \p solution = \p rhs.
+   *
+   * @details The method is intended for the solution of a small dense linear system in which A is an NxN matrix, the
+   * right-hand-side and the solution are matrices of size NxM.
+   * It employs lapack method dgetr.
    *
-   * @note this function first applies a matrix permutation and then calls the row major version of the function.
+   * @param [in]  A GEOSX array2d.
+   * @param [in]  rhs GEOSX array2d.
+   * @param [out] solution GEOSX array2d.
    */
   static void solveLinearSystem( MatColMajor< real64 const > const & A,
-                                 Vec< real64 const > const & rhs,
-                                 Vec< real64 > const & solution );
+                                 MatColMajor< real64 const > const & rhs,
+                                 MatColMajor< real64 > const & solution );
+
+  /**
+   * @copydoc solveLinearSystem( MatColMajor< real64 const > const &, MatColMajor< real64 const > const &, MatColMajor< const > const & )
+   *
+   * @note this function will allocate space to reorder the solution into column major form.
+   */
+  static void solveLinearSystem( MatRowMajor< real64 const > const & A,
+                                 MatRowMajor< real64 const > const & rhs,
+                                 MatRowMajor< real64 > const & solution );
+
+  /**
+   * @brief Solves the linear system ;
+   * \p A  \p solution = \p rhs.
+   *
+   * @details The method is intended for the solution of a small dense linear system in which A is an NxN matrix, the
+   * right-hand-side and the solution are matrices of size NxM.
+   * This solves the system in-place without allocating extra memory for the matrix or the solution. This means
+   * that at on exit the matrix is modified replaced by the LU factors and the right hand side vector is
+   * replaced by the solution.
+   * It employs lapack method dgetr.
+   *
+   * @param [in/out]  A GEOSX array2d. The matrix. On exit this will be replaced by the factorisation of A
+   * @param [in/out]  rhs GEOSX array1d. The right hand side. On exit this will be the solution
+   */
+  static void solveLinearSystem( MatColMajor< real64 > const & A,
+                                 MatColMajor< real64 > const & rhs );
+
+  /**
+   * @copydoc solveLinearSystem( MatColMajor< real64 > const &, MatRowMajor< real64 > const & )
+   *
+   * @note this function will allocate space to reorder the solution into column major form.
+   */
+  static void solveLinearSystem( MatRowMajor< real64 > const & A,
+                                 MatRowMajor< real64 > const & rhs );
 
   /**
    * @brief Vector copy;
diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt b/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt
index 1b7d0eefb3d..114dee90aac 100644
--- a/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt
+++ b/src/coreComponents/denseLinearAlgebra/unitTests/CMakeLists.txt
@@ -1,5 +1,6 @@
 set( serial_tests
-     testBlasLapack.cpp )
+     testBlasLapack.cpp
+     testSolveLinearSystem.cpp )
 
 set( dependencyList gtest denseLinearAlgebra )
 
diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testSolveLinearSystem.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testSolveLinearSystem.cpp
new file mode 100644
index 00000000000..48cf7b65dc8
--- /dev/null
+++ b/src/coreComponents/denseLinearAlgebra/unitTests/testSolveLinearSystem.cpp
@@ -0,0 +1,455 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019-     GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file testSolveLinearSystem.cpp
+ */
+
+// Source includes
+#include "denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp"
+
+#include "gtest/gtest.h"
+
+using namespace geos;
+
+constexpr int MAX_SIZE = 20;
+constexpr real64 machinePrecision = 1.0e2 * LvArray::NumericLimits< real64 >::epsilon;
+
+// Test matrices
+enum TestMatrixType
+{
+  LAPLACE,
+  SAMPLING,
+  GRCAR
+};
+
+template< TestMatrixType MATRIX_TYPE >
+struct TestMatrix {};
+
+template<>
+struct TestMatrix< TestMatrixType::LAPLACE >
+{
+  template< int USD >
+  static void create( arraySlice2d< real64, USD > const & A )
+  {
+    int const N = LvArray::integerConversion< int >( A.size( 0 ) );
+    GEOS_ASSERT( 2 <= N );
+    LvArray::forValuesInSlice( A, []( real64 & a ){ a = 0.0; } );
+    for( int i = 0; i < N-1; i++ )
+    {
+      A( i+1, i ) = -1.0;
+      A( i+1, i+1 ) = 2.0;
+      A( i, i+1 ) = -1.0;
+    }
+    A( 0, 0 ) = 2.0;
+  }
+};
+
+template<>
+struct TestMatrix< TestMatrixType::SAMPLING >
+{
+  template< int USD >
+  static void create( arraySlice2d< real64, USD > const & A )
+  {
+    int const N = LvArray::integerConversion< int >( A.size( 0 ) );
+    GEOS_ASSERT( 2 <= N );
+
+    real64 minajj = LvArray::NumericLimits< real64 >::max;
+    for( int j = 1; j <= N; j++ )
+    {
+      real64 ajj = 0.0;
+      for( int i = 1; i <= N; i++ )
+      {
+        if( i != j )
+        {
+          A( i-1, j-1 ) = static_cast< real64 >(i)/static_cast< real64 >(i-j);
+          ajj += A( i-1, j-1 );
+        }
+      }
+      A( j-1, j-1 ) = ajj;
+      minajj = LvArray::math::min( ajj, minajj );
+    }
+    for( int j = 0; j < N; j++ )
+    {
+      A( j, j ) -= minajj;
+    }
+  }
+};
+
+template<>
+struct TestMatrix< TestMatrixType::GRCAR >
+{
+  template< int USD >
+  static void create( arraySlice2d< real64, USD > const & A )
+  {
+    int const N = LvArray::integerConversion< int >( A.size( 0 ) );
+    GEOS_ASSERT( 4 < N );
+    LvArray::forValuesInSlice( A, []( real64 & a ){ a = 0.0; } );
+    for( int i = 0; i < N-1; i++ )
+    {
+      A( i+1, i ) = -1.0;
+    }
+    for( int c = 0; c < 4; c++ )
+    {
+      for( int i = 0; i < N-c; i++ )
+      {
+        A( i, i+c ) = 1.0;
+      }
+    }
+  }
+};
+
+// Randomly reorder a matrix
+template< int USD >
+void random_permutation( arraySlice2d< real64, USD > const & A )
+{
+  int const N = LvArray::integerConversion< int >( A.size( 0 ) );
+  GEOS_ASSERT( 2 <= N );
+
+  StackArray< real64, 1, MAX_SIZE > ordering( N );
+  BlasLapackLA::vectorRand( ordering.toSlice() );
+
+  for( int j = 0; j < N; j++ )
+  {
+    int const i = static_cast< int >(ordering[j]*N);
+    if( i != j )
+    {
+      for( int k = 0; k < N; ++k )
+      {
+        std::swap( A( i, k ), A( j, k ));
+      }
+    }
+  }
+}
+
+// Local naive matrix-vector multiply
+template< int USD >
+void matrix_vector_multiply( arraySlice2d< real64 const, USD > const & A,
+                             arraySlice1d< real64 const > const & x,
+                             arraySlice1d< real64 > const & b )
+{
+  int const N = LvArray::integerConversion< int >( A.size( 0 ) );
+  for( int i = 0; i < N; ++i )
+  {
+    real64 bi = 0.0;
+    for( int j = 0; j < N; ++j )
+    {
+      bi += A( i, j )*x( j );
+    }
+    b( i ) = bi;
+  }
+}
+
+// Local naive matrix-matrix multiply
+template< int USD >
+void matrix_matrix_multiply( arraySlice2d< real64 const, USD > const & A,
+                             arraySlice2d< real64 const, USD > const & X,
+                             arraySlice2d< real64, USD > const & B )
+{
+  int const K = LvArray::integerConversion< int >( A.size( 1 ) );
+  int const M = LvArray::integerConversion< int >( X.size( 0 ) );
+  int const N = LvArray::integerConversion< int >( X.size( 1 ) );
+  for( int i = 0; i < M; ++i )
+  {
+    for( int j = 0; j < N; ++j )
+    {
+      real64 bij = 0.0;
+      for( int k = 0; k < K; ++k )
+      {
+        bij += A( i, k )*X( k, j );
+      }
+      B( i, j ) = bij;
+    }
+  }
+}
+
+template< typename MatrixType >
+struct ArrayType {};
+
+template<>
+struct ArrayType< Array< real64, 2, MatrixLayout::COL_MAJOR_PERM > >
+{
+  using type = Array< real64, 1 >;
+};
+template<>
+struct ArrayType< Array< real64, 2, MatrixLayout::ROW_MAJOR_PERM > >
+{
+  using type = Array< real64, 1 >;
+};
+template<>
+struct ArrayType< StackArray< real64, 2, MAX_SIZE *MAX_SIZE, MatrixLayout::COL_MAJOR_PERM > >
+{
+  using type = StackArray< real64, 1, MAX_SIZE >;
+};
+template<>
+struct ArrayType< StackArray< real64, 2, MAX_SIZE *MAX_SIZE, MatrixLayout::ROW_MAJOR_PERM > >
+{
+  using type = StackArray< real64, 1, MAX_SIZE >;
+};
+
+template< typename MatrixType >
+class LinearSolveFixture : public ::testing::Test
+{
+public:
+  using VectorType = typename ArrayType< MatrixType >::type;
+public:
+  LinearSolveFixture() = default;
+  ~LinearSolveFixture() override = default;
+
+  template< TestMatrixType TEST_MATRIX, int N >
+  void test_matrix_vector_solve() const
+  {
+    static_assert( 1 < N && N <= MAX_SIZE );
+
+    MatrixType A( N, N );
+    TestMatrix< TEST_MATRIX >::create( A.toSlice() );
+    random_permutation( A.toSlice() );
+
+    VectorType b ( N );
+    VectorType x ( N );
+    VectorType x0 ( N );
+
+    // Save selected values of A to check later
+    real64 const a00 = A( 0, 0 );
+    real64 const a01 = A( 0, 1 );
+    real64 const a10 = A( N-1, N-2 );
+    real64 const a11 = A( N-1, N-1 );
+
+    // Create actual solution
+    LvArray::forValuesInSlice( x0.toSlice(), []( real64 & a ){ a = 1.0; } );
+
+    // Multiply to get rhs
+    matrix_vector_multiply( A.toSliceConst(), x0.toSliceConst(), b.toSlice() );
+
+    // Save selected values of b to check later
+    real64 const b0 = b( 0 );
+    real64 const b1 = b( N-1 );
+
+    // Solve
+    BlasLapackLA::solveLinearSystem( A.toSliceConst(), b.toSliceConst(), x.toSlice() );
+
+    // Check solution
+    for( int i = 0; i < N; ++i )
+    {
+      EXPECT_NEAR( x( i ), x0( i ), machinePrecision );
+    }
+
+    // Check that we have not destroyed A
+    EXPECT_NEAR( A( 0, 0 ), a00, machinePrecision );
+    EXPECT_NEAR( A( 0, 1 ), a01, machinePrecision );
+    EXPECT_NEAR( A( N-1, N-2 ), a10, machinePrecision );
+    EXPECT_NEAR( A( N-1, N-1 ), a11, machinePrecision );
+
+    // Check that we have not destroyed b
+    EXPECT_NEAR( b( 0 ), b0, machinePrecision );
+    EXPECT_NEAR( b( N-1 ), b1, machinePrecision );
+  }
+
+  template< TestMatrixType TEST_MATRIX, int N >
+  void test_matrix_vector_solve_inplace( ) const
+  {
+    static_assert( 1 < N && N <= MAX_SIZE );
+
+    MatrixType A( N, N );
+    TestMatrix< TEST_MATRIX >::create( A.toSlice() );
+    random_permutation( A.toSlice() );
+
+    VectorType x ( N );
+    VectorType x0 ( N );
+
+    // Create actual solution
+    LvArray::forValuesInSlice( x0.toSlice(), []( real64 & a ){ a = 1.0; } );
+
+    // Multiply to get rhs
+    matrix_vector_multiply( A.toSliceConst(), x0.toSliceConst(), x.toSlice() );
+
+    // Solve
+    BlasLapackLA::solveLinearSystem( A, x.toSlice() );
+
+    // Check in place
+    for( int i = 0; i < N; ++i )
+    {
+      EXPECT_NEAR( x( i ), x0( i ), machinePrecision );
+    }
+  }
+
+  template< TestMatrixType TEST_MATRIX, int N, int M >
+  void test_matrix_matrix_solve( ) const
+  {
+    static_assert( 1 < N && N <= MAX_SIZE );
+    static_assert( 1 <= M && M <= MAX_SIZE );
+
+    MatrixType A ( N, N );
+    TestMatrix< TEST_MATRIX >::create( A.toSlice() );
+    random_permutation( A.toSlice() );
+
+    MatrixType B ( N, M );
+    MatrixType X ( N, M );
+    MatrixType X0 ( N, M );
+
+    real64 const a00 = A( 0, 0 );
+    real64 const a01 = A( 0, 1 );
+    real64 const a10 = A( N-1, N-2 );
+    real64 const a11 = A( N-1, N-1 );
+
+    // Create actual solution
+    // Populate matrix with random coefficients
+    BlasLapackLA::matrixRand( X0,
+                              BlasLapackLA::RandomNumberDistribution::UNIFORM_m1p1 );
+
+    // Multiply to get rhs
+    matrix_matrix_multiply( A.toSliceConst(), X0.toSliceConst(), B.toSlice() );
+
+    // Save selected values of B to check later
+    real64 const b00 = B( 0, 0 );
+    real64 const b01 = B( N-1, M-1 );
+
+    // Solve
+    BlasLapackLA::solveLinearSystem( A.toSliceConst(), B.toSliceConst(), X.toSlice() );
+
+    // Check
+    for( int i = 0; i < N; ++i )
+    {
+      for( int j = 0; j < M; ++j )
+      {
+        EXPECT_NEAR( X( i, j ), X0( i, j ), machinePrecision );
+      }
+    }
+
+    // Check that we have not destroyed A
+    EXPECT_NEAR( A( 0, 0 ), a00, machinePrecision );
+    EXPECT_NEAR( A( 0, 1 ), a01, machinePrecision );
+    EXPECT_NEAR( A( N-1, N-2 ), a10, machinePrecision );
+    EXPECT_NEAR( A( N-1, N-1 ), a11, machinePrecision );
+
+    // Check that we have not destroyed b
+    EXPECT_NEAR( B( 0, 0 ), b00, machinePrecision );
+    EXPECT_NEAR( B( N-1, M-1 ), b01, machinePrecision );
+  }
+
+  template< TestMatrixType TEST_MATRIX, int N, int M >
+  void test_matrix_matrix_solve_inplace( ) const
+  {
+    static_assert( 1 < N && N <= MAX_SIZE );
+    static_assert( 1 <= M && M <= MAX_SIZE );
+
+    MatrixType A ( N, N );
+    TestMatrix< TEST_MATRIX >::create( A.toSlice() );
+    random_permutation( A.toSlice() );
+
+    MatrixType X ( N, M );
+    MatrixType X0 ( N, M );
+
+    // Create actual solution
+    // Populate matrix with random coefficients
+    BlasLapackLA::matrixRand( X0,
+                              BlasLapackLA::RandomNumberDistribution::UNIFORM_m1p1 );
+
+    // Multiply to get rhs
+    matrix_matrix_multiply( A.toSliceConst(), X0.toSliceConst(), X.toSlice() );
+
+    // Solve
+    BlasLapackLA::solveLinearSystem( A.toSlice(), X.toSlice() );
+
+    // Check
+    for( int i = 0; i < N; ++i )
+    {
+      for( int j = 0; j < M; ++j )
+      {
+        EXPECT_NEAR( X( i, j ), X0( i, j ), machinePrecision );
+      }
+    }
+  }
+};
+
+using LinearSolveTypes = ::testing::Types<
+  Array< real64, 2, MatrixLayout::COL_MAJOR_PERM >,
+  Array< real64, 2, MatrixLayout::ROW_MAJOR_PERM >,
+  StackArray< real64, 2, MAX_SIZE *MAX_SIZE, MatrixLayout::COL_MAJOR_PERM >,
+  StackArray< real64, 2, MAX_SIZE *MAX_SIZE, MatrixLayout::ROW_MAJOR_PERM > >;
+
+class NameGenerator
+{
+public:
+  template< typename T >
+  static std::string GetName( int )
+  {
+    if constexpr (std::is_same_v< T, Array< real64, 2, MatrixLayout::COL_MAJOR_PERM > >) return "col-major-heap-array";
+    if constexpr (std::is_same_v< T, Array< real64, 2, MatrixLayout::ROW_MAJOR_PERM > >) return "row-major-heap-array";
+    if constexpr (std::is_same_v< T, StackArray< real64, 2, MAX_SIZE *MAX_SIZE, MatrixLayout::COL_MAJOR_PERM > >) return "col-major-stack-array";
+    if constexpr (std::is_same_v< T, StackArray< real64, 2, MAX_SIZE *MAX_SIZE, MatrixLayout::ROW_MAJOR_PERM > >) return "row-major-stack-array";
+  }
+};
+
+TYPED_TEST_SUITE( LinearSolveFixture, LinearSolveTypes, NameGenerator );
+
+TYPED_TEST( LinearSolveFixture, matrix_vector_solve_laplace )
+{
+  this->template test_matrix_vector_solve< TestMatrixType::LAPLACE, 5 >();
+  this->template test_matrix_vector_solve_inplace< TestMatrixType::LAPLACE, 5 >();
+  this->template test_matrix_vector_solve< TestMatrixType::LAPLACE, 12 >();
+  this->template test_matrix_vector_solve_inplace< TestMatrixType::LAPLACE, 12 >();
+}
+
+TYPED_TEST( LinearSolveFixture, matrix_vector_solve_grcar )
+{
+  this->template test_matrix_vector_solve< TestMatrixType::GRCAR, 6 >();
+  this->template test_matrix_vector_solve_inplace< TestMatrixType::GRCAR, 6 >();
+  this->template test_matrix_vector_solve< TestMatrixType::GRCAR, 10 >();
+  this->template test_matrix_vector_solve_inplace< TestMatrixType::GRCAR, 10 >();
+}
+
+TYPED_TEST( LinearSolveFixture, matrix_vector_solve_sampling )
+{
+  this->template test_matrix_vector_solve< TestMatrixType::SAMPLING, 3 >();
+  this->template test_matrix_vector_solve_inplace< TestMatrixType::SAMPLING, 3 >();
+  this->template test_matrix_vector_solve< TestMatrixType::SAMPLING, 20 >();
+  this->template test_matrix_vector_solve_inplace< TestMatrixType::SAMPLING, 20 >();
+}
+
+TYPED_TEST( LinearSolveFixture, matrix_matrix_solve_laplace )
+{
+  this->template test_matrix_matrix_solve< TestMatrixType::LAPLACE, 5, 1 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::LAPLACE, 5, 1 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::LAPLACE, 5, 3 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::LAPLACE, 5, 3 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::LAPLACE, 5, 5 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::LAPLACE, 5, 5 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::LAPLACE, 5, 12 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::LAPLACE, 5, 12 >();
+}
+
+TYPED_TEST( LinearSolveFixture, matrix_matrix_solve_grcar )
+{
+  this->template test_matrix_matrix_solve< TestMatrixType::GRCAR, 5, 1 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::GRCAR, 5, 1 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::GRCAR, 5, 3 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::GRCAR, 5, 3 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::GRCAR, 5, 5 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::GRCAR, 5, 5 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::GRCAR, 5, 12 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::GRCAR, 5, 12 >();
+}
+
+TYPED_TEST( LinearSolveFixture, matrix_matrix_solve_sampling )
+{
+  this->template test_matrix_matrix_solve< TestMatrixType::SAMPLING, 5, 1 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::SAMPLING, 5, 1 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::SAMPLING, 12, 3 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::SAMPLING, 5, 3 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::SAMPLING, 5, 5 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::SAMPLING, 5, 5 >();
+  this->template test_matrix_matrix_solve< TestMatrixType::SAMPLING, 5, 12 >();
+  this->template test_matrix_matrix_solve_inplace< TestMatrixType::SAMPLING, 5, 12 >();
+}