Summation over different determinants that are independently computed using CUDA

Question

Do you have any suggestions for improving the efficiency of the code below?

I believe that better optimization can be implemented in the GPU function cuKer_sum, which is located in the np.cu file.

In this code, I am performing a summation over different determinants that are independently computed for various combinations of g_nxyz and g_temp_nt values. The summation inside the kernel is a bit complex, so modifying it may be challenging. However, I would greatly appreciate any suggestions. Thank you!

Here are the current time stamps:

compiling...
running...

INIT
INIT 1.490000 s
NT_LOOP
-- QPROP 3.030000 s
-- INIT SUM 0.000000 s
-- CUKER 5.390000 s
nt 0     sum 1.2126701710009340E-06,2.0850628227617269E-09
nt 1     sum 1.2511169268362081E-05,2.9725532097971312E-08
-- FINAL SUM 0.000000 s
NT_LOOP 8.420000 s
CUFREE 0.000000 s
TOTAL EXEC: 9.91000 s
CU_DEV_RESET: 0.020000000 s

real    0m10.633s
user    0m5.749s
sys 0m4.706s
task finished!

I run the code using:

#!/usr/bin/env bash

echo compiling...
nvcc -arch=sm_70 np.cu -o np.out

echo running...

time ./np.out

echo task finished!

Here are the CUDA kernels:

1. cuKer_det

__device__ cuDoubleComplex cuKer_det(int gid, int tx,
                                     cuDoubleComplex *d_tqprop, 
                                     int *d_deg_ind,
                                     float *d_deg_c) {


    int b1, b2, b3; 
    int b1p, b2p, b3p; 
    int bt1, bt2, bt3; 
    int bt1p, bt2p, bt3p; 

    cuDoubleComplex d_A[9];
    cuDoubleComplex x1, x2, x3;
    cuDoubleComplex x123, r1x123;
    cuDoubleComplex r1;

    b1 = d_deg_ind[12 * tx];
    b2 = d_deg_ind[12 * tx + 1];
    b3 = d_deg_ind[12 * tx + 2];
    b1p = d_deg_ind[12 * tx + 3];
    b2p = d_deg_ind[12 * tx + 4];
    b3p = d_deg_ind[12 * tx + 5];
    bt1 = d_deg_ind[12 * tx + 6];
    bt2 = d_deg_ind[12 * tx + 7];
    bt3 = d_deg_ind[12 * tx + 8];
    bt1p = d_deg_ind[12 * tx + 9];
    bt2p = d_deg_ind[12 * tx + 10];
    bt3p = d_deg_ind[12 * tx + 11];

    d_A[0*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b1p*4 + bt1p];
    d_A[0*3 + 1] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b2p*4 + bt2p];
    d_A[0*3 + 2] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b3p*4 + bt3p];

    d_A[1*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b2*4*3*4 + bt2*3*4 + b1p*4 + bt1p];
    d_A[1*3 + 1] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b2*4*3*4 + bt2*3*4 + b2p*4 + bt2p];
    d_A[1*3 + 2] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b2*4*3*4 + bt2*3*4 + b3p*4 + bt3p];

    d_A[2*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b3*4*3*4 + bt3*3*4 + b1p*4 + bt1p];
    d_A[2*3 + 1] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b3*4*3*4 + bt3*3*4 + b2p*4 + bt2p];
    d_A[2*3 + 2] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b3*4*3*4 + bt3*3*4 + b3p*4 + bt3p];

    x1 = cuCmul( d_A[0*3 + 0], 
                 cuCsub( cuCmul(d_A[1*3 +1], d_A[2*3 +2]), 
                         cuCmul(d_A[1*3 +2], d_A[2*3 +1]) )  ); 
    x2 = cuCmul( d_A[0*3 + 1], 
                 cuCsub( cuCmul(d_A[1*3 +0], d_A[2*3 +2]), 
                         cuCmul(d_A[1*3 +2], d_A[2*3 +0]) )  );
    x3 = cuCmul( d_A[0*3 + 2], 
                 cuCsub( cuCmul(d_A[1*3 +0], d_A[2*3 +1]), 
                         cuCmul(d_A[1*3 +1], d_A[2*3 +0]) )  ); 

    r1 = make_cuDoubleComplex(d_deg_c[tx], 0.0);
    x123 = cuCadd( cuCsub(x1,x2), x3);
    r1x123 = cuCmul(r1, x123);
       
    return r1x123;

}

2. cuKer_sum

#define SHMEM_SIZE (g_threads)
__global__ void cuKer_sum(cuDoubleComplex *d_tqprop, 
                          cuDoubleComplex *d_sum_nxyz,
                          int *d_deg_ind, int *d_deg_where_d,
                          int *d_deg_len, int *d_where, int *d_start_deg,
                          float *d_deg_c, float *d_deg) {

    int tid = threadIdx.x;
    // int bid = blockIdx.x;
    int gid = blockIdx.x;// * blockDim.x + threadIdx.x;

    
    cuDoubleComplex sumA, temp_sum;
    temp_sum = make_cuDoubleComplex(0.0,0.0);
    cuDoubleComplex x1, x2, x3;

    int start = 0;
    int tx, k;

    __shared__ double sh_sum_re[SHMEM_SIZE];
    __shared__ double sh_sum_im[SHMEM_SIZE];
    

    int min_tx_sum = (tid) * g_count_imp_count_per_thread;
    int max_tx_sum = (tid + 1) * g_count_imp_count_per_thread;

    for (tx = min_tx_sum; tx < max_tx_sum; tx++) {
        if (tx >= g_count_imp_count) {break;}

        x1 = cuKer_det(gid, d_deg_where_d[tx], d_tqprop, d_deg_ind, d_deg_c);

        sumA = make_cuDoubleComplex(0.0,0.0);
        
        start = d_start_deg[tx];
        for (k = start; k < start + d_deg_len[tx]; k++) {
            x2 = cuKer_det(gid, d_where[k], d_tqprop, d_deg_ind, d_deg_c);

            sumA = cuCadd(sumA, cuCmul(x2, make_cuDoubleComplex(d_deg[k], 0.0)) );

        }
        x3 = cuCadd(cuCmul(sumA, x1), x3);

        temp_sum = x3;
    }

    sh_sum_re[tid] = cuCreal(temp_sum);
    sh_sum_im[tid] = cuCimag(temp_sum);
    __syncthreads();




    // Perform block reduction in shared memory
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sh_sum_re[tid] += sh_sum_re[tid + s];
            sh_sum_im[tid] += sh_sum_im[tid + s];
        }
        __syncthreads();
    }


    if (tid == 0) {
        d_sum_nxyz[gid] = make_cuDoubleComplex(sh_sum_re[tid], sh_sum_im[tid]);
    }

}

np.cu (Sorry can not provide the q_nx48_nt144 file. It is 35GB.)

#include "np.cuh"
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <complex.h>
#include <sys/mman.h>
#include <math.h>

#include <sys/time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>

#include <inttypes.h>
#include <string.h>
#include <assert.h>

// DEFINE CONST

    const uint64_t g_count_imp_symm= 140126;
    const uint64_t g_count_deg_index= 38865;
    const uint64_t g_threads= 128;

    const uint64_t g_count_imp_count= 37586;
    const uint64_t g_count_imp_count_per_thread= 294; //ceil(g_count_imp_count/g_threads)

    const uint64_t g_nt = 144;
    const uint64_t g_temp_nt = 2;

    const uint64_t g_nx = 48;
    const uint64_t g_ny = 48;
    const uint64_t g_nz = 48;
    const uint64_t g_nxyz = 110592;
    const uint64_t g_nc = 3;
    const uint64_t g_nd = 4;

    const uint64_t LEN = 48;
    const uint64_t XDIM = (LEN*LEN*LEN);
    const uint64_t ADIM = 3;
    const uint64_t PDIM = 4;
    const uint64_t NRI = 2;
    const uint64_t T = 144;
    const uint64_t LEN_PROP_T = (XDIM * T * ADIM * ADIM * PDIM * PDIM);
    const uint64_t LEN_PROP_TEMP = (XDIM * g_temp_nt * ADIM * ADIM * PDIM * PDIM);


#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) {
   if (code != cudaSuccess) 
   {
      fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}

double big_to_little(double big_endian) {
    union {
        double d;
        uint8_t bytes[8];
    } u;
    u.d = big_endian;
    for (int i = 0; i < 4; i++) {
        uint8_t tmp = u.bytes[i];
        u.bytes[i] = u.bytes[7 - i];
        u.bytes[7 - i] = tmp;
    }
    return u.d;
}


__device__ cuDoubleComplex cuKer_det(int gid, int tx,
                                     cuDoubleComplex *d_tqprop, 
                                     int *d_deg_ind,
                                     float *d_deg_c) {


    int b1, b2, b3; 
    int b1p, b2p, b3p; 
    int bt1, bt2, bt3; 
    int bt1p, bt2p, bt3p; 

    cuDoubleComplex d_A[9];
    cuDoubleComplex x1, x2, x3;
    cuDoubleComplex x123, r1x123;
    cuDoubleComplex r1;

    b1 = d_deg_ind[12 * tx];
    b2 = d_deg_ind[12 * tx + 1];
    b3 = d_deg_ind[12 * tx + 2];
    b1p = d_deg_ind[12 * tx + 3];
    b2p = d_deg_ind[12 * tx + 4];
    b3p = d_deg_ind[12 * tx + 5];
    bt1 = d_deg_ind[12 * tx + 6];
    bt2 = d_deg_ind[12 * tx + 7];
    bt3 = d_deg_ind[12 * tx + 8];
    bt1p = d_deg_ind[12 * tx + 9];
    bt2p = d_deg_ind[12 * tx + 10];
    bt3p = d_deg_ind[12 * tx + 11];

    d_A[0*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b1p*4 + bt1p];
    d_A[0*3 + 1] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b2p*4 + bt2p];
    d_A[0*3 + 2] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b3p*4 + bt3p];

    d_A[1*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b2*4*3*4 + bt2*3*4 + b1p*4 + bt1p];
    d_A[1*3 + 1] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b2*4*3*4 + bt2*3*4 + b2p*4 + bt2p];
    d_A[1*3 + 2] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b2*4*3*4 + bt2*3*4 + b3p*4 + bt3p];

    d_A[2*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b3*4*3*4 + bt3*3*4 + b1p*4 + bt1p];
    d_A[2*3 + 1] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b3*4*3*4 + bt3*3*4 + b2p*4 + bt2p];
    d_A[2*3 + 2] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b3*4*3*4 + bt3*3*4 + b3p*4 + bt3p];

    x1 = cuCmul( d_A[0*3 + 0], 
                 cuCsub( cuCmul(d_A[1*3 +1], d_A[2*3 +2]), 
                         cuCmul(d_A[1*3 +2], d_A[2*3 +1]) )  ); 
    x2 = cuCmul( d_A[0*3 + 1], 
                 cuCsub( cuCmul(d_A[1*3 +0], d_A[2*3 +2]), 
                         cuCmul(d_A[1*3 +2], d_A[2*3 +0]) )  );
    x3 = cuCmul( d_A[0*3 + 2], 
                 cuCsub( cuCmul(d_A[1*3 +0], d_A[2*3 +1]), 
                         cuCmul(d_A[1*3 +1], d_A[2*3 +0]) )  ); 

    r1 = make_cuDoubleComplex(d_deg_c[tx], 0.0);
    x123 = cuCadd( cuCsub(x1,x2), x3);
    r1x123 = cuCmul(r1, x123);
       
    return r1x123;

}

#define SHMEM_SIZE (g_threads)
__global__ void cuKer_sum(cuDoubleComplex *d_tqprop, 
                          cuDoubleComplex *d_sum_nxyz,
                          int *d_deg_ind, int *d_deg_where_d,
                          int *d_deg_len, int *d_where, int *d_start_deg,
                          float *d_deg_c, float *d_deg) {

    int tid = threadIdx.x;
    // int bid = blockIdx.x;
    int gid = blockIdx.x;// * blockDim.x + threadIdx.x;

    
    cuDoubleComplex sumA, temp_sum;
    temp_sum = make_cuDoubleComplex(0.0,0.0);
    cuDoubleComplex x1, x2, x3;

    int start = 0;
    int tx, k;

    __shared__ double sh_sum_re[SHMEM_SIZE];
    __shared__ double sh_sum_im[SHMEM_SIZE];
    

    int min_tx_sum = (tid) * g_count_imp_count_per_thread;
    int max_tx_sum = (tid + 1) * g_count_imp_count_per_thread;

    for (tx = min_tx_sum; tx < max_tx_sum; tx++) {
        if (tx >= g_count_imp_count) {break;}

        x1 = cuKer_det(gid, d_deg_where_d[tx], d_tqprop, d_deg_ind, d_deg_c);

        sumA = make_cuDoubleComplex(0.0,0.0);
        
        start = d_start_deg[tx];
        for (k = start; k < start + d_deg_len[tx]; k++) {
            x2 = cuKer_det(gid, d_where[k], d_tqprop, d_deg_ind, d_deg_c);

            sumA = cuCadd(sumA, cuCmul(x2, make_cuDoubleComplex(d_deg[k], 0.0)) );

        }
        x3 = cuCadd(cuCmul(sumA, x1), x3);

        temp_sum = x3;
    }

    sh_sum_re[tid] = cuCreal(temp_sum);
    sh_sum_im[tid] = cuCimag(temp_sum);
    __syncthreads();




    // Perform block reduction in shared memory
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sh_sum_re[tid] += sh_sum_re[tid + s];
            sh_sum_im[tid] += sh_sum_im[tid + s];
        }
        __syncthreads();
    }


    if (tid == 0) {
        d_sum_nxyz[gid] = make_cuDoubleComplex(sh_sum_re[tid], sh_sum_im[tid]);
    }

}



int main(int argc, char *argv[]) {

    clock_t code_start, code_end;
    code_start = clock();

    clock_t ini_start, ini_end;
    ini_start = clock();

    //START INIT//
        printf("\nINIT\n");

        int d_trial;
        cudaMalloc((void **)&d_trial, sizeof(int));

        //READ INDEX FILES//
        
            int *where;
            float *deg;

            where = (int*)malloc(g_count_imp_symm * sizeof(int));
            deg = (float*)malloc(g_count_imp_symm * sizeof(float));

            f_read_imp_symm(where, deg);

            int *deg_ind;
            float *deg_c;

            deg_ind = (int*)malloc(12 * g_count_deg_index * sizeof(int));
            deg_c = (float*)malloc(1 * g_count_deg_index * sizeof(float));

            f_read_deg_index(deg_ind, deg_c);

            int *deg_where_d; 
            int *deg_len;
            int *start_deg;

            deg_where_d = (int*)malloc(g_count_imp_count * sizeof(int));
            deg_len = (int*)malloc(g_count_imp_count * sizeof(int));
            start_deg = (int *)malloc(g_count_imp_count * sizeof(int));

            f_read_imp_count(deg_where_d, deg_len, start_deg);


            // device variable independent of lattice index
            int *d_deg_ind;
            int *d_deg_where_d;
            int *d_deg_len;
            int *d_where;
            int *d_start_deg;

            float *d_deg_c;
            float *d_deg;

            cudaMalloc((void **)&d_deg_ind, sizeof(int) * 12 * g_count_deg_index );
            cudaMalloc((void **)&d_deg_where_d, sizeof(int) * g_count_imp_count );
            cudaMalloc((void **)&d_deg_len, sizeof(int) * g_count_imp_count );
            cudaMalloc((void **)&d_where, sizeof(int) * g_count_imp_symm );
            cudaMalloc((void **)&d_start_deg, sizeof(int) * g_count_imp_count );

            cudaMalloc((void **)&d_deg_c, sizeof(float) * g_count_deg_index );
            cudaMalloc((void **)&d_deg, sizeof(float) * g_count_imp_symm );

            cudaMemcpy(d_deg_len, deg_len, sizeof(int) * 
                                           g_count_imp_count, 
                                           cudaMemcpyHostToDevice );
            cudaMemcpy(d_deg_where_d, deg_where_d, sizeof(int) * 
                                                   g_count_imp_count, 
                                                   cudaMemcpyHostToDevice );
            cudaMemcpy(d_where, where, sizeof(int) * 
                                       g_count_imp_symm, 
                                       cudaMemcpyHostToDevice );  

            cudaMemcpy(d_deg_ind, deg_ind, sizeof(int) * 12 * 
                                           g_count_deg_index, 
                                           cudaMemcpyHostToDevice );
            cudaMemcpy(d_start_deg, start_deg, sizeof(int) * 
                                               g_count_imp_count, 
                                               cudaMemcpyHostToDevice );
            cudaMemcpy(d_deg_c, deg_c, sizeof(float) * 
                                       g_count_deg_index, 
                                       cudaMemcpyHostToDevice );
            cudaMemcpy(d_deg, deg, sizeof(float) * 
                                   g_count_imp_symm, 
                                   cudaMemcpyHostToDevice );
            

        //DEFINE CLOCK_T
        
            clock_t qprop_start, qprop_end;
            double qprop_time = 0.0;
            
            clock_t ini_sum_start, ini_sum_end;
            double ini_sum_time = 0.0;

            clock_t final_sum_start, final_sum_end;
            double final_sum_time = 0.0;

            clock_t cuker_sum_start, cuker_sum_end;
            double cuker_sum_time = 0.0;

            clock_t tdet_gpu, tdet_gpu_end;
            clock_t tsum_nt_start, tsum_nt_end;

            double tt_sum_nt = 0.0;
            double tt_det_gpu = 0.0;

            clock_t cufree_start, cufree_end;

         
        ini_end = clock();
        double ini_time = (double)((double)(ini_end-ini_start)/CLOCKS_PER_SEC);
        
        printf("INIT %f s \n", ini_time);
    //END INIT//


    //START NT_LOOP//
        clock_t ntloop_start, ntloop_end;
        ntloop_start = clock();
        printf("NT_LOOP\n");


        //START QPROP//
            qprop_start = clock();
            
            double complex *tqprop_new;

            tqprop_new = (double complex*)malloc(LEN_PROP_TEMP *
                                        sizeof(double complex));

            // tqprop = (double complex*)malloc(LEN_PROP * sizeof(double complex));

            memset(tqprop_new, 0.0, LEN_PROP_TEMP *
                                sizeof(double complex));

            double complex *tqprop;

            int len_tqprop = g_temp_nt * g_nx * g_ny * g_nz * 
                             g_nc * g_nd * g_nc * g_nd;
            tqprop = (double complex*)malloc(len_tqprop *
                                        sizeof(double complex));


            memset(tqprop, 0.0, len_tqprop *
                                sizeof(double complex));


            //START READ PROP//
                int fd = open("q_nx48_nt144", O_RDONLY);
                void *prop;

                prop = mmap(NULL, NRI * LEN_PROP_T * sizeof(double),
                            PROT_READ, MAP_PRIVATE, fd, 0);

                if (prop == MAP_FAILED) {
                    perror("mmap");
                    exit(EXIT_FAILURE);
                }
                
                double* prop_double = (double*) prop;  // cast void pointer to double pointer

                uint64_t idx, idx_re, idx_im;
                // idx = XDIM * T * ADIM * PDIM * ADIM * PDIM;
                uint64_t i, j, k, l, m, n;


                for (j = 0; j < g_temp_nt; j++) {
                    for (i = 0; i < XDIM; i++) {
                        for (m = 0; m < ADIM; m++) { //iic
                            for (n = 0; n < PDIM; n++) { //iid
                                for (k = 0; k < ADIM; k++) { //ifc
                                    for (l = 0; l < PDIM; l++) { //ifd
                                        idx =   j * XDIM * ADIM * PDIM * ADIM * PDIM + 
                                                i * ADIM * PDIM * ADIM * PDIM + 
                                                m * PDIM * ADIM * PDIM + 
                                                n * ADIM * PDIM +
                                                k * PDIM +
                                                l;


                                        idx_re =    n * ADIM * NRI * PDIM * ADIM * T * XDIM + 
                                                    m * NRI * PDIM * ADIM * T * XDIM + 
                                                    0 * PDIM * ADIM * T * XDIM + 
                                                    l * ADIM * T * XDIM + 
                                                    k * T * XDIM +
                                                    j * XDIM +
                                                    i;
                                        idx_im =    n * ADIM * NRI * PDIM * ADIM * T * XDIM + 
                                                    m * NRI * PDIM * ADIM * T * XDIM + 
                                                    1 * PDIM * ADIM * T * XDIM + 
                                                    l * ADIM * T * XDIM + 
                                                    k * T * XDIM +
                                                    j * XDIM +
                                                    i;

                                        tqprop[idx] = big_to_little(prop_double[idx_re]) + big_to_little(prop_double[idx_im]) * I ;


                                    }
                                }

                            }

                        }
                    }
                }

                munmap(prop, NRI * LEN_PROP_T * sizeof(double));
                close(fd);
            //END READ PROP//

            
            cuDoubleComplex *d_tqprop;  
            cudaMalloc(&d_tqprop, len_tqprop *
                                  sizeof(cuDoubleComplex));
            cudaMemcpy(d_tqprop,tqprop, len_tqprop *
                                        sizeof(cuDoubleComplex),
                                        cudaMemcpyHostToDevice);        

            qprop_end = clock();
            qprop_time += (double)((double)(qprop_end - qprop_start) / CLOCKS_PER_SEC);

            printf("-- QPROP %f s \n", qprop_time);
        //END QPROP//

        //START INIT_SUM//
            ini_sum_start = clock();
            double complex *sum_nt;
            sum_nt = (double complex*)malloc(g_temp_nt * 
                                             sizeof(cuDoubleComplex));

            memset(sum_nt, 0.0, g_temp_nt *
                                sizeof(double complex));

            double complex *sum_nxyz;
            sum_nxyz = (double complex*)malloc(g_temp_nt * g_nxyz * 
                                               sizeof(cuDoubleComplex));

            memset(sum_nxyz, 0.0, g_temp_nt * g_nxyz *
                                  sizeof(double complex));

            cuDoubleComplex *d_sum_nxyz; 
            cudaMalloc(&d_sum_nxyz, g_temp_nt * g_nxyz * 
                                    sizeof(cuDoubleComplex));
            cudaMemset(&d_sum_nxyz, 0.0, g_temp_nt * g_nxyz * 
                                         sizeof(cuDoubleComplex));


            ini_sum_end = clock();
            ini_sum_time += (double)((double)(ini_sum_end - ini_sum_start) / CLOCKS_PER_SEC);
            printf("-- INIT SUM %f s \n", ini_sum_time);        
        //END INIT_SUM//


        //START CUKER//
            tdet_gpu = clock();

            dim3 block(g_threads);
            dim3 grid(g_nxyz * g_temp_nt);

            //START CUKER_SUM//
                cuker_sum_start = clock();
                printf("Going inside cuKer");
                cuKer_sum <<< grid, block >>> (d_tqprop, d_sum_nxyz,
                                               d_deg_ind, d_deg_where_d, 
                                               d_deg_len, d_where, d_start_deg,
                                               d_deg_c, d_deg);
                cudaDeviceSynchronize();

                cuker_sum_end = clock();
                cuker_sum_time += (double)((double)(cuker_sum_end - cuker_sum_start) / CLOCKS_PER_SEC);
            //END CUKER_SUM//


            tdet_gpu_end = clock();
            tt_det_gpu = tt_det_gpu + (double)((double)(tdet_gpu_end - tdet_gpu) / CLOCKS_PER_SEC);
            
            printf("-- CUKER %f s \n", tt_det_gpu);
        //END CUKER//


        //START FINAL SUM//
            final_sum_start = clock();

            cudaMemcpy(sum_nxyz, d_sum_nxyz, g_temp_nt * g_nxyz *
                                             sizeof(double complex),
                                             cudaMemcpyDeviceToHost );        

            tsum_nt_start = clock();
            

            for (int sit = 0; sit < g_temp_nt; sit++) {
                for (int si = 0; si < g_nxyz; si++) {
                    sum_nt[sit] = sum_nt[sit] + sum_nxyz[si + sit * g_nxyz];
                }
            }
                
            tsum_nt_end = clock();
            tt_sum_nt += (double)((double)(tsum_nt_end - tsum_nt_start) / CLOCKS_PER_SEC);



            for (int sit = 0; sit < g_temp_nt; sit++) {
                printf("nt %d \t sum %.16E,%.16E \n", sit, 
                                                      creal(sum_nt[sit]), 
                                                      cimag(sum_nt[sit]));
            }    

            free(sum_nxyz);
            free(sum_nt);
            
            final_sum_end = clock();
            final_sum_time += (double)((double)(final_sum_end - final_sum_start) / CLOCKS_PER_SEC);  
            printf("-- FINAL SUM %f s \n", final_sum_time);
        //END FINAL SUM//


        ntloop_end = clock();
        double ntloop_time = (double)((double)(ntloop_end-ntloop_start)/CLOCKS_PER_SEC);
        printf("NT_LOOP %f s \n", ntloop_time);
    //END NT_LOOP//
    
    
    //START CUFREE//
        cufree_start = clock();

        cudaFree(d_deg_ind);
        cudaFree(d_deg_where_d);
        cudaFree(d_deg_len);
        cudaFree(d_where);
        cudaFree(d_start_deg);

        cudaFree(d_deg_c);
        cudaFree(d_deg);

        cudaFree(d_tqprop);
        cudaFree(d_sum_nxyz);

        cufree_end = clock();    

        double cufree_time = (double)((double)(cufree_end-cufree_start)/CLOCKS_PER_SEC);
        printf("CUFREE %f s \n", cufree_time);
    //END CUFREE//


    code_end = clock();
    //PRINT CPU TIME//
    
    double code_time = (double)((double)(code_end-code_start)/CLOCKS_PER_SEC);
    printf("TOTAL EXEC: %.5f s \n", code_time);
                

    //PRINT CU_DEV_RESET//
        clock_t cu_reset_start, cu_reset_end;
        cu_reset_start = clock();

        cudaDeviceReset();

        cu_reset_end = clock();

        double cu_reset_time = (double)((double)(cu_reset_end-cu_reset_start)/CLOCKS_PER_SEC);
        printf("CU_DEV_RESET: %4.9f s \n", cu_reset_time);


    return 0;

}    

np.cuh

#ifndef __NP_CUH_
#define __NP_CUH_

//for memset
#include <cstring>

#include <cuda.h>
#include <cuComplex.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include "read_files.h"



#endif

read_files.h

#ifndef __READ_FILES_H_
#define __READ_FILES_H_

#include <iostream>
#include <stdio.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <ctime>
#include <stdint.h>

#include <assert.h>
#include <cassert>

#include <sys/mman.h>

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>

#include <inttypes.h>


extern const uint64_t g_count_imp_symm;
extern const uint64_t g_count_deg_index;
extern const uint64_t g_count_imp_count;

extern const uint64_t g_nt;
extern const uint64_t g_nx;
extern const uint64_t g_ny;
extern const uint64_t g_nz;
extern const uint64_t g_nxyz;
extern const uint64_t g_nc;
extern const uint64_t g_nd;

extern const uint64_t msp;
extern const uint64_t g_msp;



int U = 3, D = 3;
int itl = 144, nc = 3, ns = 4, nri = 2, mdim = 4;
int nd = 4;

int nt = 144, nx = 48, ny = 48, nz = 48;

#define mx 48
#define my 48
#define mz 48
#define mt 144
#define msp mx*my*mz
#define ncs nc*ns



double  ******q;
FILE *quark1;

void f_read_imp_symm(int *where, float *deg);
void f_read_deg_index(int *deg_ind, float *deg_c);
void f_read_imp_count(int *deg_where_d, int *deg_len);

void f_read_prop();
void f_prop_compute(int a);
double f_read_8B_double ( FILE* fp);

void f_read_imp_symm(int *where, float *deg) {

    /* A: barrier 1- read imp_symm*/
    FILE *fptr1;
    fptr1 = fopen("imp_symm.txt", "r");

    for (int i = 0; i < g_count_imp_symm; i++) {
        fscanf(fptr1, "%d\n", &(where[i]));
        fscanf(fptr1, "%f\n", &(deg[i]));
    }
    fclose(fptr1);


}

void f_read_deg_index(int *deg_ind, float *deg_c) {
    FILE *fptr2;
    fptr2 = fopen("deg_index.txt", "r");

    for (int i = 0; i < g_count_deg_index; i++) {
        fscanf(fptr2, "%d %d %d %d %d %d %d %d %d %d %d %d %f\n",
               &(deg_ind[12 * i]),
               &(deg_ind[12 * i +  1]),
               &(deg_ind[12 * i +  2]),
               &(deg_ind[12 * i +  3]),
               &(deg_ind[12 * i +  4]),
               &(deg_ind[12 * i +  5]),
               &(deg_ind[12 * i +  6]),
               &(deg_ind[12 * i +  7]),
               &(deg_ind[12 * i +  8]),
               &(deg_ind[12 * i +  9]),
               &(deg_ind[12 * i +  10]),
               &(deg_ind[12 * i +  11]),
               &(deg_c[i]));
    }
    fclose(fptr2);


}

void f_read_imp_count(int *deg_where_d, int *deg_len, int *start_deg) {
    FILE *fptr2;
    fptr2 = fopen("imp_count.txt", "r");

    for (int i = 0; i < g_count_imp_count; i++) {
        fscanf(fptr2, "%d\n%d\n",
               &(deg_where_d[i]), &(deg_len[i]));
    }
    fclose(fptr2);


    int sum = 0;
    for (int i = 0; i < g_count_imp_count; i++)
    {
        start_deg[i] = sum;
        sum = sum + deg_len[i];
    }


}





#endif

A smaller -- CUKER 5.390000 s time.

Right now -- CUKER time stamp scales linearly with increasing g_temp_nt. Something that does not scale with g_temp_nt in np.cu will be great. In this code we change g_temp_nt which is currently set to 2 to g_nt = 144.

Answer 1

General Observations

This code isn't ready to be optimized, because it is not currently maintainable. Code optimization for performance is something that is done only when the code is production ready and some portion of the code is too slow. Using the current technology the best way to optimize C or C++ code is to use the optimizing feature of the compiler (the -O compiler switch). That should be the first step in optimizing the code since the compilers generally do a better job of optimization than programmers can do by hand.

To properly optimize code one needs to profile the code to find out where the bottle necks are. Profiling the code will tell you which functions are using the most time. Since this program really only has 3 functions that would be very difficult to determine where the bottlenecks are.

When this question was first posted, it had both C and C++ as language tags. As user Haris in a comment, the tags should be either one or the other but not both. I removed the C++ language tag because the code uses malloc() and free() as well as C language Input and Output (fopen() and flcose()). There is only one true C++ include file, iostream but no C++ I/O is used. The #include <iostream> should be removed from the code.

This question might also be considered off-topic because there are functions called that are not defined within the scope of the code such as :

• void f_read_prop();
• void f_prop_compute(int a);
• double f_read_8B_double(FILE* fp);

It also isn't clear that the functions defined in read_files.h are still in use since the large data file is now memory mapped.

Suggestions, break the code up into more functions, test the code and then post a follow up question with a link back to this one if you still need to optimize the code.

Programming and Maintenance Issues

Avoid Global Variables

It is very difficult to read, write, debug and maintain programs that use global variables. Global variables can be modified by any function within the program and therefore require each function to be examined before making changes in the code. In C and C++ global variables impact the namespace and they can cause linking errors if they are defined in multiple files. The answers in this stackoverflow question provide a fuller explanation.

There are a global variable define in read_files.h:

double****** q;
FILE* quark1;

If the header file read_files.h is included by more than one file this will cause linking errors because these variable will be defined in multiple files.

Header Files Should Only Contain Programming Interfaces

In the C programming language header files generally only contain typedefs, #define CONSTANTS, Enums and function prototypes, executable code should not be included in C header files, since the maintainance of the executable code would force recompile of all files that include that file. While executable code needs to be maintained, public interfaces generally remain static.

Test for Possible Memory Allocation Errors

In modern high-level languages such as C++, memory allocation errors throw an exception that the programmer can catch. This is not the case in the C programming language. While it is rare in modern computers because there is so much memory, memory allocation can fail, especially if the code is working in a limited memory application such as embedded control systems. In the C programming language when memory allocation fails, the functions malloc(), calloc() and realloc() return NULL. Referencing any memory address through a NULL pointer results in undefined behavior (UB).

Since there is a input file of 32G which is being memory mapped, memory allocation might fail. Some computers max out at 32G.

Possible unknown behavior in this case can be a memory page error (in Unix this would be call Segmentation Violation), corrupted data in the program and in very old computers it could even cause the computer to reboot (corruption of the stack pointer).

To prevent this undefined behavior a best practice is to always follow the memory allocation statement with a test that the pointer that was returned is not NULL.

    int *deg_ind = malloc(12 * g_count_deg_index * sizeof(*deg_ind));
    if (!deg_ind)
    {
        fprintf(stderr, "malloc of deg_ind failed, exiting program\n");
        return EXIT_FAILURE;
    }
    float *deg_c = malloc(1 * g_count_deg_index * sizeof(*deg_c));
    if (!deg_c)
    {
        fprintf(stderr, "malloc of deg_c failed, exiting program\n");
        return EXIT_FAILURE;
    }

Convention When Using Memory Allocation in C

When using malloc(), calloc() or realloc() in C a common convention is to sizeof(*PTR) rather sizeof(PTR_TYPE), this make the code easier to maintain and less error prone, since less editing is required if the type of the pointer changes. This is shown in the example above.

Complexity

The function main() is 374 lines of code, this is too complex (does too much). As programs grow in size the use of main() should be limited to calling functions that parse the command line, calling functions that set up for processing, calling functions that execute the desired function of the program, and calling functions to clean up after the main portion of the program.

A best practice in programming is to limit the size of one function or sub-routing to a single screen in an IDE or editor, the reason for this is that it is very difficult to keep track of what is going on in a function if it is larger than this. This code is very difficult to read due to the lack of functions, and that means it is very difficult to debug and maintain. One screen is generally 55 to 60 lines of code.

There is also a programming principle called the Single Responsibility Principle that applies here. The Single Responsibility Principle states:

that every module, class, or function should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by that module, class or function.

Duplication of Header Files

The code includes both assert.h and cassert, the cassert header file is assert.h. Use one or the other. In since this is primarily C code use assert.h.

Answer 2

Disclaimer

I don't know CUDA coding very well so I may miss some things.

This second answer is in response to questions on my first answer from the Original Poster.

The cuKer_sum() function can probably benefit from these observations as well.

Code that contains commented out lines is generally not considered ready for review.

Don't Ignore Return Values

Most CUDA functions, such as cudaMalloc() return an integer value that indicates whether the CUDA call was successful or not. There should be a test of this value to decide whether to continue the program or to terminate due to errors. Here is an example from a Stack Overflow Answer:

void *fixed_cudaMalloc(size_t len)
{
    void *p;
    if (cudaMalloc(&p, len) == success_code) return p;
    return 0;
}

If the code doesn't test the return value there may be undefined behavior in the code.

Possible Syntax Error in the Code

The function cudaMalloc() expects a pointer to a pointer so that the device memory (GPU memory) can be allocated, the following 2 lines of code in main() may be causing problems (undefined behavior) because the variable d_trial is an integer value rather than a pointer to an integer.

        int d_trial;
        cudaMalloc((void **)&d_trial, sizeof(int));

Programming and Maintenance Issues in cuKer_det()

The function cuker_det() is also too complex (does too much). This function could be much shorter by reducing the redundancy in the code using loops and functions.

DRY Code

There is a programming principle called the Don't Repeat Yourself Principle sometimes referred to as DRY code. If you find yourself repeating the same code multiple times it is better to encapsulate it in a function. If it is possible to loop through the code that can reduce repetition as well.

Based on the usage the following variable declarations in cuKer_det():

    int b1, b2, b3;
    int b1p, b2p, b3p;
    int bt1, bt2, bt3;
    int bt1p, bt2p, bt3p;

should probably be array declarations instead:

    int b[3];
    int bp[3];
    int bt[3];
    int btp[3];

#define BP_OFFSET 3
#define BT_OFFSET 6
#define BTP_OFFSET 9
 
    for (size_t i = 0; i < 3; i++)
    {
        b[i] = d_deg_ind[(12 * tx) + i];
        bp[i] = d_deg_ind[(12 * tx) + i + BP_OFFSET];
        bt[i] = d_deg_ind[(12 * tx) + i + BT_OFFSET];
        btp[i] = d_deg_ind[(12 * tx) + i + BT_OFFSET];
    }

The assignments to the d_A array are also repetitive and could be in a function or a function called in a loop.

Write Self Documenting Code

It isn't clear in the code what the variables b1, b2, b3, b1p ... are, I or anyone else who needs to maintain this code would be completely lost as to what this function is doing. This is also true of the variables d_deg_ind, d_tqprop, d_A etc. The function names cuKer_det() and cuKer_sum() are only slightly more informative (What is cuKer?). If you went back to this code in 6 months to a year, you might not understand it anymore, someone that didn't originally write the code has no chance of understanding.

Horizontal Spacing

To make code more readable it is common to put spaces between operators and operands in expressions.

    d_A[0*3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
                            b1*4*3*4 + bt1*3*4 + b1p*4 + bt1p];

is more difficult to read than

    d_A[0 * 3 + 0] = d_tqprop[gid * g_nc * g_nc * g_nd * g_nd +
        b1 * 4 * 3 * 4 + bt1 * 3 * 4 + b1p * 4 + bt1p];

It isn't clear why there is multiplication of constants in the above expressions, why is there b1 * 4 * 3 *4, bt1 * 3 * 4. Why not just b1 * 48 and bt1 * 12.

Magic Numbers

At the top of the program a number of constants are defined, and this is very good, but the practice is not maintained through out the code. One of the most important things a coder should do is be consistent. The function cuKer_det() has a lot of raw numbers in it (12, 1 through 12, 4 * 3 * 4, 4 * 3, 4), it might be better to create symbolic constants for them to make the code more readable and easier to maintain. These numbers may be used in many places and being able to change them by editing only one line makes maintenance easier.

Numeric constants in code are sometimes referred to as Magic Numbers, because there is no obvious meaning for them. There is a discussion of this on stackoverflow.

A common practice in the C programming language is that constants should be ALL CAPITALS, so that they can easily be found.

What cuker_det() Might Look Like Following DRY Code Recommendations:

Since it isn't completely clear what cuKer_det() is doing, the following rewrite of the function may be incorrect. It is definitely possible to reduce the number of lines in the code and the code complexity in the function.

__device__ cuDoubleComplex cuKer_det(int gid, int tx,
    cuDoubleComplex* d_tqprop,
    int* d_deg_ind,
    float* d_deg_c)
{
    int b[3];
    int bp[3];
    int bt[3];
    int btp[3];

#define BP_OFFSET 3
#define BT_OFFSET 6
#define BTP_OFFSET 9

    for (size_t i = 0; i < 3; i++)
    {
        b[i] = d_deg_ind[(12 * tx) + i];
        bp[i] = d_deg_ind[(12 * tx) + i + BP_OFFSET];
        bt[i] = d_deg_ind[(12 * tx) + i + BT_OFFSET];
        btp[i] = d_deg_ind[(12 * tx) + i + BT_OFFSET];
    }
    cuDoubleComplex d_A[9];

    size_t g_ncsqr = 9;
    size_t g_ndsqr = 16;

    for (size_t i = 0; i < 3; i++)
    {
        for (size_t j = 0; j < 3; j++)
        {
            d_A[(i * 3) + j] = d_tqprop[gid * g_ncsqr * g_ndsqr +
                b[i] * 4 * 3 * 4 + bt[i] * 3 * 4 + bp[j] * 4 + btp[j]];
        }
    }

    cuDoubleComplex x1 = cuCmul(d_A[0 * 3 + 0],
        cuCsub(cuCmul(d_A[1 * 3 + 1], d_A[2 * 3 + 2]),
            cuCmul(d_A[1 * 3 + 2], d_A[2 * 3 + 1])));
    cuDoubleComplex x2 = cuCmul(d_A[0 * 3 + 1],
        cuCsub(cuCmul(d_A[1 * 3 + 0], d_A[2 * 3 + 2]),
            cuCmul(d_A[1 * 3 + 2], d_A[2 * 3 + 0])));
    cuDoubleComplex x3 = cuCmul(d_A[0 * 3 + 2],
        cuCsub(cuCmul(d_A[1 * 3 + 0], d_A[2 * 3 + 1]),
            cuCmul(d_A[1 * 3 + 1], d_A[2 * 3 + 0])));

    cuDoubleComplex r1 = make_cuDoubleComplex(d_deg_c[tx], 0.0);
    cuDoubleComplex x123 = cuCadd(cuCsub(x1, x2), x3);
    cuDoubleComplex r1x123 = cuCmul(r1, x123);

    return r1x123;
}

Answer 3

improving the efficiency of the code

I'll only address that goal.

Move redundant calculations

Example:

// for (l = 0; l < PDIM; l++) { //ifd
//   idx = j * XDIM * ADIM * PDIM * ADIM * PDIM + 
//         i * ADIM * PDIM * ADIM * PDIM + 
//         m * PDIM * ADIM * PDIM + 
//         n * ADIM * PDIM +
//         k * PDIM +
//         l;
uint64_t sum = j * XDIM * ADIM * PDIM * ADIM * PDIM + 
               i * ADIM * PDIM * ADIM * PDIM + 
               m * PDIM * ADIM * PDIM + 
               n * ADIM * PDIM +
               k * PDIM;
for (l = 0; l < PDIM; l++) { //ifd
  idx = sum + l;

This change could be extended to moving parts of the computation to prior parts of each nested for loop.

A Good Compile may do this anyways, yet worth profiling this alternative.

Use size_t instead of uint64_t for array indexing

size_t may be narrower than uint64_t. For array sizing and indexing, math wider than size_t is usually not needed.

Minor: Use const

When referenced data is not modified, use const to let the compiler know that the data is not changed. Good compilers will see this anyways, yet lesser compilers will optimize better.

Review all user functions for this.

Example

// __device__ cuDoubleComplex cuKer_det(
//     int gid, int tx, 
//     cuDoubleComplex *d_tqprop, int *d_deg_ind, float *d_deg_c) {
__device__ cuDoubleComplex cuKer_det(
    int gid, int tx, 
    const cuDoubleComplex *d_tqprop, const int *d_deg_ind, const float *d_deg_c) {

Advanced: Use restrict

This effectively asserts that the pointer does not reference data that is overlapped by others somehow. This may help a lot or not much at all. Use this for all pointers that are meant to point to non-overlapping data.

    // cuDoubleComplex *d_sum_nxyz
    cuDoubleComplex * restrict d_sum_nxyz

It is the difference between memcpy() and memmove() and allows for select optimizations the compiler cannot assume.

Use float

Unless double precision and range needed, consider float, float complex. Be sure then to use float constants and float function calls then too.

I'd expect this to make code 4x as fast.