diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index e906b6cb..0271f936 100644
--- a/benchmark/CMakeLists.txt
+++ b/benchmark/CMakeLists.txt
@@ -1,13 +1,7 @@
-add_library(tuvx_benchmark_utils OBJECT)
-
-target_compile_features(tuvx_benchmark_utils PUBLIC cxx_std_11)
-
-# link benchmark library to sources in this file
-target_sources(tuvx_benchmark_utils PUBLIC benchmark_tridiagonal_solver.cpp)
-target_link_libraries(tuvx_benchmark_utils PUBLIC musica::tuvx LAPACK::LAPACK
-                                                  benchmark::benchmark)
-target_include_directories(tuvx_benchmark_utils PUBLIC ${LAPACK_INCLUDE_DIRS})
-
 add_executable(benchmark_tridiagonal_solver benchmark_tridiagonal_solver.cpp)
+target_include_directories(
+  benchmark_tridiagonal_solver PUBLIC ${OpenBLAS_INCLUDE_DIRS}
+                                      ${LAPACK_INCLUDE_DIRS})
+
 target_link_libraries(benchmark_tridiagonal_solver
                       PUBLIC LAPACK::LAPACK benchmark::benchmark musica::tuvx)
diff --git a/benchmark/benchmark_tridiagonal_solver.cpp b/benchmark/benchmark_tridiagonal_solver.cpp
index d4da4103..f17a89bb 100644
--- a/benchmark/benchmark_tridiagonal_solver.cpp
+++ b/benchmark/benchmark_tridiagonal_solver.cpp
@@ -1,10 +1,13 @@
-#include "mkl_lapacke.h"
 
 #include <tuvx/linear_algebra/linear_algebra.hpp>
 
 #include <benchmark/benchmark.h>
 
-#include <mkl_lapack.h>
+#ifdef TUVX_COMPILE_WITH_INTEL
+  #include <mkl_lapacke.h>
+#elif TUVX_COMPILE_WITH_GCC
+  #include <lapacke.h>
+#endif
 
 using namespace tuvx;
 typedef TridiagonalMatrix<double> trid_matd;
@@ -13,7 +16,9 @@ typedef std::vector<double> vecd;
 typedef TridiagonalMatrix<float> trid_matf;
 typedef std::vector<float> vecf;
 
-const std::size_t system_size = 1000;
+const std::size_t system_size = 1e6;
+
+const bool diagonally_dominant = false;
 
 static void BM_LAPACKE_SINGLE_PRECISISON(benchmark::State& state)
 {
@@ -21,13 +26,14 @@ static void BM_LAPACKE_SINGLE_PRECISISON(benchmark::State& state)
   vecf b(system_size);
   trid_matf A(system_size);
 
-  FillRandom<float>(A);
-  FillRandom<float>(x);
-  b = Dot<float>(A, x);
-
-  // Perform setup here
   for (auto _ : state)
   {
+    state.PauseTiming();
+    FillRandom<float>(A, diagonally_dominant);
+    FillRandom<float>(x);
+    b = Dot<float>(A, x);
+    state.ResumeTiming();
+
     LAPACKE_sgtsv(
         LAPACK_ROW_MAJOR,
         system_size,
@@ -46,13 +52,13 @@ static void BM_LAPACKE_DOUBLE_PRECISISON(benchmark::State& state)
   vecd b(system_size);
   trid_matd A(system_size);
 
-  FillRandom<double>(A);
-  FillRandom<double>(x);
-  b = Dot<double>(A, x);
-
-  // Perform setup here
   for (auto _ : state)
   {
+    state.PauseTiming();
+    FillRandom<double>(A, diagonally_dominant);
+    FillRandom<double>(x);
+    b = Dot<double>(A, x);
+    state.ResumeTiming();
     LAPACKE_dgtsv(
         LAPACK_ROW_MAJOR,
         system_size,
@@ -70,15 +76,14 @@ static void BM_TUVX_DOUBLE_PRECISISON(benchmark::State& state)
   vecd x(system_size);
   vecd b(system_size);
   trid_matd A(system_size);
-  vecd x_approx(system_size);
-  FillRandom<double>(A);
-  FillRandom<double>(x);
-  b = Dot<double>(A, x);
 
-  // Perform setup here
   for (auto _ : state)
   {
-    x_approx = Solve<double>(A, b);
+    state.PauseTiming();
+    FillRandom<double>(A, diagonally_dominant);
+    FillRandom<double>(x);
+    state.ResumeTiming();
+    Solve<double>(A, b);
   }
 }
 
@@ -87,18 +92,20 @@ static void BM_TUVX_SINGLE_PRECISISON(benchmark::State& state)
   vecf x(system_size);
   vecf b(system_size);
   trid_matf A(system_size);
-  vecf x_approx(system_size);
-  FillRandom<float>(A);
-  FillRandom<float>(x);
-  b = Dot<float>(A, x);
 
   // Perform setup here
   for (auto _ : state)
   {
-    x_approx = Solve<float>(A, b);
+    state.PauseTiming();
+    FillRandom<float>(A, diagonally_dominant);
+    FillRandom<float>(x);
+    b = Dot<float>(A, x);
+    state.ResumeTiming();
+    Solve<float>(A, b);
   }
 }
-// Register the function as a benchmark
+
+// Registering benchmarks
 BENCHMARK(BM_LAPACKE_DOUBLE_PRECISISON);
 BENCHMARK(BM_LAPACKE_SINGLE_PRECISISON);
 BENCHMARK(BM_TUVX_DOUBLE_PRECISISON);
diff --git a/include/tuvx/linear_algebra/linear_algebra.hpp b/include/tuvx/linear_algebra/linear_algebra.hpp
index 33bfd888..7c93d108 100644
--- a/include/tuvx/linear_algebra/linear_algebra.hpp
+++ b/include/tuvx/linear_algebra/linear_algebra.hpp
@@ -37,7 +37,7 @@ namespace tuvx
   /// @brief Fills a matrix with uniformly distributed random values.
   /// @param A tridiagonal matrix to allocate and fill
   template<typename T>
-  void FillRandom(TridiagonalMatrix<T> &A);
+  void FillRandom(TridiagonalMatrix<T> &A, bool diagonally_dominant = false);
 
   /// @fn print vector function
   /// @brief displays the data stored inside a std::vector
@@ -53,12 +53,12 @@ namespace tuvx
   void Print(const TridiagonalMatrix<T> &A);
 
   /// @fn Tridiagonal Linear System Solver
-  /// @brief Thomas' algorithm for solving tridiagonal linear system (A x = b)
+  /// @brief Thomas' algorithm for solving tridiagonal linear system (A x = b).
+  /// store solution in b.
   /// @param A tridiagonal coeffcient matrix
   /// @param b right hand side vector of the tridiagonal system.
-  /// @returns x solution that solves A x = b.
   template<typename T>
-  std::vector<T> Solve(TridiagonalMatrix<T> A, std::vector<T> b);
+  void Solve(TridiagonalMatrix<T> &A, std::vector<T> &b);
 
   /// @fn Tridiagonal Matrix - vector dot product
   /// @brief Specialized dot product function for tridiagonal matrices.
diff --git a/include/tuvx/linear_algebra/linear_algebra.inl b/include/tuvx/linear_algebra/linear_algebra.inl
index 3d67353e..ce75196e 100644
--- a/include/tuvx/linear_algebra/linear_algebra.inl
+++ b/include/tuvx/linear_algebra/linear_algebra.inl
@@ -28,7 +28,7 @@ namespace tuvx
   }
 
   template<typename T>
-  inline std::vector<T> Solve(TridiagonalMatrix<T> A, std::vector<T> b)
+  inline void Solve(TridiagonalMatrix<T> &A, std::vector<T> &b)
   {
     T temp;
     std::size_t N = b.size();
@@ -49,7 +49,6 @@ namespace tuvx
         break;
       }
     }
-    return b;
   }
 
   template<typename T>
@@ -65,19 +64,23 @@ namespace tuvx
   }
 
   template<typename T>
-  inline void FillRandom(TridiagonalMatrix<T> &A)
+  inline void FillRandom(TridiagonalMatrix<T> &A, bool diagonally_dominant)
   {
     FillRandom<T>(A.main_diagonal_);
     FillRandom<T>(A.lower_diagonal_);
     FillRandom<T>(A.upper_diagonal_);
-    // make diagonally dominant (diagonal value greater than sum of its row)
-    std::size_t i = 0;
-    A.main_diagonal_[i] += A.upper_diagonal_[i] + (T)500;
-    for (i = 1; i < A.size_ - 1; i++)
+
+    if (diagonally_dominant)
     {
-      A.main_diagonal_[i] += A.lower_diagonal_[i - 1] + A.upper_diagonal_[i] + (T)500;
+      // make diagonally dominant (diagonal value greater than sum of its row)
+      std::size_t i = 0;
+      A.main_diagonal_[i] += A.upper_diagonal_[i];
+      for (i = 1; i < A.size_ - 1; i++)
+      {
+        A.main_diagonal_[i] += A.lower_diagonal_[i - 1] + A.upper_diagonal_[i];
+      }
+      A.main_diagonal_[i] += A.lower_diagonal_[i - 1];
     }
-    A.main_diagonal_[i] += A.lower_diagonal_[i - 1] + (T)500;
   }
 
   template<typename T>
diff --git a/test/unit/linear_algebra/test_tridiagonal_solver.cpp b/test/unit/linear_algebra/test_tridiagonal_solver.cpp
index cc9c578d..12838b0e 100644
--- a/test/unit/linear_algebra/test_tridiagonal_solver.cpp
+++ b/test/unit/linear_algebra/test_tridiagonal_solver.cpp
@@ -19,7 +19,7 @@ typedef std::vector<float> vecf;
 const double tol_dp = std::numeric_limits<double>::epsilon();
 const float tol_sp = std::numeric_limits<float>::epsilon();
 
-const std::size_t number_of_runs = 100;
+const std::size_t number_of_runs = 20;
 
 std::size_t sizes[5] = { 1000, 10000, 100000, 1000000, 10000000 };
 
@@ -44,9 +44,9 @@ TEST(TridiagSolveTest, SinglePrecision)
       FillRandom<float>(x);
       b = Dot<float>(A, x);
 
-      vecf x_approx = Solve<float>(A, b);
+      Solve<float>(A, b);
 
-      errors[i] += ComputeError<float>(x, x_approx);
+      errors[i] += ComputeError<float>(x, b);
     }
 
     errors[i] /= number_of_runs;
@@ -75,9 +75,9 @@ TEST(TridiagSolveTest, DoublePrecision)
       FillRandom<double>(x);
       b = Dot<double>(A, x);
 
-      vecd x_approx = Solve<double>(A, b);
+      Solve<double>(A, b);
 
-      errors[i] += ComputeError<double>(x, x_approx);
+      errors[i] += ComputeError<double>(x, b);
     }
     errors[i] /= number_of_runs;
     EXPECT_LE(errors[i], tol_dp);