Repository: ssciwr/hammingdist (language: python)

Issue:
Add CUDA implementation of distance function for NVIDIA GPUs
possible strategies
- naive implementation assuming everything fits into gpu ram
  - if we restrict to linux and sm60+ we can use [unified memory](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#unified-memory-programming) which allows memory oversubscription
  - should be relatively straightforward to implement
  - likely very poor scaling for large datasets
- assuming the gene vectors fit into gpu ram
  - split distances matrix into nxn sub-matrices, where one sub-matrix fits in gpu ram along with the gene vectors
  - calculate one sub-matrix at a time & copy each to cpu when done
  - should scale better
- assuming even gene vectors don't fit into gpu ram
  - also need to block the gene vectors, then sum all contributions


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 54c4dd0..68049b5 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -65,4 +65,4 @@ jobs:
       - name: run tests
         shell: bash
         working-directory: ${{runner.workspace}}/build
-        run: ctest
+        run: ctest -j2 --rerun-failed --output-on-failure
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index 58eae2b..7e3ee38 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -27,9 +27,11 @@ jobs:
         with:
           submodules: "recursive"
       - uses: pypa/cibuildwheel@v2.16
+        env:
+          CIBW_MANYLINUX_X86_64_IMAGE: sameli/manylinux2014_x86_64_cuda_11.8
       - uses: actions/upload-artifact@v4
         with:
-          name: wheels-${{ matrix.os }}
+          name: artifacts-${{ matrix.os }}
           path: ./wheelhouse/*.whl
 
   upload_pypi:
@@ -42,7 +44,7 @@ jobs:
     steps:
       - uses: actions/download-artifact@v4
         with:
-          pattern: wheels-*
+          pattern: artifacts-*
           merge-multiple: true
           path: dist
 
diff --git a/.gitmodules b/.gitmodules
index f88f44d..0e6d3ab 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -13,3 +13,6 @@
 [submodule "Catch2"]
 	path = ext/Catch2
 	url = https://github.com/catchorg/Catch2.git
+[submodule "ext/fmt"]
+	path = ext/fmt
+	url = https://github.com/fmtlib/fmt.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8919e23..69984d6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,8 +1,8 @@
-cmake_minimum_required(VERSION 3.11..3.27)
+cmake_minimum_required(VERSION 3.23..3.27)
 
 project(
   hammingdist
-  VERSION 0.21.0
+  VERSION 1.0.0
   LANGUAGES CXX)
 
 include(CTest)
@@ -39,6 +39,14 @@ set(HAMMING_WITH_NEON
     no
     CACHE BOOL "Enable NEON optimized code on Arm64 CPUs")
 
+set(HAMMING_WITH_CUDA
+    no
+    CACHE BOOL "Build with CUDA support (nvidia gpu)")
+
+if(HAMMING_WITH_CUDA)
+  enable_language(CUDA)
+endif()
+
 # Add git submodules
 add_subdirectory(ext)
 
diff --git a/README.md b/README.md
index 4719e67..061a51e 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-A small C++ tool to calculate pairwise distances between gene sequences given in fasta format.
+A small and fast C++ tool to calculate pairwise distances between gene sequences given in fasta format.
 
 [![DOI](https://zenodo.org/badge/308676358.svg)](https://zenodo.org/badge/latestdoi/308676358)
 [![pypi releases](https://img.shields.io/pypi/v/hammingdist.svg)](https://pypi.org/project/hammingdist)
@@ -117,9 +117,34 @@ distance = hammingdist.distance("ACGTX", "AAGTX", include_x=True)
 
 ## OpenMP on linux
 
-The latest versions of hammingdist on linux are now built with OpenMP (multithreading) support.
-If this causes any issues, you can install a previous version of hammingdist without OpenMP support:
+On linux hammingdist is built with OpenMP (multithreading) support, and will automatically make use of all available CPU threads.
 
-```bash
-pip install hammingdist==0.11.0
+## CUDA on linux
+
+On linux hammingdist is also built with CUDA (Nvidia GPU) support.
+To use the GPU instead of the CPU, set `use_gpu=True` when calling `from_fasta`:
+
+```python
+import hammingdist
+
+data = hammingdist.from_fasta("example.fasta", use_gpu=True)
+```
+
+Additionally, the lower triangular matrix file can now be directly constructed from the fasta file
+using the GPU with the `from_fasta_to_lower_triangular` function.
+This avoids storing the entire distances matrix in memory and interleaves computation on the GPU with disk I/O on the CPU,
+which means it requires less RAM and runs faster.
+
+```python
+import hammingdist
+
+hammingdist.from_fasta_to_lower_triangular('input_fasta.txt', 'output_lower_triangular.txt', use_gpu=True)
 ```
+
+![overview](plots/speed.png)
+
+## Performance history
+
+A rough measure of the impact of the different performance improvements in hammingdist:
+
+![overview](plots/overview.png)
diff --git a/ext/CMakeLists.txt b/ext/CMakeLists.txt
index 6a00470..897253e 100644
--- a/ext/CMakeLists.txt
+++ b/ext/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory(fmt)
+
 if(HAMMING_BUILD_PYTHON)
   add_subdirectory(pybind11)
 endif()
diff --git a/ext/Catch2 b/ext/Catch2
index 182c910..f981c9c 160000
--- a/ext/Catch2
+++ b/ext/Catch2
@@ -1,1 +1,1 @@
-Subproject commit 182c910b4b63ff587a3440e08f84f70497e49a81
+Subproject commit f981c9cbcac07a2690e5a86767eba490b5465463
diff --git a/ext/fmt b/ext/fmt
new file mode 160000
index 0000000..f5e5435
--- /dev/null
+++ b/ext/fmt
@@ -0,0 +1,1 @@
+Subproject commit f5e54359df4c26b6230fc61d38aa294581393084
diff --git a/include/hamming/distance_avx2.hh b/include/hamming/distance_avx2.hh
index 8fc6cdc..7823a52 100644
--- a/include/hamming/distance_avx2.hh
+++ b/include/hamming/distance_avx2.hh
@@ -1,5 +1,4 @@
-#ifndef _HAMMING_DISTANCE_AVX2_HH
-#define _HAMMING_DISTANCE_AVX2_HH
+#pragma once
 
 #include <cstdint>
 #include <vector>
@@ -12,5 +11,3 @@ int distance_avx2(const std::vector<GeneBlock> &a,
                   const std::vector<GeneBlock> &b);
 
 }
-
-#endif
diff --git a/include/hamming/distance_avx512.hh b/include/hamming/distance_avx512.hh
index 919fd31..a21d414 100644
--- a/include/hamming/distance_avx512.hh
+++ b/include/hamming/distance_avx512.hh
@@ -1,5 +1,4 @@
-#ifndef _HAMMING_DISTANCE_AVX512_HH
-#define _HAMMING_DISTANCE_AVX512_HH
+#pragma once
 
 #include <cstdint>
 #include <vector>
@@ -12,5 +11,3 @@ int distance_avx512(const std::vector<GeneBlock> &a,
                     const std::vector<GeneBlock> &b);
 
 }
-
-#endif
diff --git a/include/hamming/distance_cuda.hh b/include/hamming/distance_cuda.hh
new file mode 100644
index 0000000..c7c7e52
--- /dev/null
+++ b/include/hamming/distance_cuda.hh
@@ -0,0 +1,27 @@
+#pragma once
+
+#include <cstdint>
+#include <iostream>
+#include <vector>
+
+#include "hamming/hamming_impl_types.hh"
+
+namespace hamming {
+
+bool distance_cuda_have_device();
+
+int distance_cuda(const std::vector<GeneBlock> &a,
+                  const std::vector<GeneBlock> &b);
+
+// for now explicit function def for each choice of integer type
+std::vector<uint8_t>
+distances_cuda_8bit(const std::vector<std::vector<GeneBlock>> &data);
+
+std::vector<uint16_t>
+distances_cuda_16bit(const std::vector<std::vector<GeneBlock>> &data);
+
+void distances_cuda_to_lower_triangular(
+    const std::vector<std::vector<GeneBlock>> &data,
+    const std::string &filename);
+
+} // namespace hamming
diff --git a/include/hamming/distance_sse2.hh b/include/hamming/distance_sse2.hh
index 6920226..ac5d2c7 100644
--- a/include/hamming/distance_sse2.hh
+++ b/include/hamming/distance_sse2.hh
@@ -1,5 +1,4 @@
-#ifndef _HAMMING_DISTANCE_SSE2_HH
-#define _HAMMING_DISTANCE_SSE2_HH
+#pragma once
 
 #include <cstdint>
 #include <vector>
@@ -12,5 +11,3 @@ int distance_sse2(const std::vector<GeneBlock> &a,
                   const std::vector<GeneBlock> &b);
 
 }
-
-#endif
diff --git a/include/hamming/hamming.hh b/include/hamming/hamming.hh
index 61af9b7..06f90a0 100644
--- a/include/hamming/hamming.hh
+++ b/include/hamming/hamming.hh
@@ -1,13 +1,12 @@
-#ifndef _HAMMING_HH
-#define _HAMMING_HH
+#pragma once
 
 #include "hamming/hamming_impl.hh"
 #include "hamming/hamming_types.hh"
+#include "hamming/hamming_utils.hh"
 #include <cmath>
 #include <fstream>
 #include <sstream>
 #include <string>
-#include <unordered_map>
 #include <vector>
 
 namespace hamming {
@@ -19,10 +18,11 @@ inline std::size_t uint_sqrt(std::size_t x) {
 template <typename DistIntType> struct DataSet {
   explicit DataSet(std::vector<std::string> &data, bool include_x = false,
                    bool clear_input_data = false,
-                   std::vector<std::size_t> &&indices = {})
+                   std::vector<std::size_t> &&indices = {},
+                   bool use_gpu = false)
       : nsamples(data.size()), sequence_indices(std::move(indices)) {
     validate_data(data);
-    result = distances<DistIntType>(data, include_x, clear_input_data);
+    result = distances<DistIntType>(data, include_x, clear_input_data, use_gpu);
   }
 
   explicit DataSet(const std::string &filename) {
@@ -67,35 +67,7 @@ template <typename DistIntType> struct DataSet {
   }
 
   void dump_lower_triangular(const std::string &filename) {
-    std::ofstream stream(filename);
-#ifdef HAMMING_WITH_OPENMP
-    std::size_t block_size = 200;
-#pragma omp parallel for ordered schedule(static, 1)
-    for (std::size_t i_start = 1; i_start < nsamples; i_start += block_size) {
-      std::stringstream line;
-      std::size_t i_end = i_start + block_size;
-      if (i_end > nsamples) {
-        i_end = nsamples;
-      }
-      for (std::size_t i = i_start; i < i_end; ++i) {
-        std::size_t offset{i * (i - 1) / 2};
-        for (std::size_t j = 0; j + 1 < i; ++j) {
-          line << static_cast<int>(result[offset + j]) << ",";
-        }
-        line << static_cast<int>(result[offset + i - 1]) << "\n";
-      }
-#pragma omp ordered
-      stream << line.str();
-    }
-#else
-    std::size_t k = 0;
-    for (std::size_t i = 1; i < nsamples; ++i) {
-      for (std::size_t j = 0; j + 1 < i; ++j) {
-        stream << static_cast<int>(result[k++]) << ",";
-      }
-      stream << static_cast<int>(result[k++]) << "\n";
-    }
-#endif
+    write_lower_triangular(filename, result);
   }
 
   void dump_sequence_indices(const std::string &filename) {
@@ -127,67 +99,52 @@ template <typename DistIntType> struct DataSet {
   std::vector<std::size_t> sequence_indices{};
 };
 
-DataSet<DefaultDistIntType> from_stringlist(std::vector<std::string> &data);
+DataSet<DefaultDistIntType> from_stringlist(std::vector<std::string> &data,
+                                            bool include_x = false,
+                                            bool use_gpu = false);
 
 DataSet<DefaultDistIntType> from_csv(const std::string &filename);
 
-DataSet<DefaultDistIntType> from_lower_triangular(const std::string &filename);
-
 template <typename DistIntType>
-DataSet<DistIntType>
-from_fasta(const std::string &filename, bool include_x = false,
-           bool remove_duplicates = false, std::size_t n = 0) {
-  std::vector<std::string> data;
-  data.reserve(n);
-  if (n == 0) {
-    n = std::numeric_limits<std::size_t>::max();
-    data.reserve(65536);
-  }
-  std::unordered_map<std::string, std::size_t> map_seq_to_index;
-  std::vector<std::size_t> sequence_indices{};
-  // Initializing the stream
+DataSet<DistIntType> from_lower_triangular(const std::string &filename) {
+  std::vector<DistIntType> distances;
   std::ifstream stream(filename);
-  std::size_t count = 0;
-  std::size_t count_unique = 0;
   std::string line;
-  // skip first header
-  std::getline(stream, line);
-  while (count < n && !stream.eof()) {
-    std::string seq{};
-    while (std::getline(stream, line) && line[0] != '>') {
-      seq.append(line);
-    }
-    if (remove_duplicates) {
-      auto result = map_seq_to_index.emplace(std::move(seq), count_unique);
-      if (result.second) {
-        ++count_unique;
-      }
-      sequence_indices.push_back(result.first->second);
-    } else {
-      data.push_back(std::move(seq));
-    }
-    ++count;
-  }
-  if (remove_duplicates) {
-    // copy each unique sequence to the vector of strings
-    data.resize(count_unique);
-    for (auto &key_value_pair : map_seq_to_index) {
-      data[key_value_pair.second] = key_value_pair.first;
+  while (std::getline(stream, line)) {
+    std::istringstream s(line);
+    std::string d;
+    while (s.good()) {
+      std::getline(s, d, ',');
+      distances.push_back(safe_int_cast<DistIntType>(std::stoi(d)));
     }
   }
+  return DataSet(std::move(distances));
+}
+
+template <typename DistIntType>
+DataSet<DistIntType> from_fasta(const std::string &filename,
+                                bool include_x = false,
+                                bool remove_duplicates = false,
+                                std::size_t n = 0, bool use_gpu = false) {
+  auto [data, sequence_indices] = read_fasta(filename, remove_duplicates, n);
   return DataSet<DistIntType>(data, include_x, true,
-                              std::move(sequence_indices));
+                              std::move(sequence_indices), use_gpu);
 }
 
+void from_fasta_to_lower_triangular(const std::string &input_filename,
+                                    const std::string &output_filename,
+                                    bool remove_duplicates = false,
+                                    std::size_t n = 0, bool use_gpu = false);
+
 ReferenceDistIntType distance(const std::string &seq0, const std::string &seq1,
                               bool include_x = false);
+
 std::vector<ReferenceDistIntType>
 fasta_reference_distances(const std::string &reference_sequence,
                           const std::string &fasta_file,
                           bool include_x = false);
+
 std::vector<std::size_t> fasta_sequence_indices(const std::string &fasta_file,
                                                 std::size_t n = 0);
 
 } // namespace hamming
-
-#endif
diff --git a/include/hamming/hamming_impl.hh b/include/hamming/hamming_impl.hh
index 4e0a97e..6f1814f 100644
--- a/include/hamming/hamming_impl.hh
+++ b/include/hamming/hamming_impl.hh
@@ -1,35 +1,37 @@
-#ifndef _HAMMING_IMPL_HH
-#define _HAMMING_IMPL_HH
+#pragma once
 
 #include <array>
-#if !(defined(__aarch64__) || defined(_M_ARM64))
-#include <cpuinfo_x86.h>
-#endif
+#include <charconv>
+#include <chrono>
 #include <cstdint>
+#include <iostream>
 #include <limits>
 #include <string>
 #include <vector>
 #ifdef HAMMING_WITH_OPENMP
 #include <omp.h>
 #endif
-#ifdef HAMMING_WITH_NEON
-#include "hamming/distance_neon.hh"
-#endif
-#ifdef HAMMING_WITH_SSE2
-#include "hamming/distance_sse2.hh"
-#endif
-#ifdef HAMMING_WITH_AVX2
-#include "hamming/distance_avx2.hh"
-#endif
-#ifdef HAMMING_WITH_AVX512
-#include "hamming/distance_avx512.hh"
+#ifdef HAMMING_WITH_CUDA
+#include "hamming/distance_cuda.hh"
 #endif
-
 #include "hamming/hamming_impl_types.hh"
 #include "hamming/hamming_types.hh"
 
 namespace hamming {
 
+inline bool cuda_gpu_available() {
+#ifdef HAMMING_WITH_CUDA
+  return distance_cuda_have_device();
+#else
+  return false;
+#endif
+}
+
+typedef int (*distance_func_ptr)(const std::vector<GeneBlock> &,
+                                 const std::vector<GeneBlock> &);
+
+distance_func_ptr get_fastest_supported_distance_func();
+
 std::array<GeneBlock, 256> lookupTable(bool include_x = false);
 
 template <typename DistIntType> DistIntType safe_int_cast(int x) {
@@ -52,21 +54,53 @@ std::vector<SparseData> to_sparse_data(const std::vector<std::string> &data,
 std::vector<std::vector<GeneBlock>>
 to_dense_data(const std::vector<std::string> &data);
 
+std::pair<std::vector<std::string>, std::vector<std::size_t>>
+read_fasta(const std::string &filename, bool remove_duplicates = false,
+           std::size_t n = 0);
+
 std::vector<GeneBlock> from_string(const std::string &str);
 
 template <typename DistIntType>
 std::vector<DistIntType> distances(std::vector<std::string> &data,
-                                   bool include_x, bool clear_input_data) {
+                                   bool include_x, bool clear_input_data,
+                                   bool use_gpu) {
   std::vector<DistIntType> result((data.size() - 1) * data.size() / 2, 0);
+  auto start_time = std::chrono::high_resolution_clock::now();
+  auto print_timing = [&start_time](const std::string &event,
+                                    bool final = false) {
+    std::cout << "# hammingdist :: ..." << event << " completed in "
+              << std::chrono::duration_cast<std::chrono::milliseconds>(
+                     std::chrono::high_resolution_clock::now() - start_time)
+                     .count()
+              << " ms.";
+    if (!final) {
+      std::cout << "..";
+    }
+    std::cout << std::endl;
+    start_time = std::chrono::high_resolution_clock::now();
+  };
+
+#ifndef HAMMING_WITH_CUDA
+  if (use_gpu) {
+    throw std::runtime_error("hammingdist was not compiled with GPU support, "
+                             "please set use_gpu=False");
+  }
+#endif
   auto sparse = to_sparse_data(data, include_x);
   std::size_t nsamples{data.size()};
   std::size_t sample_length{data[0].size()};
 
+  if (include_x && use_gpu) {
+    throw std::runtime_error("use_gpu=True cannot be used if include_x=True, "
+                             "please set use_gpu=False");
+  }
+
   // if X is included, we have to use the sparse distance function
   bool use_sparse = include_x;
   // otherwise, use heuristic to choose distance function: if < 0.5% of values
-  // differ from reference genome, use sparse distance function
-  if (!include_x) {
+  // differ from reference genome, and we're using the CPU, use sparse distance
+  // function
+  if (!include_x && !use_gpu) {
     constexpr double sparse_threshold{0.005};
     std::size_t n_diff{0};
     for (const auto &s : sparse) {
@@ -77,11 +111,14 @@ std::vector<DistIntType> distances(std::vector<std::string> &data,
     use_sparse = frac_diff < sparse_threshold;
   }
   if (use_sparse) {
+    std::cout << "# hammingdist :: Using CPU with sparse distance function..."
+              << std::endl;
     if (clear_input_data) {
       data.clear();
     }
+    print_timing("pre-processing");
 #ifdef HAMMING_WITH_OPENMP
-#pragma omp parallel for
+#pragma omp parallel for default(none) shared(result, sparse, nsamples)
 #endif
     for (std::size_t i = 0; i < nsamples; ++i) {
       std::size_t offset{i * (i - 1) / 2};
@@ -90,52 +127,45 @@ std::vector<DistIntType> distances(std::vector<std::string> &data,
             safe_int_cast<DistIntType>(distance_sparse(sparse[i], sparse[j]));
       }
     }
+    print_timing("distance calculation", true);
     return result;
   }
 
-  // otherwise use fastest supported dense distance function
+  // otherwise use the fastest supported dense distance function
   auto dense = to_dense_data(data);
   if (clear_input_data) {
     data.clear();
   }
-  int (*distance_func)(const std::vector<GeneBlock> &a,
-                       const std::vector<GeneBlock> &b) = distance_cpp;
 
-#if defined(__aarch64__) || defined(_M_ARM64)
-#ifdef HAMMING_WITH_NEON
-  distance_func = distance_neon;
-#endif
-#else
-  const auto features = cpu_features::GetX86Info().features;
-#ifdef HAMMING_WITH_SSE2
-  if (features.sse2) {
-    distance_func = distance_sse2;
-  }
-#endif
-#ifdef HAMMING_WITH_AVX2
-  if (features.avx2) {
-    distance_func = distance_avx2;
-  }
-#endif
-#ifdef HAMMING_WITH_AVX512
-  if (features.avx512bw) {
-    distance_func = distance_avx512;
+#ifdef HAMMING_WITH_CUDA
+  if (use_gpu) {
+    std::cout << "# hammingdist :: Using GPU..." << std::endl;
+    print_timing("pre-processing");
+    if constexpr (sizeof(DistIntType) == 1) {
+      return distances_cuda_8bit(dense);
+    } else if constexpr (sizeof(DistIntType) == 2) {
+      return distances_cuda_16bit(dense);
+    } else {
+      throw std::runtime_error("No GPU implementation available");
+    }
   }
 #endif
-#endif
 
+  auto distance_func{get_fastest_supported_distance_func()};
+  print_timing("pre-processing");
 #ifdef HAMMING_WITH_OPENMP
-#pragma omp parallel for schedule(static, 1)
+#pragma omp parallel for schedule(static, 1) default(none)                     \
+    shared(result, dense, nsamples, distance_func)
 #endif
   for (std::size_t i = 0; i < nsamples; ++i) {
     std::size_t offset{i * (i - 1) / 2};
-    for (std::size_t j = 0; j < i; ++j)
+    for (std::size_t j = 0; j < i; ++j) {
       result[offset + j] =
           safe_int_cast<DistIntType>(distance_func(dense[i], dense[j]));
+    }
   }
+  print_timing("distance calculation", true);
   return result;
 }
 
 } // namespace hamming
-
-#endif
diff --git a/include/hamming/hamming_impl_types.hh b/include/hamming/hamming_impl_types.hh
index d75af12..69d4d44 100644
--- a/include/hamming/hamming_impl_types.hh
+++ b/include/hamming/hamming_impl_types.hh
@@ -1,5 +1,4 @@
-#ifndef _HAMMING_IMPL_TYPES_HH
-#define _HAMMING_IMPL_TYPES_HH
+#pragma once
 
 #include <array>
 #include <cstdint>
@@ -14,5 +13,3 @@ constexpr GeneBlock mask_gene0{0x0f};
 constexpr GeneBlock mask_gene1{0xf0};
 
 } // namespace hamming
-
-#endif
diff --git a/include/hamming/hamming_types.hh b/include/hamming/hamming_types.hh
index 8fdf069..6d6b39c 100644
--- a/include/hamming/hamming_types.hh
+++ b/include/hamming/hamming_types.hh
@@ -1,5 +1,4 @@
-#ifndef _HAMMING_TYPES_HH
-#define _HAMMING_TYPES_HH
+#pragma once
 
 #include <algorithm>
 #include <array>
@@ -15,5 +14,3 @@ using ReferenceDistIntType = uint32_t;
 using DefaultDistIntType = uint8_t;
 
 } // namespace hamming
-
-#endif
diff --git a/include/hamming/hamming_utils.hh b/include/hamming/hamming_utils.hh
new file mode 100644
index 0000000..849a296
--- /dev/null
+++ b/include/hamming/hamming_utils.hh
@@ -0,0 +1,112 @@
+#pragma once
+
+#include <fmt/core.h>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+#ifdef HAMMING_WITH_OPENMP
+#include <omp.h>
+#endif
+
+namespace hamming {
+
+template <typename IntDistType>
+void append_int(std::string &str, IntDistType value) {
+  str.append(fmt::to_string(value));
+}
+
+template <typename DistIntType>
+std::string
+lower_triangular_lines(const std::vector<DistIntType> &partial_distances,
+                       std::size_t i_start, std::size_t i_end, std::size_t i0,
+                       std::size_t j0, std::size_t iN, std::size_t jN) {
+  std::string lines{};
+  std::size_t k{[i_start, i0, j0]() {
+    if (i_start == i0) {
+      // first row, no offset required
+      return std::size_t{0};
+    }
+    // offset for elements from first row and any previous full rows
+    return i0 - j0 + i_start * (i_start - 1) / 2 - (i0 + 1) * i0 / 2;
+  }()};
+  for (std::size_t i = i_start; i < i_end; ++i) {
+    // first row starts at j0, other rows are full and start at 0
+    std::size_t j_first{i == i0 ? j0 : 0};
+    std::size_t j_last{i - 1};
+    char final_char{'\n'};
+    if (i == iN && j_last != jN) {
+      // this is an incomplete final row of partial_distances: no newline at end
+      j_last = jN;
+      final_char = ',';
+    }
+    // do all but last element for this line
+    for (std::size_t j = j_first; j < j_last; ++j) {
+      append_int(lines, partial_distances[k++]);
+      lines.push_back(',');
+    }
+    // do last element for this line with ',' or '\n' as appropriate
+    append_int(lines, partial_distances[k++]);
+    lines.push_back(final_char);
+  }
+  return lines;
+}
+
+static std::size_t row_from_index(std::size_t index) {
+  return static_cast<std::size_t>(
+      floor(sqrt(2.0 * static_cast<double>(index) + 0.5) + 0.5));
+}
+
+static std::size_t col_from_index(std::size_t index, std::size_t row) {
+  return index - row * (row - 1) / 2;
+}
+
+template <typename DistIntType>
+void partial_write_lower_triangular(
+    const std::string &filename,
+    const std::vector<DistIntType> &partial_distances,
+    std::size_t distances_offset, std::size_t n_partial_distances) {
+  if (n_partial_distances == 0) {
+    return;
+  }
+  // infer row/col (i0, j0) of first element of partial_distances
+  std::size_t i0{row_from_index(distances_offset)};
+  std::size_t j0{col_from_index(distances_offset, i0)};
+  // infer row/col (iN, jN) of last element of partial_distances
+  std::size_t last_element_index{distances_offset + n_partial_distances - 1};
+  std::size_t iN{row_from_index(last_element_index)};
+  std::size_t jN{col_from_index(last_element_index, iN)};
+  auto flag = [distances_offset]() {
+    if (distances_offset == 0) {
+      // overwrite any existing data
+      return std::ios_base::out;
+    }
+    // append to existing data
+    return std::ios_base::out | std::ios_base::app | std::ios_base::ate;
+  }();
+  std::ofstream output_file_stream(filename, flag);
+#ifdef HAMMING_WITH_OPENMP
+  constexpr std::size_t max_lines_per_thread{200};
+#pragma omp parallel for schedule(static, 1) ordered default(none)             \
+    shared(partial_distances, output_file_stream, i0, j0, iN, jN)
+  for (std::size_t i_start = i0; i_start <= iN;
+       i_start += max_lines_per_thread) {
+    std::size_t i_end{std::min(i_start + max_lines_per_thread, iN + 1)};
+    auto line = lower_triangular_lines(partial_distances, i_start, i_end, i0,
+                                       j0, iN, jN);
+#pragma omp ordered
+    output_file_stream << line;
+  }
+#else
+  output_file_stream << lower_triangular_lines(partial_distances, i0, iN + 1,
+                                               i0, j0, iN, jN);
+#endif
+}
+
+template <typename DistIntType>
+void write_lower_triangular(const std::string &filename,
+                            const std::vector<DistIntType> &distances) {
+  partial_write_lower_triangular(filename, distances, 0, distances.size());
+}
+
+} // namespace hamming
diff --git a/plots/overview.plt b/plots/overview.plt
new file mode 100644
index 0000000..1134bb9
--- /dev/null
+++ b/plots/overview.plt
@@ -0,0 +1,10 @@
+# to generate: gnuplot overview.plt
+set terminal png enhanced crop
+set output 'overview.png'
+set title 'Hammingdist performance history (million genomes / second)'
+set logscale y
+set ylabel "million genomes / second"
+set xrange [-0.3:4.3]
+set yrange [0.005:200]
+set bmargin 5
+p 'overview.txt' u (1e3/$2):xtic(1) w lp t "" ps 2 pt 7 lc 1 lt 0
diff --git a/plots/overview.png b/plots/overview.png
new file mode 100644
index 0000000..2ed878a
Binary files /dev/null and b/plots/overview.png differ
diff --git a/plots/overview.txt b/plots/overview.txt
new file mode 100644
index 0000000..f36efe5
--- /dev/null
+++ b/plots/overview.txt
@@ -0,0 +1,8 @@
+# time in ns to calculate distance between two genomes of length 32k:
+# note: fairly approximate numbers - see comment for origin of each one
+
+"Python"  54000  # (very rough guess: original dataset was I think ~40k genomes, they said it took many hours to run, so assuming 40k*20k distances in 12 hours we get (12*60*60*1e9)/(40000*20000) ns per distance calc)
+"C++"  3266  # c++ (bench_distance_cpp/32768 benchmark result on hgscomp01)
+"C++\nAVX512"  835  # c++/SIMD-AVX512 (bench_distance_avx512/32768 benchmark result on hgscomp01)
+"C++\nAVX512\nOpenMP\n(52 cores)"  42  # c++/SIMD-AVX512/OpenMP 52-cores (here I just divided above value by 19.7x which was the relative speed-up of a full from_fasta calc using 5000 genomes from the real data with 52 vs 1 core on hgscomp01)
+"C++\nCUDA\n(A100 GPU)" 11 # for 100k genomes from_fasta_to_lower_triangular took 60s, divided by 100000*(100000-1)/2 number of distance calcs, multiplied by 1e9 for s -> ns
diff --git a/plots/speed.plt b/plots/speed.plt
new file mode 100644
index 0000000..121b945
--- /dev/null
+++ b/plots/speed.plt
@@ -0,0 +1,13 @@
+# to generate: gnuplot speed.plt
+set terminal png enhanced crop
+set output 'speed.png'
+set title 'Hammingdist lower triangular output speed'
+set xrange [0:330000]
+set yrange [0:100]
+set ylabel "Million distance matrix coefficients per second"
+set xlabel "Sample count"
+a = 0.5 * 1e-6
+p \
+'speed.txt' u 1:(a*($1**2)/$4) w lp t "from-fasta-to-lower-triangular GPU (A100)", \
+'speed.txt' u 1:(a*($1**2)/$3) w lp t "from-fasta/dump-lower-triangular GPU (A100)", \
+'speed.txt' u 1:(a*($1**2)/$2) w lp t "from-fasta/dump-lower-triangular CPU (52-core Xeon)"
diff --git a/plots/speed.png b/plots/speed.png
new file mode 100644
index 0000000..2e6df31
Binary files /dev/null and b/plots/speed.png differ
diff --git a/plots/speed.txt b/plots/speed.txt
new file mode 100644
index 0000000..4ccaf80
--- /dev/null
+++ b/plots/speed.txt
@@ -0,0 +1,15 @@
+# Time in seconds to go from n 32k-length random genomes to a lower-triangular distances matrix file
+# CPU: 52 core Xeon: from-fasta / dump-lower-triangular
+# GPU: A100: from-fasta / dump-lower-triangular
+# GPU-stream: A100: from-fasta-to-lower-triangular (interleaves I/O and computation)
+# n    CPU        GPU                  GPU-stream
+2000   0.653878   0.5127140621189028   0.13557552301790565
+3000   0.964759   0.7554650978418067   0.23105294106062502
+5000   1.63097    1.3445364689687267   0.46112991601694375
+10000  3.88273    3.1115588791435584   1.34703719697427
+20000  10.779     7.625964550068602    4.461474176961929
+30000  20.7542    13.989239579997957   9.265670311986469
+50000  59.7482    32.871332183945924   20.969019818934612
+100000 185.941    107.00290298496839   81.37370045098942
+200000 735.223    361.97915480996016   271.80367124499753
+300000 1575.75    780.524687159108     610.4017879960593
diff --git a/pyproject.toml b/pyproject.toml
index a1c2119..6bcd607 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "hammingdist"
-version = "0.21.0"
+version = "1.0.0"
 description = "A fast tool to calculate Hamming distances"
 readme = "README.md"
 license = {text = "MIT"}
@@ -37,21 +37,24 @@ Issues = "https://github.com/ssciwr/hammingdist/issues"
 [project.optional-dependencies]
 test = ["pytest", "numpy"]
 
+[tool.scikit-build]
+cmake.verbose = true
+
 [tool.scikit-build.cmake.define]
 BUILD_TESTING = "OFF"
 HAMMING_BUILD_BENCHMARKS = "OFF"
 HAMMING_BUILD_PYTHON = "ON"
 
 [tool.cibuildwheel]
-skip = "*-manylinux_i686"
-test-skip = "pp* *-musllinux*"
+skip = "*-manylinux_i686 *-musllinux*"
+test-skip = "pp*"
 test-extras = "test"
 test-command = "pytest {project}/python/tests -v"
 environment = { BLAS="None", LAPACK="None", ATLAS="None" }
 build-verbosity = 3
 
 [tool.cibuildwheel.linux]
-environment = { CMAKE_ARGS="-DHAMMING_WITH_OPENMP=ON" }
+environment = { CMAKE_ARGS="-DHAMMING_WITH_OPENMP=ON -DHAMMING_WITH_CUDA=ON" }
 
 [[tool.cibuildwheel.overrides]]
 select = "*-macosx_arm64*"
diff --git a/python/hammingdist.cc b/python/hammingdist.cc
index 6c82f19..b78205f 100644
--- a/python/hammingdist.cc
+++ b/python/hammingdist.cc
@@ -59,7 +59,7 @@ PYBIND11_MODULE(hammingdist, m) {
         "(full matrix expected)");
   m.def("from_fasta", &from_fasta<uint8_t>, py::arg("filename"),
         py::arg("include_x") = false, py::arg("remove_duplicates") = false,
-        py::arg("n") = 0,
+        py::arg("n") = 0, py::arg("use_gpu") = false,
         "Creates a dataset by reading from a fasta file (assuming all "
         "sequences have equal length). Maximum value of an element in the "
         "distances matrix: 255. Distances that would have been larger than "
@@ -67,13 +67,26 @@ PYBIND11_MODULE(hammingdist, m) {
         "than this see `from_fasta_large` instead.");
   m.def("from_fasta_large", &from_fasta<uint16_t>, py::arg("filename"),
         py::arg("include_x") = false, py::arg("remove_duplicates") = false,
-        py::arg("n") = 0,
+        py::arg("n") = 0, py::arg("use_gpu") = false,
         "Creates a dataset by reading from a fasta file (assuming all "
         "sequences have equal length). Maximum value of an element in the "
         "distances matrix: 65535");
-  m.def("from_lower_triangular", &from_lower_triangular,
+  m.def("from_fasta_to_lower_triangular", &from_fasta_to_lower_triangular,
+        py::arg("fasta_filename"), py::arg("output_filename"),
+        py::arg("remove_duplicates") = false, py::arg("n") = 0,
+        py::arg("use_gpu") = true,
+        "Construct lower triangular distances matrix output file from the "
+        "fasta file,"
+        "requires an NVIDIA GPU. Maximum value of an element in "
+        "the distances matrix: 65535");
+  m.def("from_lower_triangular", &from_lower_triangular<uint8_t>,
         "Creates a dataset by reading already computed distances from lower "
-        "triangular format");
+        "triangular format. Maximum value of an element in the distances "
+        "matrix: 255.");
+  m.def("from_lower_triangular_large", &from_lower_triangular<uint16_t>,
+        "Creates a dataset by reading already computed distances from lower "
+        "triangular format. Maximum value of an element in the distances "
+        "matrix: 65535.");
   m.def("distance", &distance, py::arg("seq0"), py::arg("seq1"),
         py::arg("include_x") = false,
         "Calculate the distance between seq0 and seq1");
@@ -98,6 +111,8 @@ PYBIND11_MODULE(hammingdist, m) {
       "constructing the distances matrix."
       "For each genome in the input fasta file it gives the index of the "
       "corresponding row in the distances matrix which excludes duplicates");
+  m.def("cuda_gpu_available", &cuda_gpu_available,
+        "True if a GPU that supports CUDA is available");
 }
 
 } // namespace hamming
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 6f2d66c..5215b04 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,15 +1,17 @@
 # Build hamming library
-add_library(hamming STATIC hamming.cc hamming_impl.cc)
+add_library(hamming STATIC hamming.cc hamming_impl.cc hamming_utils.cc)
 target_include_directories(hamming PUBLIC ../include)
+target_include_directories(hamming PRIVATE .)
 target_link_libraries(hamming PUBLIC CpuFeatures::cpu_features)
+target_link_libraries(hamming PUBLIC fmt::fmt-header-only)
 if(HAMMING_WITH_OPENMP)
   find_package(OpenMP REQUIRED)
   target_compile_definitions(hamming PUBLIC HAMMING_WITH_OPENMP)
   target_link_libraries(hamming PUBLIC OpenMP::OpenMP_CXX)
 endif()
 
-# compile optional SIMD code as separate libraries which can be used at runtime
-# if there is CPU support
+# compile optional SIMD/CUDA code as separate libraries which can be used at
+# runtime if there is hardware support
 if(HAMMING_WITH_SSE2)
   target_compile_definitions(hamming PUBLIC HAMMING_WITH_SSE2)
   add_library(distance_sse2 STATIC distance_sse2.cc)
@@ -41,9 +43,27 @@ if(HAMMING_WITH_NEON)
   target_link_libraries(hamming PRIVATE distance_neon)
 endif()
 
+if(HAMMING_WITH_CUDA)
+  target_compile_definitions(hamming PUBLIC HAMMING_WITH_CUDA)
+  add_library(distance_cuda STATIC distance_cuda.cu)
+  target_link_libraries(distance_cuda PUBLIC fmt::fmt-header-only)
+  target_include_directories(distance_cuda PRIVATE .)
+  target_include_directories(distance_cuda PUBLIC ../include)
+  target_link_libraries(hamming PRIVATE distance_cuda)
+  set_target_properties(distance_cuda PROPERTIES CUDA_ARCHITECTURES "all")
+  if(HAMMING_WITH_OPENMP)
+    target_compile_definitions(distance_cuda PUBLIC HAMMING_WITH_OPENMP)
+    target_link_libraries(distance_cuda PUBLIC OpenMP::OpenMP_CXX)
+    # This is required to make nvcc pass the -fopenmp option to gcc:
+    target_compile_options(distance_cuda PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
+                                                 -Xcompiler=-fopenmp>)
+  endif()
+endif()
+
 # Build library benchmarks
 if(HAMMING_BUILD_BENCHMARKS)
-  add_executable(bench bench.cc hamming_bench.cc hamming_impl_bench.cc)
+  add_executable(bench bench.cc hamming_bench.cc hamming_impl_bench.cc
+                       hamming_utils_bench.cc)
   if(HAMMING_WITH_SSE2)
     target_sources(bench PRIVATE distance_sse2_bench.cc)
     target_link_libraries(bench PRIVATE distance_sse2)
@@ -60,13 +80,17 @@ if(HAMMING_BUILD_BENCHMARKS)
     target_sources(bench PRIVATE distance_neon_bench.cc)
     target_link_libraries(bench PRIVATE distance_neon)
   endif()
+  if(HAMMING_WITH_CUDA)
+    target_sources(bench PRIVATE distance_cuda_bench.cc)
+    target_link_libraries(bench PRIVATE distance_cuda)
+  endif()
   target_link_libraries(bench PRIVATE hamming benchmark::benchmark
                                       CpuFeatures::cpu_features)
 endif()
 
 # Build tests
 if(BUILD_TESTING)
-  include(../ext/Catch2/contrib/Catch.cmake)
+  include(../ext/Catch2/extras/Catch.cmake)
   add_executable(tests tests.cc hamming_t.cc hamming_impl_t.cc)
   if(HAMMING_WITH_SSE2)
     target_sources(tests PRIVATE distance_sse2_t.cc)
@@ -84,7 +108,11 @@ if(BUILD_TESTING)
     target_sources(tests PRIVATE distance_neon_t.cc)
     target_link_libraries(tests PRIVATE distance_neon)
   endif()
-  target_link_libraries(tests PRIVATE hamming Catch2::Catch2
+  if(HAMMING_WITH_CUDA)
+    target_sources(tests PRIVATE distance_cuda_t.cc)
+    target_link_libraries(tests PRIVATE distance_cuda)
+  endif()
+  target_link_libraries(tests PRIVATE hamming Catch2::Catch2WithMain
                                       CpuFeatures::cpu_features)
-  catch_discover_tests(tests)
+  catch_discover_tests(tests EXTRA_ARGS "--allow-running-no-tests")
 endif()
diff --git a/src/bench.cc b/src/bench.cc
index f3fede7..c35a30c 100644
--- a/src/bench.cc
+++ b/src/bench.cc
@@ -29,10 +29,11 @@ void randomize_n(std::string &str, std::size_t n, std::mt19937 &gen) {
   }
 }
 
-std::vector<std::string> make_stringlist(int64_t n, std::mt19937 &gen) {
+std::vector<std::string> make_stringlist(int64_t n, int64_t n_samples,
+                                         std::mt19937 &gen) {
   std::vector<std::string> v;
-  v.reserve(n);
-  for (int64_t row = 0; row < n; ++row) {
+  v.reserve(n_samples);
+  for (int64_t row = 0; row < n_samples; ++row) {
     v.push_back(make_string(n, gen));
   }
   return v;
diff --git a/src/bench.hh b/src/bench.hh
index 510ec42..52498cb 100644
--- a/src/bench.hh
+++ b/src/bench.hh
@@ -1,5 +1,4 @@
-#ifndef HAMMING_BENCH_HH
-#define HAMMING_BENCH_HH
+#pragma once
 
 #include <benchmark/benchmark.h>
 #include <random>
@@ -9,11 +8,29 @@
 namespace hamming {
 
 std::string make_string(int64_t n, std::mt19937 &gen, bool include_dash = true);
+
 void randomize_n(std::string &str, std::size_t n, std::mt19937 &gen);
-std::vector<std::string> make_stringlist(int64_t n, std::mt19937 &gen);
+
+std::vector<std::string> make_stringlist(int64_t n, int64_t n_samples,
+                                         std::mt19937 &gen);
+
 void write_fasta(const std::string &filename, const std::string &seq,
                  std::size_t n_seq, std::mt19937 &gen);
 
-} // namespace hamming
+template <typename DistIntType>
+std::vector<DistIntType> make_distances(int64_t n, std::mt19937 &gen) {
+  std::vector<DistIntType> v{};
+  auto n_elements{n * (n - 1) / 2};
+  v.reserve(n_elements);
+  std::uniform_int_distribution<> distrib(
+      0, std::numeric_limits<DistIntType>::max());
+  for (int64_t i = 0; i < n_elements; ++i) {
+    v.push_back(distrib(gen));
+  }
+  return v;
+}
+
+const char *const benchmark_tmp_output_file{"tmp.bench.output"};
+const char *const benchmark_tmp_input_file{"tmp.bench.input"};
 
-#endif
+} // namespace hamming
diff --git a/src/cuda_mem.hh b/src/cuda_mem.hh
new file mode 100644
index 0000000..7f25cf5
--- /dev/null
+++ b/src/cuda_mem.hh
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <stdexcept>
+#include <vector>
+
+template <typename T> T *CheckedCudaMalloc(std::size_t n) {
+  T *a{nullptr};
+  std::size_t sz{sizeof(T) * n};
+  if (cudaError err{cudaMalloc(&a, sz)}; err != cudaSuccess) {
+    throw std::runtime_error(cudaGetErrorString(err));
+  }
+  return a;
+}
+
+template <typename T>
+void CheckedCopyToDevice(T *dest, const std::vector<T> &src) {
+  std::size_t count{sizeof(T) * src.size()}; // count in bytes
+  if (cudaError err{
+          cudaMemcpy(dest, src.data(), count, cudaMemcpyHostToDevice)};
+      err != cudaSuccess) {
+    throw std::runtime_error(cudaGetErrorString(err));
+  }
+}
+
+template <typename T>
+void CheckedCopyToHost(T *dest, const T *src, std::size_t n_elements) {
+  std::size_t count{sizeof(T) * n_elements}; // count in bytes
+  if (cudaError err{cudaMemcpy(dest, src, count, cudaMemcpyDeviceToHost)};
+      err != cudaSuccess) {
+    throw std::runtime_error(cudaGetErrorString(err));
+  }
+}
diff --git a/src/distance_cuda.cu b/src/distance_cuda.cu
new file mode 100644
index 0000000..9ddfcb5
--- /dev/null
+++ b/src/distance_cuda.cu
@@ -0,0 +1,233 @@
+#include "cuda_mem.hh"
+#include "hamming/distance_cuda.hh"
+#include "hamming/hamming_utils.hh"
+#include <chrono>
+#include <cuda/std/limits>
+
+namespace hamming {
+
+template <typename DistIntType>
+__global__ void Dist(DistIntType *partial_distances, const std::uint8_t *genes,
+                     std::uint64_t distances_offset,
+                     unsigned int geneBlocksPerSample) {
+  // Calculates all gridDim.x entries of the partial_distances array.
+  //
+  // The full distances array is a flat nsamples * (nsamples - 1) / 2 element
+  // array that contains the lower-triangular elements of the nSamples x
+  // nSamples partial_distances matrix.
+  //
+  // This kernel is provided with partial_distances, which should have
+  // gridDim.x entries, and which is filled with values corresponding to
+  // entries in the full distances array with an offset of distances_offset.
+  //
+
+  // this array is shared between the threads in this block
+  // and must be large enough to store one int per thread
+  extern __shared__ int s[];
+
+  // index in the partial_distances array where we'll put the result from this
+  // block of threads
+  uint64_t distancesIndex{static_cast<uint64_t>(blockIdx.x)};
+  // index of this value in the full distances array
+  uint64_t trueDistancesIndex{distancesIndex + distances_offset};
+  // infer indices of the two genes corresponding to this distances index
+  uint64_t distancesRowIndex{static_cast<std::size_t>(
+      floor(sqrt(2.0 * static_cast<double>(trueDistancesIndex) + 0.5) + 0.5))};
+  uint64_t distancesColIndex{trueDistancesIndex -
+                             distancesRowIndex * (distancesRowIndex - 1) / 2};
+  uint64_t uint32sPerSample{geneBlocksPerSample / 4};
+  uint64_t geneAIndex{distancesRowIndex * uint32sPerSample};
+  uint64_t geneBIndex{distancesColIndex * uint32sPerSample};
+
+  unsigned int threadIndex{threadIdx.x};
+  // calculate partial sum for each thread and store in shared memory s
+  int r0{0};
+  int r1{0};
+  int r2{0};
+  int r3{0};
+  // NOTE: this cast is only safe if genes is 32-bit aligned AND we each sample
+  // in genes is padded such that the first element of each sample is also
+  // 32-bit aligned!
+  const uint32_t *genes_as_uint32{reinterpret_cast<const uint32_t *>(genes)};
+  uint mask_lower{0x0f0f0f0f};
+  uint mask_upper{0xf0f0f0f0};
+  // NOTE: this loop is also only correct if the length of each sample in genes
+  // is a multiple of 8, which we do by padding the samples with '-'.
+  for (int j = 2 * threadIndex; j < uint32sPerSample; j += 2 * blockDim.x) {
+    auto c0{genes_as_uint32[geneAIndex + j] & genes_as_uint32[geneBIndex + j]};
+    auto c1{genes_as_uint32[geneAIndex + j + 1] &
+            genes_as_uint32[geneBIndex + j + 1]};
+    r0 += __popc(__vseteq4(c0 & mask_lower, 0u));
+    r1 += __popc(__vseteq4(c0 & mask_upper, 0u));
+    r2 += __popc(__vseteq4(c1 & mask_lower, 0u));
+    r3 += __popc(__vseteq4(c1 & mask_upper, 0u));
+  }
+  s[threadIndex] = r0 + r1 + r2 + r3;
+  // synchronise shared memory s between all threads in this block
+  __syncthreads();
+  // sum elements of s using reduction until partial sums are stored in
+  // the first 64 elements of s
+  for (int offset = blockDim.x / 2; offset > 32; offset >>= 1) {
+    if (threadIndex < offset) {
+      s[threadIndex] += s[threadIndex + offset];
+    }
+    __syncthreads();
+  }
+  if (threadIndex < 32) {
+    // one more reduction in each of the 32 threads in this warp
+    int sum{s[threadIndex] + s[threadIndex + 32]};
+    // now sum the values of sum within this warp
+    constexpr unsigned int FULL_MASK{0xffffffff};
+    for (int offset = 16; offset > 0; offset /= 2) {
+      sum += __shfl_down_sync(FULL_MASK, sum, offset);
+    }
+    if (threadIndex == 0) {
+      auto maxDist{cuda::std::numeric_limits<DistIntType>::max()};
+      partial_distances[distancesIndex] = sum > maxDist ? maxDist : sum;
+    }
+  }
+}
+
+template <typename DistIntType>
+std::vector<DistIntType>
+distances_cuda(const std::vector<std::vector<GeneBlock>> &data,
+               const std::string &filename = {}) {
+  std::vector<DistIntType> distances{};
+  std::size_t timing_gpu_ms = 0;
+  std::size_t timing_io_ms = 0;
+  auto timing0{std::chrono::high_resolution_clock::now()};
+  bool output_to_vector{false};
+  if (filename.empty()) {
+    output_to_vector = true;
+  }
+  std::size_t nSamples{data.size()};
+  std::size_t nDistances{nSamples * (nSamples - 1) / 2};
+  std::size_t geneBlocksPerSample{data[0].size()};
+  // 2^31-1 is limit on number of CUDA blocks in x-dim, which corresponds to
+  // 0.5/1GB of distances data for each chunk. For large datasets I/O becomes
+  // the bottleneck and the larger the chunk the faster the I/O tends to be.
+  std::size_t nPartialDistances{std::min(nDistances, 2147483647ul)};
+
+  if (output_to_vector) {
+    // need to store all distances on host
+    distances.resize(nDistances);
+  } else {
+    // only need to store a single block of partial distances on host
+    distances.resize(nPartialDistances);
+  }
+  // allocate memory for genes on device
+  // one gene is 30k chars -> 15k bytes in dense format
+  // so 1 million samples -> 15GB
+  auto *genes{CheckedCudaMalloc<GeneBlock>(nSamples * geneBlocksPerSample)};
+  // copy genes to device
+  for (std::size_t i = 0; i < data.size(); ++i) {
+    CheckedCopyToDevice(genes + i * geneBlocksPerSample, data[i]);
+  }
+
+  // allocate memory for partial distances matrix on device
+  auto *partial_distances{CheckedCudaMalloc<DistIntType>(nPartialDistances)};
+  // keep track of how many distance elements are available to write to disk
+  std::size_t available_distance_elements{0};
+  // keep track of where in the full distances array these elements should go
+  std::size_t distances_offset{0};
+  // GPU timing
+  cudaEvent_t start, stop;
+  cudaEventCreate(&start);
+  cudaEventCreate(&stop);
+  // I/O timing
+  auto timing_io = std::chrono::high_resolution_clock::now();
+  while (distances_offset + available_distance_elements < nDistances) {
+    // use nThreadsPerBlock in x dim of block
+    uint nThreadsPerBlock{128};
+    dim3 threadsPerBlock{nThreadsPerBlock, 1, 1};
+    // use up to nPartialDistances blocks, one block per distance element
+    dim3 numBlocks{static_cast<uint>(std::min(nPartialDistances,
+                                              nDistances - distances_offset)),
+                   1, 1};
+    cudaEventRecord(start);
+    timing_io = std::chrono::high_resolution_clock::now();
+    // launch a kernel with shared memory of size int[nThreadsPerBlock] -
+    // this call returns immediately and the kernel runs asynchronously on the
+    // GPU
+    Dist<<<numBlocks, threadsPerBlock, nThreadsPerBlock * sizeof(int)>>>(
+        partial_distances, genes, distances_offset, geneBlocksPerSample);
+    cudaEventRecord(stop);
+    if (auto err = cudaGetLastError(); err != cudaSuccess) {
+      throw std::runtime_error(cudaGetErrorString(err));
+    }
+
+    if (!output_to_vector) {
+      // write previous kernel's output (if any) to disk using CPU while the new
+      // kernel is running on the GPU to interleave I/O with computation
+      partial_write_lower_triangular(filename, distances, distances_offset,
+                                     available_distance_elements);
+    }
+    timing_io_ms += std::chrono::duration_cast<std::chrono::milliseconds>(
+                        std::chrono::high_resolution_clock::now() - timing_io)
+                        .count();
+    // copy partial_distances from GPU to distances vector on HOST - this call
+    // waits until the kernel has completed before copying the memory.
+    CheckedCopyToHost(distances.data() +
+                          (output_to_vector ? distances_offset : 0),
+                      partial_distances, numBlocks.x);
+    cudaEventSynchronize(stop);
+    float milliseconds = 0;
+    cudaEventElapsedTime(&milliseconds, start, stop);
+    timing_gpu_ms += static_cast<std::size_t>(milliseconds);
+    distances_offset += available_distance_elements;
+    available_distance_elements = numBlocks.x;
+  }
+  if (!output_to_vector) {
+    // write final kernel's output to disk
+    timing_io = std::chrono::high_resolution_clock::now();
+    partial_write_lower_triangular(filename, distances, distances_offset,
+                                   available_distance_elements);
+    timing_io_ms += std::chrono::duration_cast<std::chrono::milliseconds>(
+                        std::chrono::high_resolution_clock::now() - timing_io)
+                        .count();
+    distances.clear();
+  }
+  // free data on gpu
+  cudaFree(genes);
+  cudaFree(partial_distances);
+  std::cout << "# hammingdist :: ...distance calculation completed in "
+            << std::chrono::duration_cast<std::chrono::milliseconds>(
+                   std::chrono::high_resolution_clock::now() - timing0)
+                   .count()
+            << " ms (GPU: " << timing_gpu_ms << " / IO: " << timing_io_ms
+            << ")." << std::endl;
+  return distances;
+}
+
+std::vector<uint8_t>
+distances_cuda_8bit(const std::vector<std::vector<GeneBlock>> &data) {
+  return distances_cuda<uint8_t>(data, {});
+}
+
+std::vector<uint16_t>
+distances_cuda_16bit(const std::vector<std::vector<GeneBlock>> &data) {
+  return distances_cuda<uint16_t>(data, {});
+}
+
+void distances_cuda_to_lower_triangular(
+    const std::vector<std::vector<GeneBlock>> &data,
+    const std::string &filename) {
+  distances_cuda<uint16_t>(data, filename);
+}
+
+int distance_cuda(const std::vector<GeneBlock> &a,
+                  const std::vector<GeneBlock> &b) {
+  // wrapper for testing cuda kernel with existing distance API
+  std::vector<std::vector<GeneBlock>> data{a, b};
+  return distances_cuda<int>(data, {})[0];
+}
+
+bool distance_cuda_have_device() {
+  int nDevices = 0;
+  if (cudaError_t err{cudaGetDeviceCount(&nDevices)}; err != cudaSuccess) {
+    return false;
+  }
+  return nDevices > 0;
+}
+
+} // namespace hamming
diff --git a/src/distance_cuda_bench.cc b/src/distance_cuda_bench.cc
new file mode 100644
index 0000000..ae6ab5c
--- /dev/null
+++ b/src/distance_cuda_bench.cc
@@ -0,0 +1,19 @@
+#include "bench.hh"
+#include "hamming/distance_cuda.hh"
+#include "hamming/hamming.hh"
+
+using namespace hamming;
+
+static void bench_distance_cuda(benchmark::State &state) {
+  std::mt19937 gen(12345);
+  int64_t n{state.range(0)};
+  auto s1{from_string(make_string(n, gen))};
+  auto s2{from_string(make_string(n, gen))};
+  int d{0};
+  for (auto _ : state) {
+    d += distance_cuda(s1, s2);
+  }
+  state.SetComplexityN(n);
+}
+
+BENCHMARK(bench_distance_cuda)->Range(4096, 4194304)->Complexity();
diff --git a/src/distance_cuda_t.cc b/src/distance_cuda_t.cc
new file mode 100644
index 0000000..dbb5b55
--- /dev/null
+++ b/src/distance_cuda_t.cc
@@ -0,0 +1,62 @@
+#include "hamming/distance_cuda.hh"
+#include "tests.hh"
+
+using namespace hamming;
+
+TEST_CASE("distance_cuda() returns all return zero for identical vectors",
+          "[impl][distance][gpu]") {
+  std::mt19937 gen(12345);
+  for (int n :
+       {1,       2,      3,      4,      5,      6,      7,      8,
+        9,       10,     11,     12,     13,     14,     15,     16,
+        17,      18,     19,     20,     31,     32,     33,     63,
+        64,      65,     127,    128,    129,    254,    255,    256,
+        256,     511,    512,    513,    1023,   1024,   1025,   2047,
+        2048,    2049,   4095,   4096,   4097,   8191,   8192,   8193,
+        32767,   32768,  32769,  65535,  65536,  65537,  131071, 131072,
+        131073,  262143, 262144, 262145, 524287, 524288, 524289, 1048575,
+        1048576, 1048577}) {
+    CAPTURE(n);
+    auto g1{make_gene_vector(n, gen)};
+    REQUIRE(distance_cuda(g1, g1) == 0);
+  }
+}
+
+TEST_CASE("distance_cuda() all return n for n A's and n G's",
+          "[impl][distance][gpu]") {
+  for (int n :
+       {1,       2,      3,      4,      5,      6,      7,      8,
+        9,       10,     11,     12,     13,     14,     15,     16,
+        17,      18,     19,     20,     31,     32,     33,     63,
+        64,      65,     127,    128,    129,    254,    255,    256,
+        256,     511,    512,    513,    1023,   1024,   1025,   2047,
+        2048,    2049,   4095,   4096,   4097,   8191,   8192,   8193,
+        32767,   32768,  32769,  65535,  65536,  65537,  131071, 131072,
+        131073,  262143, 262144, 262145, 524287, 524288, 524289, 1048575,
+        1048576, 1048577}) {
+    CAPTURE(n);
+    auto g1 = from_string(std::string(n, 'A'));
+    auto g2 = from_string(std::string(n, 'G'));
+    REQUIRE(distance_cuda(g1, g2) == n);
+  }
+}
+
+TEST_CASE("distance_cuda() returns same as distance_cpp() for random vectors",
+          "[impl][distance][gpu]") {
+  std::mt19937 gen(123451);
+  for (int n :
+       {1,       2,      3,      4,      5,      6,      7,      8,
+        9,       10,     11,     12,     13,     14,     15,     16,
+        17,      18,     19,     20,     31,     32,     33,     63,
+        64,      65,     127,    128,    129,    254,    255,    256,
+        256,     511,    512,    513,    1023,   1024,   1025,   2047,
+        2048,    2049,   4095,   4096,   4097,   8191,   8192,   8193,
+        32767,   32768,  32769,  65535,  65536,  65537,  131071, 131072,
+        131073,  262143, 262144, 262145, 524287, 524288, 524289, 1048575,
+        1048576, 1048577}) {
+    CAPTURE(n);
+    auto g1{make_gene_vector(n, gen)};
+    auto g2{make_gene_vector(n, gen)};
+    REQUIRE(distance_cuda(g1, g2) == distance_cpp(g1, g2));
+  }
+}
diff --git a/src/hamming.cc b/src/hamming.cc
index ff3c6aa..75b26a7 100644
--- a/src/hamming.cc
+++ b/src/hamming.cc
@@ -7,34 +7,45 @@
 #include <fstream>
 #include <iostream>
 #include <limits>
-#include <sstream>
 #include <string>
 #include <unordered_map>
 #include <vector>
 
 namespace hamming {
 
-DataSet<DefaultDistIntType> from_stringlist(std::vector<std::string> &data) {
-  return DataSet<DefaultDistIntType>(data);
+DataSet<DefaultDistIntType> from_stringlist(std::vector<std::string> &data,
+                                            bool include_x, bool use_gpu) {
+  return DataSet<DefaultDistIntType>(data, include_x, false, {}, use_gpu);
 }
 
 DataSet<DefaultDistIntType> from_csv(const std::string &filename) {
   return DataSet<DefaultDistIntType>(filename);
 }
 
-DataSet<DefaultDistIntType> from_lower_triangular(const std::string &filename) {
-  std::vector<DefaultDistIntType> distances;
-  std::ifstream stream(filename);
-  std::string line;
-  while (std::getline(stream, line)) {
-    std::istringstream s(line);
-    std::string d;
-    while (s.good()) {
-      std::getline(s, d, ',');
-      distances.push_back(safe_int_cast<DefaultDistIntType>(std::stoi(d)));
-    }
+void from_fasta_to_lower_triangular(const std::string &input_filename,
+                                    const std::string &output_filename,
+                                    bool remove_duplicates, std::size_t n,
+                                    bool use_gpu) {
+  if (use_gpu) {
+    std::cout << "# hammingdist :: Using GPU..." << std::endl;
+  }
+  auto start_time = std::chrono::high_resolution_clock::now();
+  auto [data, sequence_indices] =
+      read_fasta(input_filename, remove_duplicates, n);
+  auto dense_data = to_dense_data(data);
+  std::cout << "# hammingdist :: ...pre-processing completed in "
+            << std::chrono::duration_cast<std::chrono::milliseconds>(
+                   std::chrono::high_resolution_clock::now() - start_time)
+                   .count()
+            << " ms..." << std::endl;
+#ifdef HAMMING_WITH_CUDA
+  if (use_gpu) {
+    distances_cuda_to_lower_triangular(dense_data, output_filename);
+    return;
   }
-  return DataSet(std::move(distances));
+#endif
+  throw std::runtime_error(
+      "from_fasta_to_lower_triangular is currently only available on GPU");
 }
 
 ReferenceDistIntType distance(const std::string &seq0, const std::string &seq1,
diff --git a/src/hamming_bench.cc b/src/hamming_bench.cc
index 2f4ff28..148b845 100644
--- a/src/hamming_bench.cc
+++ b/src/hamming_bench.cc
@@ -7,13 +7,15 @@
 
 using namespace hamming;
 
+constexpr int64_t sampleLength{30000};
+
 static void bench_from_stringlist(benchmark::State &state) {
 #ifdef HAMMING_WITH_OPENMP
   omp_set_num_threads(1);
 #endif
   std::mt19937 gen(12345);
   int64_t n{state.range(0)};
-  auto v{make_stringlist(n, gen)};
+  auto v{make_stringlist(sampleLength, n, gen)};
   for (auto _ : state) {
     from_stringlist(v);
   }
@@ -24,17 +26,40 @@ static void bench_from_stringlist(benchmark::State &state) {
 static void bench_from_stringlist_omp(benchmark::State &state) {
   omp_set_num_threads(state.range(0));
   std::mt19937 gen(12345);
-  auto v{make_stringlist(8192, gen)};
+  auto v{make_stringlist(sampleLength, 1024, gen)};
   for (auto _ : state) {
     from_stringlist(v);
   }
 }
 #endif
 
+static void bench_from_stringlist_gpu(benchmark::State &state) {
+  std::mt19937 gen(12345);
+  int64_t n{state.range(0)};
+  auto v{make_stringlist(sampleLength, n, gen)};
+  for (auto _ : state) {
+    from_stringlist(v, false, true);
+  }
+  state.SetComplexityN(n);
+}
+
+static void bench_from_fasta_to_lower_triangular_gpu(benchmark::State &state) {
+  std::mt19937 gen(12345);
+  int64_t n{state.range(0)};
+  std::string fasta_file{benchmark_tmp_input_file};
+  std::string lt_file{benchmark_tmp_output_file};
+  auto reference_seq{make_string(sampleLength, gen, true)};
+  write_fasta(fasta_file, reference_seq, state.range(0), gen);
+  for (auto _ : state) {
+    from_fasta_to_lower_triangular(fasta_file, lt_file, false);
+  }
+  state.SetComplexityN(n);
+}
+
 static void bench_fasta_reference_distances(benchmark::State &state) {
   std::mt19937 gen(12345);
-  std::string fasta_file{"fasta.txt"};
-  auto reference_seq{make_string(30000, gen, true)};
+  std::string fasta_file{benchmark_tmp_input_file};
+  auto reference_seq{make_string(sampleLength, gen, true)};
   write_fasta(fasta_file, reference_seq, state.range(0), gen);
   std::vector<ReferenceDistIntType> distances;
   for (auto _ : state) {
@@ -44,7 +69,7 @@ static void bench_fasta_reference_distances(benchmark::State &state) {
 
 BENCHMARK(bench_from_stringlist)
     ->RangeMultiplier(2)
-    ->Range(128, 8192)
+    ->Range(16, 1024)
     ->Complexity();
 #ifdef HAMMING_WITH_OPENMP
 BENCHMARK(bench_from_stringlist_omp)
@@ -55,7 +80,17 @@ BENCHMARK(bench_from_stringlist_omp)
     ->Arg(12)
     ->Arg(24);
 #endif
+#ifdef HAMMING_WITH_CUDA
+BENCHMARK(bench_from_stringlist_gpu)
+    ->RangeMultiplier(2)
+    ->Range(16, 8192)
+    ->Complexity();
+BENCHMARK(bench_from_fasta_to_lower_triangular_gpu)
+    ->RangeMultiplier(2)
+    ->Range(16, 16384)
+    ->Complexity();
+#endif
 BENCHMARK(bench_fasta_reference_distances)
     ->RangeMultiplier(2)
-    ->Range(16, 32384)
+    ->Range(16, 1024)
     ->Complexity();
diff --git a/src/hamming_impl.cc b/src/hamming_impl.cc
index 2ba768f..e2109c9 100644
--- a/src/hamming_impl.cc
+++ b/src/hamming_impl.cc
@@ -1,6 +1,23 @@
 #include "hamming/hamming_impl.hh"
 #include <algorithm>
+#if !(defined(__aarch64__) || defined(_M_ARM64))
+#include <cpuinfo_x86.h>
+#endif
 #include <stdexcept>
+#include <unordered_map>
+#ifdef HAMMING_WITH_SSE2
+#include "hamming/distance_sse2.hh"
+#endif
+#ifdef HAMMING_WITH_AVX2
+#include "hamming/distance_avx2.hh"
+#endif
+#ifdef HAMMING_WITH_AVX512
+#include "hamming/distance_avx512.hh"
+#endif
+#ifdef HAMMING_WITH_NEON
+#include "hamming/distance_neon.hh"
+#endif
+
 namespace hamming {
 
 // bit meaning:
@@ -26,6 +43,40 @@ std::array<GeneBlock, 256> lookupTable(bool include_x) {
   return lookup;
 }
 
+distance_func_ptr get_fastest_supported_distance_func() {
+  std::string simd_str = "no";
+  distance_func_ptr distance_func{distance_cpp};
+#if defined(__aarch64__) || defined(_M_ARM64)
+#ifdef HAMMING_WITH_NEON
+  distance_func = distance_neon;
+  simd_str = "NEON";
+#endif
+#else
+  const auto features = cpu_features::GetX86Info().features;
+#ifdef HAMMING_WITH_SSE2
+  if (features.sse2) {
+    distance_func = distance_sse2;
+    simd_str = "SSE2";
+  }
+#endif
+#ifdef HAMMING_WITH_AVX2
+  if (features.avx2) {
+    distance_func = distance_avx2;
+    simd_str = "AVX2";
+  }
+#endif
+#ifdef HAMMING_WITH_AVX512
+  if (features.avx512bw) {
+    distance_func = distance_avx512;
+    simd_str = "AVX512";
+  }
+#endif
+#endif
+  std::cout << "# hammingdist :: Using CPU with " << simd_str
+            << " SIMD extensions..." << std::endl;
+  return distance_func;
+}
+
 void validate_data(const std::vector<std::string> &data) {
   if (data.empty() || data[0].empty()) {
     throw std::runtime_error("Error: Empty sequence");
@@ -143,11 +194,55 @@ to_dense_data(const std::vector<std::string> &data) {
   return dense;
 }
 
+std::pair<std::vector<std::string>, std::vector<std::size_t>>
+read_fasta(const std::string &filename, bool remove_duplicates, std::size_t n) {
+  std::pair<std::vector<std::string>, std::vector<std::size_t>>
+      data_and_sequence_indices;
+  auto &[data, sequence_indices] = data_and_sequence_indices;
+  data.reserve(n);
+  if (n == 0) {
+    n = std::numeric_limits<std::size_t>::max();
+    data.reserve(65536);
+  }
+  std::unordered_map<std::string, std::size_t> map_seq_to_index;
+  // Initializing the stream
+  std::ifstream stream(filename);
+  std::size_t count = 0;
+  std::size_t count_unique = 0;
+  std::string line;
+  // skip first header
+  std::getline(stream, line);
+  while (count < n && !stream.eof()) {
+    std::string seq{};
+    while (std::getline(stream, line) && line[0] != '>') {
+      seq.append(line);
+    }
+    if (remove_duplicates) {
+      auto result = map_seq_to_index.emplace(std::move(seq), count_unique);
+      if (result.second) {
+        ++count_unique;
+      }
+      sequence_indices.push_back(result.first->second);
+    } else {
+      data.push_back(std::move(seq));
+    }
+    ++count;
+  }
+  if (remove_duplicates) {
+    // copy each unique sequence to the vector of strings
+    data.resize(count_unique);
+    for (auto &key_value_pair : map_seq_to_index) {
+      data[key_value_pair.second] = key_value_pair.first;
+    }
+  }
+  return data_and_sequence_indices;
+}
+
 std::vector<GeneBlock> from_string(const std::string &str) {
   alignas(16) std::vector<GeneBlock> r;
   auto lookup = lookupTable();
   std::size_t n_full_blocks{str.size() / 2};
-  r.reserve(1 + n_full_blocks);
+  r.reserve(1 + n_full_blocks + 3);
   auto iter_str = str.cbegin();
   for (std::size_t i_block = 0; i_block < n_full_blocks; ++i_block) {
     r.push_back(lookup[*iter_str] & mask_gene0);
@@ -160,6 +255,11 @@ std::vector<GeneBlock> from_string(const std::string &str) {
     r.push_back(lookup[*iter_str] & mask_gene0);
     r.back() |= (lookup['-'] & mask_gene1);
   }
+  // pad to ensure 64-bit alignment
+  while (8 * (r.size() / 8) != r.size()) {
+    r.push_back(lookup['-']);
+  }
+
   return r;
 }
 
diff --git a/src/hamming_t.cc b/src/hamming_t.cc
index a88d6d9..0b94e61 100644
--- a/src/hamming_t.cc
+++ b/src/hamming_t.cc
@@ -9,6 +9,12 @@ using namespace hamming;
 static constexpr std::array<char, 4> valid_chars{'A', 'C', 'G', 'T'};
 static constexpr std::array<char, 6> invalid_chars{' ', 'N', '*',
                                                    '?', 'a', '.'};
+static std::vector<bool> valid_use_gpu_values() {
+  if (cuda_gpu_available()) {
+    return {false, true};
+  }
+  return {false};
+}
 
 static int dist1(char c1, char c2) {
   std::vector<std::string> v{std::string{c1}, std::string{c2}};
@@ -16,7 +22,7 @@ static int dist1(char c1, char c2) {
 }
 
 static int dist2(char c1, char c2) {
-  return distance(std::string{c1}, std::string{c2});
+  return static_cast<int>(distance(std::string{c1}, std::string{c2}));
 }
 
 TEST_CASE("distance between two equal valid characters is 0", "[distance]") {
@@ -130,17 +136,26 @@ TEST_CASE("from_fasta single line sequences", "[hamming]") {
   of << ">seq1\n";
   of << "ACGTGTCGTTTCGACGAGTCG\n";
   of.close();
-  for (bool remove_duplicates : {false, true}) {
-    for (bool include_x : {false, true}) {
-      CAPTURE(include_x);
-      CAPTURE(remove_duplicates);
-      for (int n : {0, 2, 3, 8}) {
-        auto d =
-            from_fasta<uint8_t>(tmp_file_name, include_x, remove_duplicates, n);
-        REQUIRE(d[{0, 0}] == 0);
-        REQUIRE(d[{0, 1}] == 2);
-        REQUIRE(d[{1, 0}] == 2);
-        REQUIRE(d[{1, 1}] == 0);
+  for (auto use_gpu : valid_use_gpu_values()) {
+    for (bool remove_duplicates : {false, true}) {
+      for (bool include_x : {false, true}) {
+        CAPTURE(include_x);
+        CAPTURE(remove_duplicates);
+        CAPTURE(use_gpu);
+        if (use_gpu && include_x) {
+          // include_x cannot be used with use_gpu
+          REQUIRE_THROWS(from_fasta<uint8_t>(tmp_file_name, include_x,
+                                             remove_duplicates, 1, use_gpu));
+        } else {
+          for (int n : {0, 2, 3, 8}) {
+            auto d = from_fasta<uint8_t>(tmp_file_name, include_x,
+                                         remove_duplicates, n, use_gpu);
+            REQUIRE(d[{0, 0}] == 0);
+            REQUIRE(d[{0, 1}] == 2);
+            REQUIRE(d[{1, 0}] == 2);
+            REQUIRE(d[{1, 1}] == 0);
+          }
+        }
       }
     }
   }
@@ -166,21 +181,31 @@ TEST_CASE("from_fasta single line sequences with duplicates", "[hamming]") {
   of << "ACGTGTCGTATCGACGTGTCG\n";
   of.close();
   std::vector<std::size_t> sequence_indices{0, 1, 1, 1, 0, 2};
-  for (int n : {0, 6, 22}) {
-    for (bool include_x : {false, true}) {
-      CAPTURE(include_x);
-      auto d = from_fasta<uint8_t>(tmp_file_name, include_x, true, n);
-      REQUIRE(d.nsamples == 3);
-      REQUIRE(d.sequence_indices == sequence_indices);
-      REQUIRE(d[{0, 0}] == 0);
-      REQUIRE(d[{0, 1}] == 2);
-      REQUIRE(d[{0, 2}] == 1);
-      REQUIRE(d[{1, 0}] == 2);
-      REQUIRE(d[{1, 1}] == 0);
-      REQUIRE(d[{1, 2}] == 2);
-      REQUIRE(d[{2, 0}] == 1);
-      REQUIRE(d[{2, 1}] == 2);
-      REQUIRE(d[{2, 2}] == 0);
+  for (auto use_gpu : valid_use_gpu_values()) {
+    for (int n : {0, 6, 22}) {
+      for (bool include_x : {false, true}) {
+        CAPTURE(include_x);
+        CAPTURE(use_gpu);
+        if (use_gpu && include_x) {
+          // include_x cannot be used with use_gpu
+          REQUIRE_THROWS(
+              from_fasta<uint8_t>(tmp_file_name, include_x, true, n, use_gpu));
+        } else {
+          auto d =
+              from_fasta<uint8_t>(tmp_file_name, include_x, true, n, use_gpu);
+          REQUIRE(d.nsamples == 3);
+          REQUIRE(d.sequence_indices == sequence_indices);
+          REQUIRE(d[{0, 0}] == 0);
+          REQUIRE(d[{0, 1}] == 2);
+          REQUIRE(d[{0, 2}] == 1);
+          REQUIRE(d[{1, 0}] == 2);
+          REQUIRE(d[{1, 1}] == 0);
+          REQUIRE(d[{1, 2}] == 2);
+          REQUIRE(d[{2, 0}] == 1);
+          REQUIRE(d[{2, 1}] == 2);
+          REQUIRE(d[{2, 2}] == 0);
+        }
+      }
     }
   }
   std::remove(tmp_file_name);
@@ -200,16 +225,23 @@ TEST_CASE("from_fasta multi-line sequences", "[hamming]") {
   of << "ACGTGTCGTGTCGACGTGTCG\n";
   of << "ACGTGTCGTGTCG\n";
   of.close();
-  for (bool remove_duplicates : {false, true}) {
-    for (bool include_x : {false, true}) {
-      CAPTURE(include_x);
-      for (int n : {0, 2, 3, 8}) {
-        auto d =
-            from_fasta<uint8_t>(tmp_file_name, include_x, remove_duplicates, 2);
-        REQUIRE(d[{0, 0}] == 0);
-        REQUIRE(d[{0, 1}] == 2);
-        REQUIRE(d[{1, 0}] == 2);
-        REQUIRE(d[{1, 1}] == 0);
+  for (auto use_gpu : valid_use_gpu_values()) {
+    for (bool remove_duplicates : {false, true}) {
+      for (bool include_x : {false, true}) {
+        CAPTURE(include_x);
+        CAPTURE(use_gpu);
+        if (use_gpu && include_x) {
+          // include_x cannot be used with use_gpu
+          REQUIRE_THROWS(from_fasta<uint8_t>(tmp_file_name, include_x,
+                                             remove_duplicates, 2, use_gpu));
+        } else {
+          auto d = from_fasta<uint8_t>(tmp_file_name, include_x,
+                                       remove_duplicates, 2, use_gpu);
+          REQUIRE(d[{0, 0}] == 0);
+          REQUIRE(d[{0, 1}] == 2);
+          REQUIRE(d[{1, 0}] == 2);
+          REQUIRE(d[{1, 1}] == 0);
+        }
       }
     }
   }
@@ -242,7 +274,7 @@ TEST_CASE("invalid input data: inconsistent sequence lengths", "[invalid]") {
 
 TEST_CASE("from_csv reproduces correct data", "[hamming]") {
   std::mt19937 gen(12345);
-  std::vector<std::string> data(10);
+  std::vector<std::string> data(107);
   for (auto &d : data)
     d = make_test_string(201, gen);
 
@@ -266,17 +298,21 @@ TEMPLATE_TEST_CASE("distance integer saturates instead of overflowing",
   auto n_max{static_cast<std::size_t>(std::numeric_limits<TestType>::max())};
   std::mt19937 gen(12345);
   std::vector<std::string> data(2);
-  for (auto n : {n_max, n_max + 1, n_max + 99}) {
-    CAPTURE(n);
-    data[0] = std::string(n, 'A');
-    data[1] = std::string(n, 'T');
-    DataSet<TestType> dataSet(data);
-    REQUIRE(dataSet[{0, 1}] == n_max);
-    REQUIRE(dataSet[{1, 0}] == n_max);
+  for (auto use_gpu : valid_use_gpu_values()) {
+    for (auto n : {n_max, n_max + 1, n_max + 99}) {
+      CAPTURE(use_gpu);
+      CAPTURE(n);
+      data[0] = std::string(n, 'A');
+      data[1] = std::string(n, 'T');
+      DataSet<TestType> dataSet(data, false, false, {}, use_gpu);
+      REQUIRE(dataSet[{0, 1}] == n_max);
+      REQUIRE(dataSet[{1, 0}] == n_max);
+    }
   }
 }
 
-TEST_CASE("from_lower_triangular reproduces correct data", "[hamming]") {
+TEMPLATE_TEST_CASE("from_lower_triangular reproduces correct data", "[hamming]",
+                   uint8_t, uint16_t) {
   std::mt19937 gen(12345);
   char tmp_file_name[L_tmpnam];
   REQUIRE(std::tmpnam(tmp_file_name) != nullptr);
@@ -286,11 +322,11 @@ TEST_CASE("from_lower_triangular reproduces correct data", "[hamming]") {
     for (auto &d : data)
       d = make_test_string(24, gen);
 
-    DataSet<uint8_t> ref(data);
+    DataSet<TestType> ref(data);
     REQUIRE(ref.nsamples == n);
     ref.dump_lower_triangular(std::string(tmp_file_name));
 
-    auto restore = from_lower_triangular(std::string(tmp_file_name));
+    auto restore = from_lower_triangular<TestType>(std::string(tmp_file_name));
     REQUIRE(ref.nsamples == restore.nsamples);
     for (std::size_t i = 0; i < ref.nsamples; ++i) {
       for (std::size_t j = 0; j < ref.nsamples; ++j) {
@@ -300,3 +336,101 @@ TEST_CASE("from_lower_triangular reproduces correct data", "[hamming]") {
   }
   std::remove(tmp_file_name);
 }
+
+TEST_CASE("from_stringlist GPU and CPU implementations give consistent results",
+          "[hamming][gpu]") {
+  if (!cuda_gpu_available()) {
+    SKIP("No CUDA gpu available");
+  }
+  std::mt19937 gen(12345);
+  for (bool include_x : {false}) {
+    for (int n : {1, 5, 13, 32, 89, 185, 497, 1092}) {
+      for (int n_samples : {2, 3, 4, 7, 11, 32, 33, 257, 689}) {
+        std::vector<std::string> stringlist;
+        stringlist.reserve(n_samples);
+        for (std::size_t i = 0; i < n_samples; ++i) {
+          stringlist.push_back(make_test_string(n, gen, include_x));
+        }
+        CAPTURE(include_x);
+        auto d_cpu{from_stringlist(stringlist, include_x, false)};
+        auto d_gpu{from_stringlist(stringlist, include_x, true)};
+        for (std::size_t i = 0; i < n_samples; ++i) {
+          for (std::size_t j = 0; j < n_samples; ++j) {
+            REQUIRE(d_cpu[{i, j}] == d_gpu[{i, j}]);
+          }
+        }
+      }
+    }
+  }
+}
+
+TEMPLATE_TEST_CASE(
+    "from_fasta GPU and CPU implementations give consistent results",
+    "[hamming][gpu]", uint8_t, uint16_t) {
+  if (!cuda_gpu_available()) {
+    SKIP("No CUDA gpu available");
+  }
+  char tmp_file_name[L_tmpnam];
+  std::mt19937 gen(12345);
+  REQUIRE(std::tmpnam(tmp_file_name) != nullptr);
+  CAPTURE(tmp_file_name);
+  for (bool remove_duplicates : {false, true}) {
+    for (bool include_x : {false}) {
+      for (int n : {17, 88, 381, 1023}) {
+        for (int n_samples : {2, 3, 4, 7, 11, 127, 128, 255, 256, 257, 703}) {
+          write_test_fasta(tmp_file_name, n, n_samples, gen, include_x);
+          CAPTURE(include_x);
+          auto d_cpu{from_fasta<TestType>(tmp_file_name, include_x,
+                                          remove_duplicates, 0, false)};
+          auto d_gpu{from_fasta<TestType>(tmp_file_name, include_x,
+                                          remove_duplicates, 0, true)};
+          for (std::size_t i = 0; i < n_samples; ++i) {
+            for (std::size_t j = 0; j < n_samples; ++j) {
+              REQUIRE(d_cpu[{i, j}] == d_gpu[{i, j}]);
+            }
+          }
+        }
+      }
+    }
+  }
+  std::remove(tmp_file_name);
+}
+
+TEST_CASE("from_fasta_to_lower_triangular GPU consistent with CPU from_fasta",
+          "[hamming][gpu]") {
+  if (!cuda_gpu_available()) {
+    SKIP("No CUDA gpu available");
+  }
+  std::mt19937 gen(12345);
+  char tmp_fasta_file_name[L_tmpnam];
+  REQUIRE(std::tmpnam(tmp_fasta_file_name) != nullptr);
+  CAPTURE(tmp_fasta_file_name);
+  char tmp_lt_file_name[L_tmpnam];
+  REQUIRE(std::tmpnam(tmp_lt_file_name) != nullptr);
+  CAPTURE(tmp_lt_file_name);
+  for (bool remove_duplicates : {false, true}) {
+    for (int n : {17, 88, 381, 1023}) {
+      for (int n_samples :
+           {2, 3, 4, 7, 11, 127, 128, 255, 256, 257, 703, 1012}) {
+        CAPTURE(remove_duplicates);
+        CAPTURE(n);
+        CAPTURE(n_samples);
+        write_test_fasta(tmp_fasta_file_name, n, n_samples, gen, false);
+        auto d_cpu{from_fasta<uint16_t>(tmp_fasta_file_name, false,
+                                        remove_duplicates, 0, false)};
+        from_fasta_to_lower_triangular(tmp_fasta_file_name, tmp_lt_file_name,
+                                       remove_duplicates, 0, true);
+        auto d_gpu{from_lower_triangular<uint16_t>(tmp_lt_file_name)};
+        for (std::size_t i = 0; i < n_samples; ++i) {
+          for (std::size_t j = 0; j < n_samples; ++j) {
+            CAPTURE(i);
+            CAPTURE(j);
+            REQUIRE(d_cpu[{i, j}] == d_gpu[{i, j}]);
+          }
+        }
+      }
+    }
+  }
+  std::remove(tmp_fasta_file_name);
+  std::remove(tmp_lt_file_name);
+}
diff --git a/src/hamming_utils.cc b/src/hamming_utils.cc
new file mode 100644
index 0000000..312c032
--- /dev/null
+++ b/src/hamming_utils.cc
@@ -0,0 +1,1 @@
+#include "hamming/hamming_utils.hh"
diff --git a/src/hamming_utils_bench.cc b/src/hamming_utils_bench.cc
new file mode 100644
index 0000000..7fae605
--- /dev/null
+++ b/src/hamming_utils_bench.cc
@@ -0,0 +1,78 @@
+#include "bench.hh"
+#include "hamming/hamming.hh"
+#include "hamming/hamming_utils.hh"
+#ifdef HAMMING_WITH_OPENMP
+#include <omp.h>
+#endif
+
+using namespace hamming;
+
+static void bench_distance_cpp(benchmark::State &state) {
+#ifdef HAMMING_WITH_OPENMP
+  omp_set_num_threads(1);
+#endif
+  std::mt19937 gen(12345);
+  int64_t n{state.range(0)};
+  auto s1{from_string(make_string(n, gen))};
+  auto s2{from_string(make_string(n, gen))};
+  int d{0};
+  for (auto _ : state) {
+    d += distance_cpp(s1, s2);
+  }
+  state.SetComplexityN(n);
+}
+
+static void bench_distance_sparse(benchmark::State &state) {
+#ifdef HAMMING_WITH_OPENMP
+  omp_set_num_threads(1);
+#endif
+  std::mt19937 gen(12345);
+  int64_t n{state.range(0)};
+  auto s1{make_string(n, gen, false)};
+  auto s2{s1};
+  // make ~0.5% of s2 elements differ from s1
+  randomize_n(s2, n / 200, gen);
+  auto sparse = to_sparse_data({s1, s2}, false);
+  int d{0};
+  for (auto _ : state) {
+    d += distance_sparse(sparse[0], sparse[1]);
+  }
+  state.SetComplexityN(n);
+}
+
+static void bench_partial_write_lower_triangular(benchmark::State &state) {
+  int64_t n{state.range(0)};
+  std::mt19937 gen(12345);
+  auto v{make_distances<uint16_t>(n, gen)};
+  for (auto _ : state) {
+    partial_write_lower_triangular(benchmark_tmp_output_file, v, 0, v.size());
+  }
+  state.SetComplexityN(n);
+}
+
+#ifdef HAMMING_WITH_OPENMP
+static void bench_partial_write_lower_triangular_omp(benchmark::State &state) {
+  omp_set_num_threads(state.range(0));
+  std::mt19937 gen(12345);
+  auto v{make_distances<uint16_t>(16384, gen)};
+  for (auto _ : state) {
+    partial_write_lower_triangular(benchmark_tmp_output_file, v, 0, v.size());
+  }
+}
+#endif
+
+BENCHMARK(bench_distance_sparse)->Range(4096, 4194304)->Complexity();
+
+BENCHMARK(bench_distance_cpp)->Range(4096, 4194304)->Complexity();
+
+BENCHMARK(bench_partial_write_lower_triangular)->Range(2, 32768)->Complexity();
+
+#ifdef HAMMING_WITH_OPENMP
+BENCHMARK(bench_partial_write_lower_triangular_omp)
+    ->Arg(1)
+    ->Arg(2)
+    ->Arg(4)
+    ->Arg(8)
+    ->Arg(12)
+    ->Arg(24);
+#endif


Test patch:
diff --git a/python/tests/test_hammingdist.py b/python/tests/test_hammingdist.py
index 2f3dfd9..581ba1f 100644
--- a/python/tests/test_hammingdist.py
+++ b/python/tests/test_hammingdist.py
@@ -3,6 +3,8 @@
 import random
 import pytest
 
+gpu_options = [False, True] if hammingdist.cuda_gpu_available() else [False]
+
 
 def write_fasta_file(filename, sequences):
     with open(filename, "w") as f:
@@ -35,7 +37,8 @@ def check_output_sizes(dat, n_in, n_out, tmp_out_file, fasta_sequence_indices=No
 @pytest.mark.parametrize(
     "from_fasta_func", [hammingdist.from_fasta, hammingdist.from_fasta_large]
 )
-def test_from_fasta(from_fasta_func, tmp_path):
+@pytest.mark.parametrize("use_gpu", gpu_options)
+def test_from_fasta(from_fasta_func, use_gpu, tmp_path):
     sequences = [
         "ACGTGTCGTGTCGACGTGTCG",
         "ACGTGTCGTTTCGACGAGTCG",
@@ -48,34 +51,36 @@ def test_from_fasta(from_fasta_func, tmp_path):
     output_file = str(tmp_path / "out.txt")
     write_fasta_file(fasta_file, sequences)
 
-    data = hammingdist.from_fasta(fasta_file)
+    data = from_fasta_func(fasta_file, use_gpu=use_gpu)
     check_output_sizes(data, 6, 6, output_file)
 
-    data = hammingdist.from_fasta(fasta_file, n=5)
+    data = from_fasta_func(fasta_file, n=5, use_gpu=use_gpu)
     check_output_sizes(data, 5, 5, output_file)
 
-    data = hammingdist.from_fasta(fasta_file, include_x=True)
+    data = from_fasta_func(fasta_file, include_x=True, use_gpu=False)
     check_output_sizes(data, 6, 6, output_file)
 
     fasta_sequence_indices = hammingdist.fasta_sequence_indices(fasta_file)
-    data = hammingdist.from_fasta(fasta_file, remove_duplicates=True)
+    data = from_fasta_func(fasta_file, remove_duplicates=True, use_gpu=use_gpu)
     check_output_sizes(data, 6, 4, output_file, fasta_sequence_indices)
 
     fasta_sequence_indices = hammingdist.fasta_sequence_indices(fasta_file)
-    data = hammingdist.from_fasta(fasta_file, remove_duplicates=True, include_x=True)
+    data = from_fasta_func(
+        fasta_file, remove_duplicates=True, include_x=True, use_gpu=False
+    )
     check_output_sizes(data, 6, 4, output_file, fasta_sequence_indices)
 
     fasta_sequence_indices = hammingdist.fasta_sequence_indices(fasta_file, n=2)
-    data = hammingdist.from_fasta(fasta_file, include_x=True, n=2)
+    data = from_fasta_func(fasta_file, include_x=True, n=2, use_gpu=False)
     check_output_sizes(data, 2, 2, output_file, fasta_sequence_indices)
 
     fasta_sequence_indices = hammingdist.fasta_sequence_indices(fasta_file, n=3)
-    data = hammingdist.from_fasta(fasta_file, remove_duplicates=True, n=3)
+    data = from_fasta_func(fasta_file, remove_duplicates=True, n=3, use_gpu=use_gpu)
     check_output_sizes(data, 3, 2, output_file, fasta_sequence_indices)
 
     fasta_sequence_indices = hammingdist.fasta_sequence_indices(fasta_file, n=5)
-    data = hammingdist.from_fasta(
-        fasta_file, remove_duplicates=True, n=5, include_x=True
+    data = from_fasta_func(
+        fasta_file, remove_duplicates=True, n=5, include_x=True, use_gpu=False
     )
     check_output_sizes(data, 5, 3, output_file, fasta_sequence_indices)
 
@@ -125,3 +130,25 @@ def test_distance():
     assert hammingdist.distance("ACGTX", "ACCTX", include_x=True) == 1
     with pytest.raises(RuntimeError):
         hammingdist.distance("ACGT", "ACC")
+
+
+@pytest.mark.skipif(
+    not hammingdist.cuda_gpu_available(),
+    reason="No CUDA GPU available or hammingdist was compiled without CUDA support",
+)
+def test_from_fasta_to_lower_triangular(tmp_path):
+    sequences = [
+        "ACGTGTCGTGTCGACGTGTCGCAGGTGTCGACGTGTCGCAGGTGTCGACGTGTCGCAG",
+        "CCGTGTCGTGTCGACGTGTCGC-GGTGTCGACGTGTCGCAGGTGTCGACGTGTCGCAG",
+        "CAGTGT-GTGTCGACGTGTCGCAGGTGTCGACGTGTCGCAGGTGTCGACGTG--GCAG",
+    ]
+    lower_triangular_dist = [[1], [2, 1]]
+    fasta_file = str(tmp_path / "fasta.txt")
+    output_file = str(tmp_path / "out.txt")
+    write_fasta_file(fasta_file, sequences)
+    hammingdist.from_fasta_to_lower_triangular(fasta_file, output_file)
+    with open(output_file) as f:
+        data = f.read().splitlines()
+    assert len(data) == 2
+    assert np.allclose(np.fromstring(data[0], sep=","), lower_triangular_dist[0])
+    assert np.allclose(np.fromstring(data[1], sep=","), lower_triangular_dist[1])
diff --git a/src/tests.cc b/src/tests.cc
index 66022f4..f933e6e 100644
--- a/src/tests.cc
+++ b/src/tests.cc
@@ -24,4 +24,13 @@ std::vector<GeneBlock> make_gene_vector(int n, std::mt19937 &gen,
   return from_string(make_test_string(n, gen, include_x));
 }
 
+void write_test_fasta(const std::string &filename, int n, std::size_t n_seq,
+                      std::mt19937 &gen, bool include_x) {
+  std::ofstream fs;
+  fs.open(filename);
+  for (std::size_t i = 0; i < n_seq; ++i) {
+    fs << ">seq" << i << "\n" << make_test_string(n, gen, include_x) << "\n";
+  }
+  fs.close();
+}
 } // namespace hamming
diff --git a/src/tests.hh b/src/tests.hh
index 9999496..88284c9 100644
--- a/src/tests.hh
+++ b/src/tests.hh
@@ -1,9 +1,10 @@
-#ifndef HAMMING_TESTS_HH
-#define HAMMING_TESTS_HH
+#pragma once
 
 #include "hamming/hamming.hh"
 #include "hamming/hamming_impl.hh"
-#include <catch2/catch.hpp>
+#include <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers.hpp>
 #include <random>
 #include <vector>
 
@@ -14,6 +15,7 @@ std::string make_test_string(int n, std::mt19937 &gen, bool include_x = false);
 std::vector<GeneBlock> make_gene_vector(int n, std::mt19937 &gen,
                                         bool include_x = false);
 
-} // namespace hamming
+void write_test_fasta(const std::string &filename, int n, std::size_t n_seq,
+                      std::mt19937 &gen, bool include_x = false);
 
-#endif
+} // namespace hamming