diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1ec9e4d8..db2d537b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,13 +1,12 @@
-################################################################################
+# ##############################################################################
 # Preamble
 
 cmake_minimum_required(VERSION 3.17)
 
 project(
   tuv-x
-  VERSION 0.9.0
-  LANGUAGES Fortran CXX C
-)
+  VERSION 0.8.0
+  LANGUAGES Fortran CXX C)
 
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${PROJECT_SOURCE_DIR}/cmake")
@@ -37,7 +36,7 @@ if(${CMAKE_VERSION} VERSION_LESS "3.21")
   endif()
 endif()
 
-################################################################################
+# ##############################################################################
 # Projet wide setup options
 
 include(CMakeDependentOption)
@@ -45,7 +44,6 @@ option(TUVX_ENABLE_MPI "Enable MPI parallel support" OFF)
 cmake_dependent_option(TUVX_ENABLE_OPENMP "Enable OpenMP support" OFF
                        "TUVX_ENABLE_MPI" OFF)
 option(TUVX_ENABLE_TESTS "Build tests" ON)
-option(TUVX_ENABLE_BENCHMARK "Build benchmark examples" ON)
 option(TUVX_ENABLE_COVERAGE "Enable code coverage output" OFF)
 option(TUVX_ENABLE_MEMCHECK "Enable memory checking in tests" OFF)
 option(TUVX_ENABLE_NC_CONFIG "Use nc-config to determine NetCDF libraries" OFF)
@@ -96,12 +94,12 @@ if(TUVX_ENABLE_TESTS)
             ${CMAKE_CURRENT_SOURCE_DIR}/examples ${CMAKE_BINARY_DIR}/examples)
 endif()
 
-################################################################################
+# ##############################################################################
 # Dependencies
 
 include(cmake/dependencies.cmake)
 
-################################################################################
+# ##############################################################################
 # TUV-x targets
 
 add_subdirectory(src)
@@ -117,11 +115,8 @@ if(NOT BUILD_SHARED_LIBS)
   set_target_properties(tuv-x PROPERTIES POSITION_INDEPENDENT_CODE ON)
 endif()
 
-target_link_libraries(tuv-x 
-  PUBLIC 
-    musica::tuvx
-    yaml-cpp::yaml-cpp
-)
+target_link_libraries(tuv-x PUBLIC musica::tuvx yaml-cpp::yaml-cpp
+                                   ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES})
 
 target_include_directories(
   tuv-x PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/include>
@@ -131,14 +126,14 @@ if(TUVX_ENABLE_OPENMP)
   target_link_libraries(tuv-x PUBLIC OpenMP::OpenMP_Fortran)
 endif()
 
-################################################################################
+# ##############################################################################
 # TUV-x docs
 
 if(TUVX_BUILD_DOCS)
   add_subdirectory(docs)
 endif()
 
-################################################################################
+# ##############################################################################
 # TUV-x tests
 
 if(PROJECT_IS_TOP_LEVEL AND TUVX_ENABLE_TESTS)
@@ -160,9 +155,10 @@ if(PROJECT_IS_TOP_LEVEL AND TUVX_ENABLE_TESTS)
 
   enable_testing()
   add_subdirectory(test)
+  add_subdirectory(benchmark)
 endif()
 
-################################################################################
+# ##############################################################################
 # Packaging
 
 # only include packaging if we are the top level project being built
@@ -170,12 +166,4 @@ if(PROJECT_IS_TOP_LEVEL)
   add_subdirectory(packaging)
 endif()
 
-################################################################################
-# benchmarking
-
-if(TUVX_ENABLE_BENCHMARK)
-  add_subdirectory(benchmark)
-endif()
-
-
-################################################################################
+# ##############################################################################
diff --git a/benchmark/benchmark_tridiagonal_solver.cpp b/benchmark/benchmark_tridiagonal_solver.cpp
index e1ea6e3f..5779728d 100644
--- a/benchmark/benchmark_tridiagonal_solver.cpp
+++ b/benchmark/benchmark_tridiagonal_solver.cpp
@@ -3,35 +3,33 @@
 
 #include <benchmark/benchmark.h>
 
+#include <random>
+#include <vector>
+
 #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;
-typedef std::vector<double> vecd;
-
-typedef TridiagonalMatrix<float> trid_matf;
-typedef std::vector<float> vecf;
-
 const std::size_t system_size = 1e6;
 
-const bool diagonally_dominant = false;
+const bool diagonally_dominant = true;
+
+const unsigned random_number_seed = 1;
 
 static void BM_LAPACKE_SINGLE_PRECISISON(benchmark::State& state)
 {
-  vecf x(system_size);
-  vecf b(system_size);
-  trid_matf A(system_size);
-
+  std::vector<float> x(system_size);
+  std::vector<float> b(system_size);
+  tuvx::TridiagonalMatrix<float> A(system_size);
+  std::mt19937 random_device(random_number_seed);
   for (auto _ : state)
   {
     state.PauseTiming();
-    FillRandom<float>(A, diagonally_dominant);
-    FillRandom<float>(x);
-    b = Dot<float>(A, x);
+    tuvx::FillRandom<float>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<float>(x, random_number_seed);
+    b = tuvx::Dot<float>(A, x);
     state.ResumeTiming();
 
     LAPACKE_sgtsv(
@@ -48,18 +46,20 @@ static void BM_LAPACKE_SINGLE_PRECISISON(benchmark::State& state)
 
 static void BM_LAPACKE_DOUBLE_PRECISISON(benchmark::State& state)
 {
-  vecd x(system_size);
-  vecd b(system_size);
-  trid_matd A(system_size);
+  std::vector<double> x(system_size);
+  std::vector<double> b(system_size);
+  tuvx::TridiagonalMatrix<double> A(system_size);
 
+  std::mt19937 random_device(random_number_seed);
   // Perform setup here
   for (auto _ : state)
   {
     state.PauseTiming();
-    FillRandom<double>(A, diagonally_dominant);
-    FillRandom<double>(x);
-    b = Dot<double>(A, x);
+    tuvx::FillRandom<double>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<double>(x, random_number_seed);
+    b = tuvx::Dot<double>(A, x);
     state.ResumeTiming();
+
     LAPACKE_dgtsv(
         LAPACK_ROW_MAJOR,
         system_size,
@@ -74,37 +74,37 @@ static void BM_LAPACKE_DOUBLE_PRECISISON(benchmark::State& state)
 
 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);
+  std::vector<double> x(system_size);
+  std::vector<double> b(system_size);
+  tuvx::TridiagonalMatrix<double> A(system_size);
+  std::vector<double> x_approx(system_size);
 
   for (auto _ : state)
   {
     state.PauseTiming();
-    FillRandom<double>(A, diagonally_dominant);
-    FillRandom<double>(x);
+    tuvx::FillRandom<double>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<double>(x, random_number_seed);
     state.ResumeTiming();
-    Solve<double>(A, b);
+    tuvx::Solve<double>(A, b);
   }
 }
 
 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);
+  std::vector<float> x(system_size);
+  std::vector<float> b(system_size);
+  tuvx::TridiagonalMatrix<float> A(system_size);
+  std::vector<float> x_approx(system_size);
 
   // Perform setup here
   for (auto _ : state)
   {
     state.PauseTiming();
-    FillRandom<float>(A, diagonally_dominant);
-    FillRandom<float>(x);
-    b = Dot<float>(A, x);
+    tuvx::FillRandom<float>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<float>(x, random_number_seed);
+    b = tuvx::Dot<float>(A, x);
     state.ResumeTiming();
-    Solve<float>(A, b);
+    tuvx::Solve<float>(A, b);
   }
 }
 
diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake
index d41d95e4..d451d662 100644
--- a/cmake/dependencies.cmake
+++ b/cmake/dependencies.cmake
@@ -1,14 +1,14 @@
 find_package(PkgConfig REQUIRED)
 include(FetchContent)
 
-################################################################################
+# ##############################################################################
 # LAPACK
 
 find_package(BLAS)
 find_package(LAPACK)
 find_package(LAPACKE)
 
-################################################################################
+# ##############################################################################
 # Memory check
 
 if(TUVX_ENABLE_MEMCHECK)
@@ -28,7 +28,7 @@ if(TUVX_ENABLE_MEMCHECK)
   endif()
 endif()
 
-################################################################################
+# ##############################################################################
 # OpenMP
 
 if(TUVX_ENABLE_OPENMP)
@@ -41,7 +41,7 @@ if(TUVX_ENABLE_OPENMP)
   endif()
 endif()
 
-################################################################################
+# ##############################################################################
 # NetCDF library
 
 find_package(PkgConfig REQUIRED)
@@ -49,7 +49,7 @@ find_package(PkgConfig REQUIRED)
 pkg_check_modules(netcdff IMPORTED_TARGET REQUIRED netcdf-fortran)
 pkg_check_modules(netcdfc IMPORTED_TARGET REQUIRED netcdf)
 
-################################################################################
+# ##############################################################################
 # yaml-cpp
 
 FetchContent_Declare(
@@ -60,7 +60,7 @@ FetchContent_Declare(
   ${FETCHCONTENT_QUIET})
 FetchContent_MakeAvailable(yaml-cpp)
 
-################################################################################
+# ##############################################################################
 # Docs
 
 if(TUVX_BUILD_DOCS)
@@ -68,15 +68,14 @@ if(TUVX_BUILD_DOCS)
   find_package(Sphinx REQUIRED)
 endif()
 
-################################################################################
+# ##############################################################################
 # google test and benchmark
 
 if(TUVX_ENABLE_BENCHMARK)
   FetchContent_Declare(
     googlebenchmark
     GIT_REPOSITORY https://github.com/google/benchmark.git
-    GIT_TAG v1.8.3
-    )
+    GIT_TAG v1.8.3)
 
   set(BENCHMARK_DOWNLOAD_DEPENDENCIES ON)
   set(BENCHMARK_ENABLE_GTEST_TESTS OFF)
@@ -105,4 +104,4 @@ if(TUVX_ENABLE_TESTS)
   set_target_properties(gtest_main PROPERTIES CXX_CLANG_TIDY "")
 endif()
 
-################################################################################
\ No newline at end of file
+# ##############################################################################
diff --git a/include/tuvx/linear_algebra/linear_algebra.hpp b/include/tuvx/linear_algebra/linear_algebra.hpp
index 7c93d108..0358cd04 100644
--- a/include/tuvx/linear_algebra/linear_algebra.hpp
+++ b/include/tuvx/linear_algebra/linear_algebra.hpp
@@ -31,13 +31,13 @@ namespace tuvx
   /// @brief Fill a std::vector with random values
   /// @param x Vector to allocate and fill
   template<typename T>
-  void FillRandom(std::vector<T> &x);
+  void FillRandom(std::vector<T> &x, const unsigned &seed);
 
   /// @fn Initialize random matrix function
   /// @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, bool diagonally_dominant = false);
+  void FillRandom(TridiagonalMatrix<T> &A, const unsigned &seed, bool diagonally_dominant = false);
 
   /// @fn print vector function
   /// @brief displays the data stored inside a std::vector
diff --git a/include/tuvx/linear_algebra/linear_algebra.inl b/include/tuvx/linear_algebra/linear_algebra.inl
index 52cf6085..385e639e 100644
--- a/include/tuvx/linear_algebra/linear_algebra.inl
+++ b/include/tuvx/linear_algebra/linear_algebra.inl
@@ -1,4 +1,5 @@
 
+#include <random>
 namespace tuvx
 {
 
@@ -52,24 +53,24 @@ namespace tuvx
   }
 
   template<typename T>
-  inline void FillRandom(std::vector<T> &x)
+  inline void FillRandom(std::vector<T> &x, const unsigned &seed)
   {
-    std::random_device dev;
-    std::mt19937 rng(dev());
     // sample from normal distribution
+
+    std::mt19937 random_device(seed);
     std::normal_distribution<double> distribution(5.0, 1.0);
     for (std::size_t i = 0; i < x.size(); i++)
     {
-      x[i] = (T)distribution(rng);
+      x[i] = (T)distribution(random_device);
     }
   }
 
   template<typename T>
-  inline void FillRandom(TridiagonalMatrix<T> &A, bool diagonally_dominant)
+  inline void FillRandom(TridiagonalMatrix<T> &A, const unsigned &seed, bool diagonally_dominant)
   {
-    FillRandom<T>(A.main_diagonal_);
-    FillRandom<T>(A.lower_diagonal_);
-    FillRandom<T>(A.upper_diagonal_);
+    FillRandom<T>(A.main_diagonal_, seed);
+    FillRandom<T>(A.lower_diagonal_, seed);
+    FillRandom<T>(A.upper_diagonal_, seed);
 
     if (diagonally_dominant)
     {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 5c9ad8f7..7981493b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,4 +1,4 @@
-################################################################################
+# ##############################################################################
 # Photo-decomp tool source
 
 # object library
@@ -44,19 +44,15 @@ endif()
 add_library(tuvx $<TARGET_OBJECTS:tuvx_object>)
 add_library(musica::tuvx ALIAS tuvx)
 
-
 if(NOT BUILD_SHARED_LIBS)
-  set_target_properties(tuvx_object
-    PROPERTIES
-    POSITION_INDEPENDENT_CODE ON
-  )
+  set_target_properties(tuvx_object PROPERTIES POSITION_INDEPENDENT_CODE ON)
 endif()
 
-set_target_properties(tuvx PROPERTIES 
-  ARCHIVE_OUTPUT_DIRECTORY ${TUVX_LIB_DIR}
-  VERSION ${PROJECT_VERSION}
-  SOVERSION ${PROJECT_VERSION_MAJOR}
-)
+set_target_properties(
+  tuvx
+  PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${TUVX_LIB_DIR}
+             VERSION ${PROJECT_VERSION}
+             SOVERSION ${PROJECT_VERSION_MAJOR})
 
 target_link_libraries(tuvx PUBLIC tuvx_object)
 
@@ -112,4 +108,4 @@ add_subdirectory(quantum_yields)
 add_subdirectory(radiative_transfer)
 add_subdirectory(spectral_weights)
 add_subdirectory(util)
-################################################################################
+# ##############################################################################
diff --git a/test/unit/linear_algebra/test_tridiagonal_solver.cpp b/test/unit/linear_algebra/test_tridiagonal_solver.cpp
index d3993562..8eaae0cd 100644
--- a/test/unit/linear_algebra/test_tridiagonal_solver.cpp
+++ b/test/unit/linear_algebra/test_tridiagonal_solver.cpp
@@ -5,6 +5,7 @@
 #include <cfloat>
 #include <cstddef>
 #include <cstdlib>
+#include <iterator>
 #include <limits>
 #include <vector>
 
@@ -14,22 +15,14 @@
   #include <lapacke.h>
 #endif
 
-using namespace tuvx;
-
-typedef TridiagonalMatrix<double> trid_matd;
-typedef std::vector<double> vecd;
-
-typedef TridiagonalMatrix<float> trid_matf;
-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 = 20;
+const std::size_t size = 10;
 const bool diagonally_dominant = true;
 
-std::size_t sizes[5] = { 1000, 10000, 100000, 1000000, 10000000 };
-
+const unsigned random_number_seed = 1;
 /// @test Tridiagonal Solver Test for single Precision Floats.
 /// @brief Generate random tridiagonal matrix $A$ and vector $x$,
 /// compute $b=A \cdot x$, and check if solution is reconstructed
@@ -38,28 +31,24 @@ std::size_t sizes[5] = { 1000, 10000, 100000, 1000000, 10000000 };
 TEST(TridiagSolveTest, SinglePrecision)
 {
   float error = 0;
-  vecd errors(5, 0);
-  for (std::size_t i = 0; i < 5; i++)
-  {
-    for (std::size_t j = 0; j < number_of_runs; j++)
-    {
-      vecf x(sizes[i]);
-      vecf b(sizes[i]);
-      trid_matf A(sizes[i]);
 
-      FillRandom<float>(A, diagonally_dominant);
-      FillRandom<float>(x);
-      b = Dot<float>(A, x);
+  for (std::size_t j = 0; j < number_of_runs; j++)
+  {
+    std::vector<float> x(size);
+    std::vector<float> b(size);
+    tuvx::TridiagonalMatrix<float> A(size);
 
-      Solve<float>(A, b);
+    tuvx::FillRandom<float>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<float>(x, random_number_seed);
+    b = tuvx::Dot<float>(A, x);
+    tuvx::Solve<float>(A, b);
 
-      errors[i] += ComputeError<float>(x, b);
-    }
+    error += tuvx::ComputeError<float>(x, b);
+  }
 
-    errors[i] /= number_of_runs;
+  error /= number_of_runs;
 
-    EXPECT_LE(errors[i], tol_sp);
-  }
+  EXPECT_LE(error, tol_sp);
 }
 
 /// @test Tridiagonal Solver Test for Double Precision Floats.
@@ -69,26 +58,22 @@ TEST(TridiagSolveTest, SinglePrecision)
 /// sizes to check consistency
 TEST(TridiagSolveTest, DoublePrecision)
 {
-  vecd errors(5, 0);
-  for (std::size_t i = 0; i < 5; i++)
+  double error = 0;
+  for (std::size_t j = 0; j < number_of_runs; j++)
   {
-    for (std::size_t j = 0; j < number_of_runs; j++)
-    {
-      vecd x(sizes[i]);
-      vecd b(sizes[i]);
-      trid_matd A(sizes[i]);
-
-      FillRandom<double>(A, diagonally_dominant);
-      FillRandom<double>(x);
-      b = Dot<double>(A, x);
-
-      Solve<double>(A, b);
-
-      errors[i] += ComputeError<double>(x, b);
-    }
-    errors[i] /= number_of_runs;
-    EXPECT_LE(errors[i], tol_dp);
+    std::vector<double> x(size);
+    std::vector<double> b(size);
+    tuvx::TridiagonalMatrix<double> A(size);
+
+    tuvx::FillRandom<double>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<double>(x, random_number_seed);
+    b = tuvx::Dot<double>(A, x);
+    tuvx::Solve<double>(A, b);
+
+    error += tuvx::ComputeError<double>(x, b);
   }
+  error /= number_of_runs;
+  EXPECT_LE(error, tol_dp);
 }
 
 /// @test LAPACKE Tridiagonal Solver Test for single Precision Floats.
@@ -98,36 +83,26 @@ TEST(TridiagSolveTest, DoublePrecision)
 /// sizes to check consistency
 TEST(LapackeTest, SinglePrecision)
 {
-  vecf errors(5, 0);
+  float error;
 
-  for (std::size_t i = 0; i < 5; i++)
+  for (std::size_t j = 0; j < number_of_runs; j++)
   {
-    for (std::size_t j = 0; j < number_of_runs; j++)
-    {
-      vecf x(sizes[i]);
-      vecf b(sizes[i]);
-      trid_matf A(sizes[i]);
-
-      FillRandom<float>(A, diagonally_dominant);
-      FillRandom<float>(x);
-      b = Dot<float>(A, x);
-
-      LAPACKE_sgtsv(
-          LAPACK_ROW_MAJOR,
-          sizes[i],
-          1,
-          A.lower_diagonal_.data(),
-          A.main_diagonal_.data(),
-          A.upper_diagonal_.data(),
-          b.data(),
-          1);
-
-      // to be written to a file
-      errors[i] += ComputeError<float>(x, b);
-    }
-    errors[i] /= number_of_runs;
-    EXPECT_LE(errors[i], tol_sp);
+    std::vector<float> x(size);
+    std::vector<float> b(size);
+    tuvx::TridiagonalMatrix<float> A(size);
+
+    tuvx::FillRandom<float>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<float>(x, random_number_seed);
+    b = tuvx::Dot<float>(A, x);
+
+    LAPACKE_sgtsv(
+        LAPACK_ROW_MAJOR, size, 1, A.lower_diagonal_.data(), A.main_diagonal_.data(), A.upper_diagonal_.data(), b.data(), 1);
+
+    // to be written to a file
+    error += tuvx::ComputeError<float>(x, b);
   }
+  error /= number_of_runs;
+  EXPECT_LE(error, tol_sp);
 }
 
 /// @test LAPACKE Tridiagonal Solver Test for double Precision Floats.
@@ -137,33 +112,22 @@ TEST(LapackeTest, SinglePrecision)
 /// sizes to check consistency
 TEST(LapackeTest, DoublePrecision)
 {
-  vecd errors(5, 0);
-  vecd times(5, 0);
-  for (std::size_t i = 0; i < 5; i++)
+  double error;
+  for (std::size_t j = 0; j < number_of_runs; j++)
   {
-    for (std::size_t j = 0; j < number_of_runs; j++)
-    {
-      vecd x(sizes[i]);
-      vecd b(sizes[i]);
-      trid_matd A(sizes[i]);
-
-      FillRandom<double>(A, diagonally_dominant);
-      FillRandom<double>(x);
-      b = Dot<double>(A, x);
-
-      LAPACKE_dgtsv(
-          LAPACK_ROW_MAJOR,
-          sizes[i],
-          1,
-          A.lower_diagonal_.data(),
-          A.main_diagonal_.data(),
-          A.upper_diagonal_.data(),
-          b.data(),
-          1);
-
-      errors[i] += ComputeError<double>(x, b);
-    }
-    errors[i] /= number_of_runs;
-    EXPECT_LE(errors[i], tol_dp);
+    std::vector<double> x(size);
+    std::vector<double> b(size);
+    tuvx::TridiagonalMatrix<double> A(size);
+
+    tuvx::FillRandom<double>(A, random_number_seed, diagonally_dominant);
+    tuvx::FillRandom<double>(x, random_number_seed);
+    b = tuvx::Dot<double>(A, x);
+
+    LAPACKE_dgtsv(
+        LAPACK_ROW_MAJOR, size, 1, A.lower_diagonal_.data(), A.main_diagonal_.data(), A.upper_diagonal_.data(), b.data(), 1);
+
+    error += tuvx::ComputeError<double>(x, b);
   }
+  error /= number_of_runs;
+  EXPECT_LE(error, tol_dp);
 }