Intra-parallel implementation (computeNext & Extend)

Makefile

@@ -16,13 +16,10 @@ CC_FLAGS=-Wall -g
 AR=ar
 AR_FLAGS=-rsc
 
-ifndef BUILD_EXAMPLES 
-BUILD_EXAMPLES=1
-endif
-
-ifndef BUILD_TOOLS 
+BUILD_EXAMPLES=0
 BUILD_TOOLS=1
-endif
+BUILD_WFA_PARALLEL=1
+
 ###############################################################################
 # Configuration rules
 ###############################################################################
@@ -41,10 +38,7 @@ ifeq ($(BUILD_EXAMPLES),1)
     APPS+=examples
 endif
 
-release: CC_FLAGS+=-O3 -march=native -flto
-release: build
-
-all: CC_FLAGS+=-O3 -march=native
+all: CC_FLAGS+=-O3 -march=native#-flto
 all: build
 
 debug: build

tools/align_benchmark/align_benchmark.c

@@ -118,6 +118,7 @@ typedef struct {
   alignment_match_funct_t wfa_match_funct;
   void* wfa_match_funct_arguments;
   uint64_t wfa_max_memory;
+  int wfa_max_threads;
   // Other algorithms parameters
   int bandwidth;
   // Misc
@@ -186,6 +187,7 @@ benchmark_args parameters = {
   .wfa_match_funct = NULL,
   .wfa_match_funct_arguments = NULL,
   .wfa_max_memory = UINT64_MAX,
+  .wfa_max_threads = 1,
   // Other algorithms parameters
   .bandwidth = -1,
   // Misc
@@ -370,6 +372,7 @@ wavefront_aligner_t* align_input_configure_wavefront(
   attributes.plot_params.resolution_points = parameters.plot;
   attributes.system.verbose = parameters.verbose;
   attributes.system.max_memory_abort = parameters.wfa_max_memory;
+  attributes.system.max_num_threads = parameters.wfa_max_threads;
   // Allocate
   return wavefront_aligner_new(&attributes);
 }
@@ -760,6 +763,7 @@ void usage() {
       "              P1 = z-drop                                               \n"
       "              P2 = steps-between-cutoffs                                \n"
       "          --wfa-max-memory <Bytes>                                      \n"
+      "          --wfa-max-threads <INT> (intra-parallelism; default=1)        \n"
       "        [Other Parameters]                                              \n"
       "          --bandwidth <INT>                                             \n"
       "        [Misc]                                                          \n"
@@ -792,6 +796,7 @@ void parse_arguments(int argc,char** argv) {
     { "wfa-heuristic-parameters", required_argument, 0, 1003 },
     { "wfa-custom-match-funct", no_argument, 0, 1004 },
     { "wfa-max-memory", required_argument, 0, 1005 },
+    { "wfa-max-threads", required_argument, 0, 1006 },
     /* Other alignment parameters */
     { "bandwidth", required_argument, 0, 2000 },
     /* Misc */
@@ -977,9 +982,12 @@ void parse_arguments(int argc,char** argv) {
       parameters.wfa_match_funct = match_function;
       parameters.wfa_match_funct_arguments = &match_function_params;
       break;
-    case 1005:
+    case 1005: // --wfa-max-memory
       parameters.wfa_max_memory = atol(optarg);
       break;
+    case 1006: // --wfa-max-threads
+      parameters.wfa_max_threads = atoi(optarg);
+      break;
     /*
      * Other alignment parameters
      */

utils/commons.h

@@ -243,6 +243,12 @@ uint64_t rand_iid(const uint64_t min,const uint64_t max);
 uint32_t nominal_prop_u32(const uint32_t base,const double factor);
 uint64_t nominal_prop_u64(const uint64_t base,const double factor);
 
+/*
+ * Inline
+ */
+#define FORCE_INLINE __attribute__((always_inline)) inline
+#define FORCE_NO_INLINE __attribute__ ((noinline))
+
 /*
  * Vectorize
  */

wavefront/Makefile

@@ -32,7 +32,10 @@ MODULES=wavefront_align \
 SRCS=$(addsuffix .c, $(MODULES))
 OBJS=$(addprefix $(FOLDER_BUILD)/, $(SRCS:.c=.o))
 
-#CC_XFLAGS=-fopt-info-vec-optimized
+ifeq ($(BUILD_WFA_PARALLEL),1) 
+PFLAGS=-DWFA_PARALLEL -fopenmp
+endif
+#CC_XFLAGS+=-fopt-info-vec-optimized
 
 ###############################################################################
 # Rules
@@ -41,6 +44,6 @@ all: $(OBJS)
 
 # General building rule
 $(FOLDER_BUILD)/%.o : %.c
-	$(CC) $(CC_FLAGS) $(CC_XFLAGS) -I$(FOLDER_ROOT) -c $< -o $@
+	$(CC) $(CC_FLAGS) $(PFLAGS) -I$(FOLDER_ROOT) -c $< -o $@
 
 

wavefront/wavefront_align.c

@@ -67,8 +67,8 @@ bool wavefront_align_reached_limits(
     wavefront_aligner_t* const wf_aligner,
     const int score) {
   // Check alignment-score limit
-  if (score >= wf_aligner->alignment_form.max_alignment_score) {
-    wf_aligner->cigar.score = wf_aligner->alignment_form.max_alignment_score;
+  if (score >= wf_aligner->system.max_alignment_score) {
+    wf_aligner->cigar.score = wf_aligner->system.max_alignment_score;
     wf_aligner->align_status.status = WF_STATUS_MAX_SCORE_REACHED;
     return true; // Stop
   }

wavefront/wavefront_aligner.c

@@ -355,7 +355,7 @@ void wavefront_aligner_set_match_funct(
 void wavefront_aligner_set_max_alignment_score(
     wavefront_aligner_t* const wf_aligner,
     const int max_alignment_score) {
-  wf_aligner->alignment_form.max_alignment_score = max_alignment_score;
+  wf_aligner->system.max_alignment_score = max_alignment_score;
 }
 void wavefront_aligner_set_max_memory(
     wavefront_aligner_t* const wf_aligner,

wavefront/wavefront_attributes.c

@@ -44,7 +44,6 @@ wavefront_aligner_attr_t wavefront_aligner_attr_default = {
         .pattern_end_free = 0,
         .text_begin_free = 0,
         .text_end_free = 0,
-        .max_alignment_score = INT_MAX, // Unlimited
     },
     // Custom matching functions
     .match_funct = NULL,           // Use default match-compare function
@@ -91,12 +90,15 @@ wavefront_aligner_attr_t wavefront_aligner_attr_default = {
     },
     // System
     .system = {
-        .check_alignment_correct = false,
+        .max_alignment_score = INT_MAX, // Unlimited
         .probe_interval_global = 2000,
         .probe_interval_compact = 6000,
         .max_memory_compact = -1,  // Automatically set based on memory-mode
         .max_memory_resident = -1, // Automatically set based on memory-mode
         .max_memory_abort = UINT64_MAX, // Unlimited
         .verbose = 0, // Quiet
+        .check_alignment_correct = false,
+        .max_num_threads = 1,           // Single thread by default
+        .min_offsets_per_thread = 1000  // Minimum WF-length to spawn a thread
     },
 };

wavefront/wavefront_attributes.h

@@ -64,8 +64,6 @@ typedef struct {
   int pattern_end_free;    // Allow free-gap at the end of the pattern
   int text_begin_free;     // Allow free-gap at the beginning of the text
   int text_end_free;       // Allow free-gap at the end of the text
-  // Limits
-  int max_alignment_score; // Maximum score allowed before quit
 } alignment_form_t;
 
 /*
@@ -92,8 +90,8 @@ typedef int (*alignment_match_funct_t)(int,int,void*);
  * Alignment system configuration
  */
 typedef struct {
-  // Debug
-  bool check_alignment_correct;  // Verify that the alignment CIGAR output is correct
+  // Limits
+  int max_alignment_score;       // Maximum score allowed before quit
   // Probing intervals
   int probe_interval_global;     // Score-ticks interval to check any limits
   int probe_interval_compact;    // Score-ticks interval to check BT-buffer compacting
@@ -108,8 +106,13 @@ typedef struct {
   //  2 - Report each sequence aligned (brief)
   //  3 - Report each sequence aligned (very verbose)
   int verbose;                   // Verbose (regulates messages during alignment)
+  // Debug
+  bool check_alignment_correct;  // Verify that the alignment CIGAR output is correct
   // Profile
   profiler_timer_t timer;        // Time alignment
+  // OS
+  int max_num_threads;           // Maximum number of threads to use to compute/extend WFs
+  int min_offsets_per_thread;    // Minimum amount of offsets to spawn a thread
 } alignment_system_t;
 
 /*

wavefront/wavefront_compute.c

@@ -415,3 +415,34 @@ void wavefront_compute_trim_ends_set(
   if (wavefront_set->out_i2wavefront) wavefront_compute_trim_ends(wf_aligner,wavefront_set->out_i2wavefront);
   if (wavefront_set->out_d2wavefront) wavefront_compute_trim_ends(wf_aligner,wavefront_set->out_d2wavefront);
 }
+/*
+ * Multithread dispatcher
+ */
+#ifdef WFA_PARALLEL
+int wavefront_compute_num_threads(
+    wavefront_aligner_t* const wf_aligner,
+    const int lo,
+    const int hi) {
+  // Parameters
+  const int max_num_threads = wf_aligner->system.max_num_threads;
+  const int min_offsets_per_thread = wf_aligner->system.min_offsets_per_thread;
+  // Compute minimum work-chunks worth spawning threads
+  const int num_chunks = WAVEFRONT_LENGTH(lo,hi)/min_offsets_per_thread;
+  const int max_workers = MIN(num_chunks,max_num_threads);
+  return MAX(max_workers,1);
+}
+void wavefront_compute_thread_limits(
+    const int thread_id,
+    const int num_theads,
+    const int lo,
+    const int hi,
+    int* const thread_lo,
+    int* const thread_hi) {
+  const int chunk_size = WAVEFRONT_LENGTH(lo,hi)/num_theads;
+  const int t_lo = lo + thread_id*chunk_size;
+  const int t_hi = (thread_id+1 == num_theads) ? hi : t_lo + chunk_size - 1;
+  *thread_lo = t_lo;
+  *thread_hi = t_hi;
+}
+#endif
+

wavefront/wavefront_compute.h

@@ -83,4 +83,23 @@ void wavefront_compute_trim_ends_set(
     wavefront_aligner_t* const wf_aligner,
     wavefront_set_t* const wavefront_set);
 
+/*
+ * Multithread dispatcher
+ */
+#ifdef WFA_PARALLEL
+int wavefront_compute_num_threads(
+    wavefront_aligner_t* const wf_aligner,
+    const int lo,
+    const int hi);
+void wavefront_compute_thread_limits(
+    const int thread_id,
+    const int num_theads,
+    const int lo,
+    const int hi,
+    int* const thread_lo,
+    int* const thread_hi);
+#else
+#define wavefront_compute_num_threads(wf_aligner,lo,hi) 1
+#endif
+
 #endif /* WAVEFRONT_COMPUTE_H_ */

wavefront/wavefront_compute_affine.c

@@ -33,6 +33,10 @@
 #include "wavefront_compute.h"
 #include "wavefront_backtrace_offload.h"
 
+#ifdef WFA_PARALLEL
+#include <omp.h>
+#endif
+
 /*
  * Compute Kernels
  */
@@ -181,8 +185,6 @@ void wavefront_compute_affine_idm_piggyback(
     if (v > pattern_length) max = WAVEFRONT_OFFSET_NULL;
     out_m[k] = max;
   }
-  // Offload backtrace
-  wavefront_backtrace_offload_affine(wf_aligner,wavefront_set,lo,hi);
 }
 /*
  * Compute next wavefront
@@ -201,18 +203,43 @@ void wavefront_compute_affine(
     wavefront_compute_allocate_output_null(wf_aligner,score); // Null s-wavefront
     return;
   }
-  // Set limits
+  // Parameters
+  const bool bt_piggyback = wf_aligner->wf_components.bt_piggyback;
   int hi, lo;
+  // Set limits
   wavefront_compute_limits(wf_aligner,&wavefront_set,&lo,&hi);
   // Allocate wavefronts
   wavefront_compute_allocate_output(wf_aligner,&wavefront_set,score,lo,hi);
   // Init wavefront ends
   wavefront_compute_init_ends(wf_aligner,&wavefront_set,lo,hi);
-  // Compute next wavefront
-  if (wf_aligner->wf_components.bt_piggyback) {
-    wavefront_compute_affine_idm_piggyback(wf_aligner,&wavefront_set,lo,hi);
+  // Multithreading dispatcher
+  const int num_threads = wavefront_compute_num_threads(wf_aligner,lo,hi);
+  if (num_threads == 1) {
+    // Compute next wavefront
+    if (bt_piggyback) {
+      wavefront_compute_affine_idm_piggyback(wf_aligner,&wavefront_set,lo,hi);
+    } else {
+      wavefront_compute_affine_idm(wf_aligner,&wavefront_set,lo,hi);
+    }
   } else {
-    wavefront_compute_affine_idm(wf_aligner,&wavefront_set,lo,hi);
+#ifdef WFA_PARALLEL
+    // Compute next wavefront in parallel
+    #pragma omp parallel num_threads(num_threads)
+    {
+      int t_lo, t_hi;
+      wavefront_compute_thread_limits(
+          omp_get_thread_num(),omp_get_num_threads(),lo,hi,&t_lo,&t_hi);
+      if (bt_piggyback) {
+        wavefront_compute_affine_idm_piggyback(wf_aligner,&wavefront_set,t_lo,t_hi);
+      } else {
+        wavefront_compute_affine_idm(wf_aligner,&wavefront_set,t_lo,t_hi);
+      }
+    }
+#endif
+  }
+  // Offload backtrace (if necessary)
+  if (bt_piggyback) {
+    wavefront_backtrace_offload_affine(wf_aligner,&wavefront_set,lo,hi);
   }
   // Trim wavefront ends
   wavefront_compute_trim_ends_set(wf_aligner,&wavefront_set);