Commit 8595f7b8 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Use empty default constructors for shared variables

parent 085ba590

src/libnptm/matrix_inverse_kernel.cu

+12 −49
Changes for src/libnptm/matrix_inverse_kernel.cu: 12 added lines, 49 removed lines.
Original line number Diff line number Diff line
@@ -32,15 +32,19 @@ struct LaunchConfig { 
};



// Global device constants

// #if defined(USE_FP32) || defined(USE_FP64)

//     // For native Float and Double

//     __attribute__((device)) dd_complex dd_zero = { 0.0, 0.0 };

//     __attribute__((device)) dd_complex dd_one  = { 1.0, 0.0 };

// #else

//     // For Double-Double (uses the original struct initialization)

#if defined(USE_FP32) || defined(USE_FP64)

    // For native Float and Double

    __attribute__((device)) const dd_complex dd_zero = { 0.0, 0.0 };

    __attribute__((device)) const dd_complex dd_one  = { 1.0, 0.0 };

#else

    // For Double-Double (uses the original struct initialization)

    // __attribute__((device)) dd_complex dd_zero = { {0.0, 0.0}, {0.0, 0.0} };

    // __attribute__((device)) dd_complex dd_one  = { {1.0, 0.0}, {0.0, 0.0} };

// #endif

    __attribute__((device)) const double DD_ZERO[4] = { 0.0, 0.0, 0.0, 0.0 };

    __attribute__((device)) const double DD_ONE[4] = { 1.0, 0.0, 0.0, 0.0 };

    #define dd_zero ((const dd_complex&)DD_ZERO[0])

    #define dd_one ((const dd_complex&)DD_ONE[0])

#endif





__global__ void lu_decompose_kernel(dd_complex* A, long long* piv, long long n) {
@@ -58,19 +62,10 @@ __global__ void lu_decompose_kernel(dd_complex* A, long long* piv, long long n) 
        dd_real max_val = dd_abs_squared(A[k * n + k]);

        for (long long i = k + 1; i < n; ++i) {

            dd_real val = dd_abs_squared(A[i * n + k]);

	    #if USE_FP32 || USE_FP64

             // Native direct comparison for Float and Double

            if (val > max_val) {

                max_val = val;

                max_idx = i;

            }

             #else

              // Comparison for Double-Double

             if (val.hi > max_val.hi || (val.hi == max_val.hi && val.lo > max_val.lo)) {

             max_val = val;

              max_idx = i;

             }

            #endif

        }

        // Swap lines if necessary

        if (max_idx != k) {
@@ -129,19 +124,10 @@ __global__ void lu_decompose_kernel(dd_complex* A, long long* piv, long long n) 
// 1. Find the pivot and swap the rows

__global__ void lu_pivot_kernel_parallel(dd_complex* A, long long* piv, long long k, long long n) {

    // 1. SHARED MEMORY FOR MAXIMUM REDUCTION

    // __shared__ dd_real s_max_val[BLOCK_SIZE];

    #if defined(USE_FP32)

    #define BUFF_SIZE (BLOCK_SIZE)

    #elif defined(USE_FP64)

    #define BUFF_SIZE (2 * BLOCK_SIZE)

    #else

    #define BUFF_SIZE (4 * BLOCK_SIZE)

    #endif

    __shared__ float b_max_val[BUFF_SIZE];

    __shared__ dd_real s_max_val[BLOCK_SIZE];

    __shared__ long long s_max_idx[BLOCK_SIZE];

    __shared__ long long final_max_idx;

    int tid = threadIdx.x;

    dd_real *s_max_val = (dd_real *)b_max_val;

    // Local initialization for each thread

    // Assume 0 as the initial minimum for abs_squared

#if defined(USE_FP32) || defined(USE_FP64)
@@ -156,19 +142,10 @@ __global__ void lu_pivot_kernel_parallel(dd_complex* A, long long* piv, long lon 
    // 2. PARALLEL SEARCH (Threads divide the lines to be checked)

    for (long long i = k + tid; i < n; i += blockDim.x) {

        dd_real val = dd_abs_squared(A[i * n + k]);

	#if USE_FP32 || USE_FP64

    // Native direct comparison for Float and Double

        if (val > local_max_val) {

            local_max_val = val;

            local_max_idx = i;

        }

     #else

    // Comparison for Double-Double

    if (val.hi > local_max_val.hi || (val.hi == local_max_val.hi && val.lo > local_max_val.lo)) {

        local_max_val = val;

        local_max_idx = i;

         }

        #endif

    }

    // Saving local maxima in shared memory

    s_max_val[tid] = local_max_val;
@@ -178,20 +155,10 @@ __global__ void lu_pivot_kernel_parallel(dd_complex* A, long long* piv, long lon 
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {

        if (tid < s) {

            dd_real val_other = s_max_val[tid + s];

         #if USE_FP32 || USE_FP64

         // Native direct comparison for Float and Double

            if (val_other > s_max_val[tid]) {

                s_max_val[tid] = val_other;

                s_max_idx[tid] = s_max_idx[tid + s];

            }

        #else

       // Comparison for Double-Double 

       if (val_other.hi > s_max_val[tid].hi || 

         (val_other.hi == s_max_val[tid].hi && val_other.lo > s_max_val[tid].lo)) {

          s_max_val[tid] = val_other;

          s_max_idx[tid] = s_max_idx[tid + s];

        }

        #endif

        }

        __syncthreads();

    }
@@ -251,8 +218,6 @@ __global__ void solve_inverse_kernel_optimized( 
    long long n

){

    long long col = blockIdx.x;

    dd_complex dd_zero(0.0, 0.0);

    dd_complex dd_one(1.0, 0.0);

    if (col >= n) return;

    // 1. Move 'y' to the global workspace (each block has its own portion of length n)

    dd_complex* g_vec = workspace + col * n;
@@ -322,7 +287,6 @@ __global__ void solve_system_kernel( 
{

    // A single linear system

    if (blockIdx.x > 0) return;

    dd_complex dd_zero(0.0, 0.0);

    dd_complex* g_vec = workspace;

    // Shared memory for reduction

    extern __shared__ dd_complex s_reduce[];
@@ -432,7 +396,6 @@ __global__ void matrix_mul_dd_kernel( 
    // Calcolo dell'indice di riga e colonna per il thread corrente

    long long row = blockIdx.y * blockDim.y + threadIdx.y;

    long long col = blockIdx.x * blockDim.x + threadIdx.x;

    dd_complex dd_zero(0.0, 0.0);



    if (row < n && col < n) {

        dd_complex acc = dd_zero; // Accumulatore inizializzato a zero (0 + 0i)