0
\$\begingroup\$

I've constructed a system that contains multiple 2D lattice and a 1D cable all in a single vector of type double. Each 2D lattice is coupled to a site along the cable. For instance, there are 800 nodes on the cable which stores the state of each 800 50-by-50 lattice that it is coupled to. The size of the vector is then 800*50*50. Using CUDA C, I implement a generalized "stencil" kernel function that computes the average of n - nearest neighbors of a given index. In this example, the stencil is generalized to be used for both the 1D cable and the 2D lattice. The code runs fine, except it takes longer than the code without CUDA! Obviously, I am not implementing it correctly. The following snippet demonstrates CUDA protocol. The OFFSETS values take the iterator to the start and end points of one of the 2D lattice that is being operated on.

first = v.begin() + OFFSET0;    
last = v.begin() + OFFSET1; 
std::copy(first,last,host_vector);
cudaMemcpy(d_in,host_vector,size,cudaMemcpyHostToDevice);
_1Dstencil<<<M,THREADS_PER_BLOCK>>>(d_in,d_out,nx,ny,nz,X_,Y_,Z_);
cudaMemcpy(host_vector,d_out,size,cudaMemcpyDeviceToHost);

The following is my implementation of the generalized kernel stencil

__global__ void _1Dstencil(double* d_in,double* d_out,int nx,int ny,int nz,int X_,int Y_,int Z_){
    int i = threadIdx.x+ blockIdx.x* blockDim.x ;       
    int i_ = (i/(Z_*Y_))%X_ , 
        j_ = (i/Z_)%Y_, 
        k_ = i%Z_ ;
    double sum = 0;
    //stencil around i  
    for(int nnx =   2*nx+1; nnx> 0; nnx--){//X      
        int start = i_-(nnx-nx-1);
        int ii = cp(start,X_,false);        
        for(int nny = 2*ny+1; nny> 0; nny--){//Y
            start = j_-(nny-ny-1);
            int jj = cp(start,Y_,false);
            for(int nnz = 2*nz+1; nnz> 0; nnz--){//Z
                    start = k_-(nnz-nz-1);
                    int kk = cp(start,Z_,false);
                    sum += *(d_in + _ND_1D(ii, jj, kk, X_, Y_, Z_));
            }
        }
    }
    d_out[i] = sum/((2*nx+1)*(2*ny+1)*(2*nz+1));
}

The functions cp(...) and _ND_1D(...) are declared as device functions which return the correct position based on flux boundaries and convert N-Dimensional coordinates to 1D coordinates, respectively. Could someone help me clarify the correct way of distributing the task on the deice? Thanks in advance!

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.