Shared Code

#include <cuda_runtime.h>

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <chrono>
#include <random>
#include <iomanip>

using namespace std;

// ============================================================
// CUDA ERROR CHECK
// // بررسی خطاهای CUDA
// ============================================================
//
// This macro checks every CUDA API call.
// // این ماکرو تمام فراخوانی‌های API مربوط به CUDA را بررسی می‌کند.
//
// If an error occurs, the program prints the CUDA error message
// and terminates.
// // اگر خطایی رخ دهد، پیام خطای CUDA نمایش داده شده و برنامه متوقف می‌شود.
// ============================================================

#define CUDA_CHECK(call)                                      \
do                                                            \
{                                                             \
    cudaError_t error = call;                                 \
                                                              \
    if (error != cudaSuccess)                                 \
    {                                                         \
        cerr << "CUDA ERROR: "                                  \
             << cudaGetErrorString(error)                    \
             << endl;                                         \
        exit(EXIT_FAILURE);                                   \
    }                                                         \
} while (0)


// ============================================================
// NEEDLEMAN-WUNSCH SCORING PARAMETERS
// // پارامترهای امتیازدهی الگوریتم Needleman-Wunsch
// ============================================================
//
// Match     = +2
// // تطابق دو کاراکتر = +2
//
// Mismatch  = -1
// // عدم تطابق دو کاراکتر = -1
//
// Gap       = -2
// // ایجاد شکاف = -2
// ============================================================

const int MATCH_SCORE    = 2;
const int MISMATCH_SCORE = -1;
const int GAP_PENALTY    = -2;


// ============================================================
// DNA GENERATOR
// // تولیدکننده رشته DNA
// ============================================================
//
// Generates reproducible DNA sequences.
// // رشته‌های DNA قابل تکرار تولید می‌کند.
//
// Fixed seeds guarantee that the same input sequences
// are generated every time.
// // استفاده از Seed ثابت باعث می‌شود در هر اجرا
// // دقیقاً همان رشته‌ها تولید شوند.
// ============================================================

string generateDNA(
    int length,
    unsigned int seed)
{
    const char bases[] =
    {
        'A',
        'C',
        'G',
        'T'
    };

    mt19937 generator(seed);

    uniform_int_distribution<int>
        distribution(0, 3);

    string sequence;

    sequence.reserve(length);

    for (int i = 0; i < length; i++)
    {
        sequence +=
            bases[
                distribution(generator)
            ];
    }

    return sequence;
}


// ============================================================
// CPU NEEDLEMAN-WUNSCH - MEMORY OPTIMIZED
// // نسخه CPU الگوریتم Needleman-Wunsch با حافظه بهینه
// ============================================================
//
// The standard Needleman-Wunsch algorithm normally requires
// O(n*m) memory when the complete matrix is stored.
// // نسخه معمولی الگوریتم Needleman-Wunsch برای ذخیره کل ماتریس
// // به O(n*m) حافظه نیاز دارد.
//
// For very large sequences this is impossible.
// // برای رشته‌های بسیار بزرگ این مقدار حافظه عملی نیست.
//
// Therefore, only two rows are stored.
// // بنابراین فقط دو ردیف از ماتریس نگهداری می‌شود.
//
// Time Complexity:
// O(n*m)
// // پیچیدگی زمانی:
// // O(n*m)
//
// Space Complexity:
// O(m)
// // پیچیدگی حافظه:
// // O(m)
// ============================================================

int needlemanWunschCPU(
    const string& seqA,
    const string& seqB)
{
    int n =
        static_cast<int>(
            seqA.length()
        );

    int m =
        static_cast<int>(
            seqB.length()
        );

    // --------------------------------------------------------
    // Allocate only two rows
    // // فقط دو ردیف از ماتریس را ایجاد می‌کنیم
    // --------------------------------------------------------

    vector<int> previousRow(
        m + 1
    );

    vector<int> currentRow(
        m + 1
    );

    // --------------------------------------------------------
    // Initialize first row
    // // مقداردهی اولیه ردیف اول
    // --------------------------------------------------------

    previousRow[0] = 0;

    for (int j = 1; j <= m; j++)
    {
        previousRow[j] =
            j * GAP_PENALTY;
    }

    // --------------------------------------------------------
    // Process all rows
    // // پردازش تمام ردیف‌ها
    // --------------------------------------------------------

    for (int i = 1; i <= n; i++)
    {
        // ----------------------------------------------------
        // Initialize first column
        // // مقداردهی اولین ستون
        // ----------------------------------------------------

        currentRow[0] =
            i * GAP_PENALTY;

        // ----------------------------------------------------
        // Calculate current row
        // // محاسبه ردیف فعلی
        // ----------------------------------------------------

        for (int j = 1; j <= m; j++)
        {
            int diagonal =
                previousRow[j - 1];

            int up =
                previousRow[j];

            int left =
                currentRow[j - 1];

            // ------------------------------------------------
            // Match / Mismatch
            // // محاسبه امتیاز تطابق یا عدم تطابق
            // ------------------------------------------------

            int diagonalScore;

            if (
                seqA[i - 1]
                ==
                seqB[j - 1]
            )
            {
                diagonalScore =
                    diagonal + MATCH_SCORE;
            }
            else
            {
                diagonalScore =
                    diagonal + MISMATCH_SCORE;
            }

            // ------------------------------------------------
            // Gap scores
            // // محاسبه امتیاز Gap
            // ------------------------------------------------

            int upScore =
                up + GAP_PENALTY;

            int leftScore =
                left + GAP_PENALTY;

            // ------------------------------------------------
            // Select maximum score
            // // انتخاب بیشترین امتیاز
            // ------------------------------------------------

            currentRow[j] =
                max(
                    diagonalScore,
                    max(
                        upScore,
                        leftScore
                    )
                );
        }

        // ----------------------------------------------------
        // Swap rows
        // // جابجایی ردیف فعلی و قبلی
        // ----------------------------------------------------

        swap(
            previousRow,
            currentRow
        );
    }

    // --------------------------------------------------------
    // Final alignment score
    // // امتیاز نهایی Alignment
    // --------------------------------------------------------

    return previousRow[m];
}


// ============================================================
// GPU KERNEL
// // هسته GPU
// ============================================================
//
// The algorithm processes one row at a time.
// // الگوریتم ردیف به ردیف پردازش می‌شود.
//
// Threads process independent cells inside a row whenever
// possible.
// // Threadها سلول‌های مستقل یک ردیف را پردازش می‌کنند.
//
// However, because Needleman-Wunsch has dependencies between
// neighboring cells, synchronization is required.
// // اما به دلیل وابستگی سلول‌ها به همسایه‌هایشان، Synchronization
// // مورد نیاز است.
// ============================================================

__global__
void needlemanWunschKernel(
    const char* seqA,
    const char* seqB,
    int* previousRow,
    int* currentRow,
    int n,
    int m,
    int i)
{
    // --------------------------------------------------------
    // Calculate global thread ID
    // // محاسبه شناسه سراسری Thread
    // --------------------------------------------------------

    int j =
        blockIdx.x * blockDim.x
        +
        threadIdx.x
        +
        1;

    // --------------------------------------------------------
    // Ignore threads outside the sequence
    // // Threadهای خارج از محدوده رشته را نادیده می‌گیریم
    // --------------------------------------------------------

    if (j > m)
    {
        return;
    }

    // --------------------------------------------------------
    // The first column is handled separately
    // // ستون اول به صورت جداگانه مقداردهی می‌شود
    // --------------------------------------------------------

    if (j == 1)
    {
        currentRow[0] =
            i * GAP_PENALTY;
    }

    // --------------------------------------------------------
    // Calculate diagonal score
    // // محاسبه امتیاز قطری
    // --------------------------------------------------------

    int diagonalValue =
        previousRow[j - 1];

    int diagonalScore;

    if (
        seqA[i - 1]
        ==
        seqB[j - 1]
    )
    {
        diagonalScore =
            diagonalValue
            +
            MATCH_SCORE;
    }
    else
    {
        diagonalScore =
            diagonalValue
            +
            MISMATCH_SCORE;
    }

    // --------------------------------------------------------
    // Calculate UP score
    // // محاسبه امتیاز از بالا
    // --------------------------------------------------------

    int upScore =
        previousRow[j]
        +
        GAP_PENALTY;

    // --------------------------------------------------------
    // Calculate LEFT score
    // // محاسبه امتیاز از چپ
    // --------------------------------------------------------

    int leftScore =
        currentRow[j - 1]
        +
        GAP_PENALTY;

    // --------------------------------------------------------
    // Select maximum
    // // انتخاب بیشترین مقدار
    // --------------------------------------------------------

    currentRow[j] =
        max(
            diagonalScore,
            max(
                upScore,
                leftScore
            )
        );
}


// ============================================================
// MAIN
// // تابع اصلی برنامه
// ============================================================

int main()
{
    cout << endl;

    cout << "============================================================\n";
    cout << " CUDA Needleman-Wunsch Sequence Alignment\n";
    cout << " Large Sequence Version\n";
    cout << "============================================================\n\n";


    // ========================================================
    // GPU INFORMATION
    // // اطلاعات GPU
    // ========================================================

    int device = 0;

    cudaDeviceProp gpu;

    CUDA_CHECK(
        cudaGetDeviceProperties(
            &gpu,
            device
        )
    );

    cout << "GPU: "
         << gpu.name
         << endl;

    cout << "Compute Capability: "
         << gpu.major
         << "."
         << gpu.minor
         << endl;

    cout << "SM Count: "
         << gpu.multiProcessorCount
         << endl;

    cout << "VRAM: "
         << fixed
         << setprecision(2)
         << (double)gpu.totalGlobalMem
            /
            (1024.0 * 1024.0 * 1024.0)
         << " GB"
         << endl;

    cout << endl;


    // ========================================================
    // GENERATE DNA SEQUENCES
    // // تولید رشته‌های DNA
    // ========================================================
    //
    // Sequence A = 1,000,000 nucleotides
    // // رشته A شامل یک میلیون نوکلئوتید است.
    //
    // Sequence B = 1,200,000 nucleotides
    // // رشته B شامل یک میلیون و دویست هزار نوکلئوتید است.
    // ========================================================

    const int SEQUENCE_A_LENGTH =
        2000;

    const int SEQUENCE_B_LENGTH =
        2400;


    cout << "Generating DNA sequences..."
         << endl;

    auto generationStart =
        chrono::high_resolution_clock::now();


    string sequenceA =
        generateDNA(
            SEQUENCE_A_LENGTH,
            12345
        );


    string sequenceB =
        generateDNA(
            SEQUENCE_B_LENGTH,
            67890
        );


    auto generationEnd =
        chrono::high_resolution_clock::now();


    double generationTime =
        chrono::duration<double, milli>(
            generationEnd
            -
            generationStart
        ).count();


    cout << "DNA generation completed."
         << endl;

    cout << "Generation Time: "
         << fixed
         << setprecision(3)
         << generationTime
         << " ms"
         << endl;

    cout << endl;


    // ========================================================
    // VALIDATE INPUT
    // // بررسی طول رشته‌ها
    // ========================================================

    int n =
        static_cast<int>(
            sequenceA.length()
        );

    int m =
        static_cast<int>(
            sequenceB.length()
        );


    cout << "Sequence A length: "
         << n
         << endl;

    cout << "Sequence B length: "
         << m
         << endl;

    cout << endl;


    if (n != SEQUENCE_A_LENGTH)
    {
        cerr
            << "ERROR: Sequence A length is invalid."
            << endl;

        return 1;
    }


    if (m != SEQUENCE_B_LENGTH)
    {
        cerr
            << "ERROR: Sequence B length is invalid."
            << endl;

        return 1;
    }


    // ========================================================
    // THEORETICAL FULL MATRIX SIZE
    // // اندازه تئوری ماتریس کامل
    // ========================================================
    //
    // A full matrix would contain:
    //
    // (n + 1) * (m + 1)
    //
    // cells.
    //
    // // ماتریس کامل دارای این تعداد سلول خواهد بود:
    // // (n + 1) * (m + 1)
    //
    // We DO NOT allocate this matrix.
    // // این ماتریس را ایجاد نمی‌کنیم.
    // ========================================================

    unsigned long long fullMatrixCells =
        static_cast<unsigned long long>(n + 1)
        *
        static_cast<unsigned long long>(m + 1);


    double theoreticalMemoryGB =
        (
            static_cast<double>(
                fullMatrixCells
            )
            *
            sizeof(int)
        )
        /
        (
            1024.0
            *
            1024.0
            *
            1024.0
            *
            1024.0
        );


    cout << "============================================================\n";
    cout << " MATRIX INFORMATION\n";
    cout << "============================================================\n";


    cout << "Full Matrix Size: "
         << n + 1
         << " x "
         << m + 1
         << endl;


    cout << "Full Matrix Cells: "
         << fullMatrixCells
         << endl;


    cout << "Theoretical Full Matrix Memory: "
         << fixed
         << setprecision(3)
         << theoreticalMemoryGB
         << " TiB"
         << endl;


    cout << endl;


    cout << "IMPORTANT:\n";

    cout << "The full matrix is NOT allocated.\n";

    cout << "Only two rows are stored.\n";


    cout << endl;


    // ========================================================
    // CPU EXECUTION
    // // اجرای الگوریتم روی CPU
    // ========================================================

    cout << "============================================================\n";
    cout << " CPU EXECUTION\n";
    cout << "============================================================\n";


    auto cpuStart =
        chrono::high_resolution_clock::now();


    int cpuFinalScore =
        needlemanWunschCPU(
            sequenceA,
            sequenceB
        );


    auto cpuEnd =
        chrono::high_resolution_clock::now();


    double cpuTime =
        chrono::duration<double, milli>(
            cpuEnd
            -
            cpuStart
        ).count();


    cout << "CPU Final Score: "
         << cpuFinalScore
         << endl;


    cout << "CPU Time: "
         << fixed
         << setprecision(3)
         << cpuTime
         << " ms"
         << endl;


    cout << endl;


    // ========================================================
    // GPU MEMORY ALLOCATION
    // // تخصیص حافظه GPU
    // ========================================================

    char* d_seqA =
        nullptr;

    char* d_seqB =
        nullptr;

    int* d_previousRow =
        nullptr;

    int* d_currentRow =
        nullptr;


    // --------------------------------------------------------
    // Allocate DNA sequence A
    // // تخصیص حافظه برای رشته A
    // --------------------------------------------------------

    CUDA_CHECK(
        cudaMalloc(
            &d_seqA,
            n * sizeof(char)
        )
    );


    // --------------------------------------------------------
    // Allocate DNA sequence B
    // // تخصیص حافظه برای رشته B
    // --------------------------------------------------------

    CUDA_CHECK(
        cudaMalloc(
            &d_seqB,
            m * sizeof(char)
        )
    );


    // --------------------------------------------------------
    // Allocate two DP rows
    // // تخصیص دو ردیف ماتریس Dynamic Programming
    // --------------------------------------------------------

    CUDA_CHECK(
        cudaMalloc(
            &d_previousRow,
            (m + 1) * sizeof(int)
        )
    );


    CUDA_CHECK(
        cudaMalloc(
            &d_currentRow,
            (m + 1) * sizeof(int)
        )
    );


    // ========================================================
    // COPY DNA DATA CPU -> GPU
    // // انتقال DNA از CPU به GPU
    // ========================================================

    CUDA_CHECK(
        cudaMemcpy(
            d_seqA,
            sequenceA.data(),
            n * sizeof(char),
            cudaMemcpyHostToDevice
        )
    );


    CUDA_CHECK(
        cudaMemcpy(
            d_seqB,
            sequenceB.data(),
            m * sizeof(char),
            cudaMemcpyHostToDevice
        )
    );


    // ========================================================
    // INITIALIZE FIRST ROW
    // // مقداردهی ردیف اول
    // ========================================================

    vector<int> firstRow(
        m + 1
    );


    firstRow[0] = 0;


    for (int j = 1; j <= m; j++)
    {
        firstRow[j] =
            j * GAP_PENALTY;
    }


    CUDA_CHECK(
        cudaMemcpy(
            d_previousRow,
            firstRow.data(),
            (m + 1) * sizeof(int),
            cudaMemcpyHostToDevice
        )
    );


    // ========================================================
    // GPU EXECUTION
    // // اجرای الگوریتم روی GPU
    // ========================================================

    cout << "============================================================\n";
    cout << " GPU EXECUTION\n";
    cout << "============================================================\n";


    const int THREADS_PER_BLOCK =
        256;


    int blocks =
        (
            m
            +
            THREADS_PER_BLOCK
            -
            1
        )
        /
        THREADS_PER_BLOCK;


    cout << "Threads per Block: "
         << THREADS_PER_BLOCK
         << endl;


    cout << "Blocks per Row: "
         << blocks
         << endl;


    cout << endl;


    auto gpuStart =
        chrono::high_resolution_clock::now();


    // ========================================================
    // PROCESS ROWS
    // // پردازش ردیف‌های ماتریس
    // ========================================================

    for (int i = 1; i <= n; i++)
    {
        needlemanWunschKernel<<<
            blocks,
            THREADS_PER_BLOCK
        >>>(
            d_seqA,
            d_seqB,
            d_previousRow,
            d_currentRow,
            n,
            m,
            i
        );


        CUDA_CHECK(
            cudaGetLastError()
        );


        // ----------------------------------------------------
        // Synchronize before processing the next row
        // // قبل از پردازش ردیف بعدی Synchronization انجام می‌شود
        // ----------------------------------------------------

        CUDA_CHECK(
            cudaDeviceSynchronize()
        );


        // ----------------------------------------------------
        // Swap previous and current rows
        // // جابجایی ردیف قبلی و ردیف فعلی
        // ----------------------------------------------------

        swap(
            d_previousRow,
            d_currentRow
        );
    }


    auto gpuEnd =
        chrono::high_resolution_clock::now();


    double gpuTime =
        chrono::duration<double, milli>(
            gpuEnd
            -
            gpuStart
        ).count();


    // ========================================================
    // COPY FINAL ROW GPU -> CPU
    // // انتقال ردیف نهایی از GPU به CPU
    // ========================================================

    vector<int> finalRow(
        m + 1
    );


    CUDA_CHECK(
        cudaMemcpy(
            finalRow.data(),
            d_previousRow,
            (m + 1) * sizeof(int),
            cudaMemcpyDeviceToHost
        )
    );


    // ========================================================
    // GPU FINAL SCORE
    // // امتیاز نهایی GPU
    // ========================================================

    int gpuFinalScore =
        finalRow[m];


    // ========================================================
    // VERIFY RESULTS
    // // بررسی صحت نتایج
    // ========================================================

    bool verificationPassed =
        (
            cpuFinalScore
            ==
            gpuFinalScore
        );


    // ========================================================
    // SPEEDUP
    // // محاسبه Speedup
    // ========================================================

    double speedup =
        cpuTime
        /
        gpuTime;


    // ========================================================
    // FINAL RESULTS
    // // نتایج نهایی
    // ========================================================

    cout << endl;

    cout << "============================================================\n";
    cout << " FINAL RESULTS\n";
    cout << "============================================================\n";


    cout << "Sequence A Length : "
         << n
         << endl;


    cout << "Sequence B Length : "
         << m
         << endl;


    cout << endl;


    cout << "Scoring Parameters\n";

    cout << "Match             : "
         << MATCH_SCORE
         << endl;


    cout << "Mismatch          : "
         << MISMATCH_SCORE
         << endl;


    cout << "Gap               : "
         << GAP_PENALTY
         << endl;


    cout << endl;


    cout << "CPU Final Score   : "
         << cpuFinalScore
         << endl;


    cout << "GPU Final Score   : "
         << gpuFinalScore
         << endl;


    cout << endl;


    cout << "CPU Time          : "
         << fixed
         << setprecision(3)
         << cpuTime
         << " ms"
         << endl;


    cout << "GPU Time          : "
         << fixed
         << setprecision(3)
         << gpuTime
         << " ms"
         << endl;


    cout << "Speedup           : "
         << fixed
         << setprecision(3)
         << speedup
         << "x"
         << endl;


    cout << endl;


    cout << "Result Verification: ";


    if (verificationPassed)
    {
        cout << "PASSED";
    }
    else
    {
        cout << "FAILED";
    }


    cout << endl;


    cout << "============================================================\n";


    // ========================================================
    // MEMORY INFORMATION
    // // اطلاعات حافظه
    // ========================================================

    double rowMemoryMB =
        (
            static_cast<double>(
                (m + 1)
                *
                sizeof(int)
            )
        )
        /
        (
            1024.0
            *
            1024.0
        );


    cout << endl;

    cout << "GPU DP Row Memory: "
         << fixed
         << setprecision(3)
         << rowMemoryMB
         << " MB per row"
         << endl;


    cout << "Full Matrix Memory Avoided: "
         << fixed
         << setprecision(3)
         << theoreticalMemoryGB
         << " TiB"
         << endl;


    // ========================================================
    // FREE GPU MEMORY
    // // آزاد کردن حافظه GPU
    // ========================================================

    CUDA_CHECK(
        cudaFree(d_seqA)
    );


    CUDA_CHECK(
        cudaFree(d_seqB)
    );


    CUDA_CHECK(
        cudaFree(d_previousRow)
    );


    CUDA_CHECK(
        cudaFree(d_currentRow)
    );


    // ========================================================
    // RESET CUDA DEVICE
    // // بازنشانی دستگاه CUDA
    // ========================================================

    CUDA_CHECK(
        cudaDeviceReset()
    );


    return 0;
}