From 93078bbebce612728dbf788a25db566ad45a816a Mon Sep 17 00:00:00 2001
From: lxy <577086410@qq.com>
Date: Fri, 13 Dec 2024 16:18:00 +0800
Subject: [PATCH] Feat: openmp parallel the pair align loop

Internal ctest changing into directory: /home/lxy/recruitment-2024-autumn/build
Test project /home/lxy/recruitment-2024-autumn/build
    Start 1: tiny case
1/3 Test #1: tiny case ........................   Passed    0.87 sec
    Start 2: small case
2/3 Test #2: small case .......................   Passed    2.41 sec
    Start 3: medium case
3/3 Test #3: medium case ......................   Passed   72.31 sec

100% tests passed, 0 tests failed out of 3

Total Test time (real) =  75.59 sec
---
 CMakeLists.txt       | 27 ++++++++++++----
 build_example.sh     |  5 +++
 env.sh               |  2 ++
 src/SmithWaterman.cc | 74 +++++++++++++++++++++++++++++++-------------
 4 files changed, 81 insertions(+), 27 deletions(-)
 create mode 100644 build_example.sh
 create mode 100644 env.sh

diff --git a/CMakeLists.txt b/CMakeLists.txt
index f104e75..55a0f64 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,10 +20,25 @@ add_test(NAME "tiny case" COMMAND ${CMAKE_BINARY_DIR}/SmithWaterman ${dataset}/t
 add_test(NAME "small case" COMMAND ${CMAKE_BINARY_DIR}/SmithWaterman ${dataset}/small_q.fna ${dataset}/small_t.fna ${dataset}/small_ref.txt)
 add_test(NAME "medium case" COMMAND ${CMAKE_BINARY_DIR}/SmithWaterman ${dataset}/medium_q.fna ${dataset}/medium_t.fna ${dataset}/medium_ref.txt)
 
-if (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC")
-  target_compile_options(SmithWaterman PRIVATE -O3 -acc -Minfo=acc)
-  find_package(OpenACC) # 若成功，执行cmake -B 会输出 "-- Found OpenACC_CXX: -acc"  
-  if(OpenACC_CXX_FOUND)
-    target_link_libraries(SmithWaterman OpenACC::OpenACC_CXX) # 否则无法正确链接acc的一系列库
+option(USE_OPENMP "Enable OpenMP support" OFF)
+option(USE_OPENACC "Enable OpenACC support" OFF)
+
+if(USE_OPENMP)
+  if (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC")
+    target_compile_options(SmithWaterman PRIVATE -O3 -mp -Minfo=mp)
+    find_package(OpenMP)
+    if(OpenMP_CXX_FOUND)
+      target_link_libraries(SmithWaterman OpenMP::OpenMP_CXX)
+    endif()
   endif()
-endif()
\ No newline at end of file
+endif()
+
+if(USE_OPENACC)
+  if (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC")
+    target_compile_options(SmithWaterman PRIVATE -O3 -acc -Minfo=acc)
+    find_package(OpenACC) # 若成功，执行cmake -B 会输出 "-- Found OpenACC_CXX: -acc"  
+    if(OpenACC_CXX_FOUND)
+      target_link_libraries(SmithWaterman OpenACC::OpenACC_CXX) # 否则无法正确链接acc的一系列库
+    endif()
+  endif()
+endif()
diff --git a/build_example.sh b/build_example.sh
new file mode 100644
index 0000000..744741f
--- /dev/null
+++ b/build_example.sh
@@ -0,0 +1,5 @@
+# Select one of the following 2 lines to generate the build files
+cmake -DUSE_OPENMP=ON -DUSE_OPENACC=OFF -B build # for OpenMP
+cmake -DUSE_OPENMP=OFF -DUSE_OPENACC=ON -B build # for OpenACC
+# Build the code
+cmake --build build/
\ No newline at end of file
diff --git a/env.sh b/env.sh
new file mode 100644
index 0000000..39d24ef
--- /dev/null
+++ b/env.sh
@@ -0,0 +1,2 @@
+spack load nvhpc
+spack load cmake@3.27.9
\ No newline at end of file
diff --git a/src/SmithWaterman.cc b/src/SmithWaterman.cc
index 2974617..3ce54a4 100644
--- a/src/SmithWaterman.cc
+++ b/src/SmithWaterman.cc
@@ -18,9 +18,28 @@ SmithWaterman::SmithWaterman(const std::string& query_seq_path,
   assert(query_seqs_size >= 1);
   target_seqs_size = target_seqs.size();
   assert(target_seqs_size >= 1);
+  #pragma acc enter data copyin(this[0:1])
+  for (auto& query_seq : query_seqs) {
+    auto seq_ptr = query_seq.sequence.data();
+    #pragma acc enter data copyin(seq_ptr[:query_seq.sequence.size()])
+  }
+  for (auto& target_seq : target_seqs) {
+    auto seq_ptr = target_seq.sequence.data();
+    #pragma acc enter data copyin(seq_ptr[:target_seq.sequence.size()])
+  }
 }
 
-SmithWaterman::~SmithWaterman() {}
+SmithWaterman::~SmithWaterman() {
+  for (auto& target_seq : target_seqs) {
+    auto seq_ptr = target_seq.sequence.data();
+    #pragma acc exit data delete(seq_ptr[0:target_seq.sequence.size()])
+  }
+  for (auto& query_seq : query_seqs) {
+    auto seq_ptr = query_seq.sequence.data();
+    #pragma acc exit data delete(seq_ptr[0:query_seq.sequence.size()])
+  }
+  #pragma acc exit data delete(this[0:1])
+}
 
 std::vector<size_t> SmithWaterman::solve() {
   // Iterate through the query sequences
@@ -65,30 +84,43 @@ void SmithWaterman::pair_align(FastaSequence& query_seq,
   // Store the highest score in each pairwise-alignment process.
   // Default to 0.
 
+  auto query_sequence_ptr = query_seq.sequence.data();
+  auto target_sequence_ptr = target_seq.sequence.data();
+  auto H_ptr = H.data();
+
   int64_t max_score = 0;
-  // Pairwise-Alignment between the two sequences
-  const int64_t max_x = query_seq_length + target_seq_length;
-  for (int64_t x = 1 + 1; x <= max_x; x++) {
-    int64_t local_max = 0;
-    for(int64_t i = 1; i <= query_seq_length; i++) {
-      int64_t j = x - i;
-      if (j >= 1 && j <= target_seq_length) {
-        int64_t index = pad * i + j;
-        // From the upper element
-        const int64_t up = H[index - pad] + gap_score;
-        // From the left element
-        const int64_t left = H[index - 1] + gap_score;
-        // From the upper-left element
-        const bool match = query_seq.sequence.at(i - 1) == target_seq.sequence.at(j - 1);
-        const int64_t upleft = H[index - pad - 1] + (match ? match_score : mismatch_score);
-        int64_t max = std::max({up, left, upleft, 0l});
-        H[index] = max;
-        local_max = std::max(local_max, max);
+  #pragma acc data present(this[:1], query_sequence_ptr[:query_seq_length], target_sequence_ptr[:target_seq_length])\
+                    create(H_ptr[0:(query_seq_length + 1) * (target_seq_length + 1)]) copy(max_score)
+  {
+    #pragma acc parallel
+    {
+      // Pairwise-Alignment between the two sequences
+      const int64_t max_x = query_seq_length + target_seq_length;
+      #pragma acc loop seq
+      for (int64_t x = 1 + 1; x <= max_x; x++) {
+        int64_t local_max = -999;
+        #pragma acc loop independent reduction(max: local_max)
+        #pragma omp parallel for reduction(max: local_max)
+        for(int64_t i = 1; i <= query_seq_length; i++) {
+          const int64_t j = x - i;
+          if (j >= 1 && j <= target_seq_length) {
+            const int64_t index = pad * i + j;
+            // From the upper element
+            const int64_t up = H_ptr[index - pad] + gap_score;
+            // From the left element
+            const int64_t left = H_ptr[index - 1] + gap_score;
+            // From the upper-left element
+            const bool match  = query_sequence_ptr[i - 1] == target_sequence_ptr[j - 1];
+            const int64_t upleft = H_ptr[index - pad - 1] + (match ? match_score : mismatch_score);
+            const int64_t max = std::max({up, left, upleft, 0l});
+            H_ptr[index] = max;
+            local_max = std::max(local_max, max);
+          }
+        }
+        max_score = std::max(max_score, local_max);
       }
     }
-    max_score = std::max(max_score, local_max);
   }
-
   max_scores.push_back(max_score);
 }