Restore SIMD prefix_sum and remove unnecessary simd_float member

Author: Gyuuul2

src/alignment/Fwbw.cpp

@@ -44,26 +44,27 @@ inline void calculate_max4(float& max, float& term1, float& term2, float& term3,
 }
 
 inline simd_float simdf32_prefixsum(simd_float a) {
-//     a = simdf32_add(a, simdi_i2fcast(simdi8_shiftl(simdf_f2icast(a), 4)));
-//     a = simdf32_add(a, simdi_i2fcast(simdi8_shiftl(simdf_f2icast(a), 8)));
-// #ifdef AVX2
-//     a = simdf32_add(a, simdi_i2fcast(simdi8_shiftl(simdf_f2icast(a), 16)));
-// #endif
-//     return a;
-    float buf[8];
-    simdf32_storeu(buf, a);
-
-    buf[1] += buf[0];
-    buf[2] += buf[1];
-    buf[3] += buf[2];
+    a = simdf32_add(a, simdi_i2fcast(simdi8_shiftl(simdf_f2icast(a), 4)));
+    a = simdf32_add(a, simdi_i2fcast(simdi8_shiftl(simdf_f2icast(a), 8)));
 #ifdef AVX2
-    buf[4] += buf[3];
-    buf[5] += buf[4];
-    buf[6] += buf[5];
-    buf[7] += buf[6];
+    a = simdf32_add(a, simdi_i2fcast(simdi8_shiftl(simdf_f2icast(a), 16)));
 #endif
+    return a;
+// Fallback scalar implementation
+//     float buf[8];
+//     simdf32_storeu(buf, a);
+
+//     buf[1] += buf[0];
+//     buf[2] += buf[1];
+//     buf[3] += buf[2];
+// #ifdef AVX2
+//     buf[4] += buf[3];
+//     buf[5] += buf[4];
+//     buf[6] += buf[5];
+//     buf[7] += buf[6];
+// #endif
 
-    return simdf32_loadu(buf);
+//     return simdf32_loadu(buf);
 }
 
 // FwBwAligner Constructor for general case: use profile scoring matrix
@@ -100,8 +101,8 @@ FwBwAligner::FwBwAligner(SubstitutionMatrix &subMat, float gapOpen, float gapExt
         exp_ge_arr[i] = exp((i * gapExtend + gapExtend) / temperature);
     }
     // Gap open and extend
-    exp_go = simdf32_set(static_cast<float>(exp(gapOpen / temperature))); 
-    exp_ge = simdf32_set(static_cast<float>(exp(gapExtend / temperature)));
+    exp_go = (static_cast<float>(exp(gapOpen / temperature))); 
+    exp_ge = (static_cast<float>(exp(gapExtend / temperature)));
     // Blosum matrix
     blosum = malloc_matrix<float>(21, 21);
     for (int i = 0; i < subMat.alphabetSize; ++i) {
@@ -149,8 +150,8 @@ FwBwAligner::FwBwAligner(float gapOpen, float gapExtend, float temperature, floa
         exp_ge_arr[i] = exp((i * gapExtend + gapExtend) / temperature);
     }
     // Gap open and extend
-    exp_go = simdf32_set(static_cast<float>(exp(gapOpen / temperature))); 
-    exp_ge = simdf32_set(static_cast<float>(exp(gapExtend / temperature)));
+    exp_go = (static_cast<float>(exp(gapOpen / temperature))); 
+    exp_ge = (static_cast<float>(exp(gapExtend / temperature)));
 
     if (backtrace != 0) {
         blosum = nullptr;
@@ -333,8 +334,8 @@ void FwBwAligner::resetParams(float newGapOpen, float newGapExtend, float newTem
     gapOpen = newGapOpen;
     gapExtend = newGapExtend;
     temperature = newTemperature;
-    exp_go = simdf32_set(static_cast<float>(exp(gapOpen / temperature))); 
-    exp_ge = simdf32_set(static_cast<float>(exp(gapExtend / temperature)));
+    exp_go = (static_cast<float>(exp(gapOpen / temperature))); 
+    exp_ge = (static_cast<float>(exp(gapExtend / temperature)));
 
     for (size_t i = 0; i < length; ++i) { 
         vj[i] = exp(((length - 1) * gapExtend + gapOpen - i * gapExtend) / temperature);
@@ -420,7 +421,7 @@ void FwBwAligner::forward() {
         std::fill(zInit[i], zInit[i] + rowsCapacity, FLT_MIN_EXP); // rowsCapacity -> tlen
     }  
     max_zm = -std::numeric_limits<float>::max(); 
-    vMax_zm = simdf32_set(max_zm);
+    simd_float vMax_zm = simdf32_set(max_zm);
     P = nullptr; // reset p. do we need this?
     for (size_t b = 0; b < blocks; ++b) {
         size_t start = b * length;
@@ -476,8 +477,8 @@ void FwBwAligner::forward() {
                 simd_float vZmPrev = simdf32_loadu(&zmBlockPrev[j]);
                 simd_float vZf = simdf32_loadu(&zfBlock[j]);
                 simd_float vZfUpdate = simdf32_add(
-                                        simdf32_mul(vZmPrev, exp_go),
-                                        simdf32_mul(vZf, exp_ge)
+                                        simdf32_mul(vZmPrev, simdf32_set(exp_go)),
+                                        simdf32_mul(vZf, simdf32_set(exp_ge))
                                         );
                 vZfUpdate = simdf32_div(vZfUpdate, vZmMaxRowBlock);
                 simdf32_storeu(&zfBlock[j], vZfUpdate);
@@ -491,11 +492,11 @@ void FwBwAligner::forward() {
                 vLastPrefixSum = simdf32_set(vCumsumZm[(VECSIZE_FLOAT - 1)]);
                 simd_float vWj = simdf32_load(&wj[j]);
                 simd_float vExp_ge_arr = simdf32_load(&exp_ge_arr[j]);
-                // simd_float vZeUpdate = simdf32_add(
-                //                         simdf32_div(vCumsumZm, vWj),
-                //                         simdf32_mul(vZeI0, vExp_ge_arr)
-                //                         );
-                simd_float vZeUpdate = simdf32_fmadd(vZeI0, vExp_ge_arr, simdf32_div(vCumsumZm, vWj));
+                simd_float vZeUpdate = simdf32_add(
+                                        simdf32_div(vCumsumZm, vWj),
+                                        simdf32_mul(vZeI0, vExp_ge_arr)
+                                        );
+                // simd_float vZeUpdate = simdf32_fmadd(vZeI0, vExp_ge_arr, simdf32_div(vCumsumZm, vWj));
                 vZeUpdate = simdf32_div(vZeUpdate, vZmMaxRowBlock);
                 simdf32_storeu(&zeBlock[j+1], vZeUpdate);
             }
@@ -693,8 +694,8 @@ void FwBwAligner::backward()  {
                 simd_float vZmPrev = simdf32_loadu(&zmBlockPrev[j]);
                 simd_float vZf = simdf32_loadu(&zfBlock[j]);
                 simd_float vZfUpdate = simdf32_add(
-                                        simdf32_mul(vZmPrev, exp_go),
-                                        simdf32_mul(vZf, exp_ge)
+                                        simdf32_mul(vZmPrev, simdf32_set(exp_go)),
+                                        simdf32_mul(vZf, simdf32_set(exp_ge))
                                         );
                 vZfUpdate = simdf32_div(vZfUpdate, vZmMaxRowBlock);
                 simdf32_storeu(&zfBlock[j], vZfUpdate);
@@ -1186,4 +1187,4 @@ int fwbw(int argc, const char **argv, const Command &command) {
     tdbr.close();
 
     return EXIT_SUCCESS;
-}
+}

src/alignment/Fwbw.h

@@ -115,14 +115,14 @@ class FwBwAligner {
     float** blosum= nullptr; // Profile true
     float* S_prev = nullptr; // backtrace true
     float* S_curr = nullptr; // backtrace true
-    simd_float exp_go;
-    simd_float exp_ge;
+    float exp_go;
+    float exp_ge;
     float max_zm;
     float sum_exp;
     // float maxP;
 
-    simd_float vMax_zm;
-    simd_float vSum_exp;
+    
+    //simd_float vSum_exp;
     size_t colSeqLen_padding;
 
     s_align alignResult;