/home/julian/svn/gpucpumd/mdgpu2/gpu_tools.cu

Go to the documentation of this file.

#include "gpu_tools.cuh"


//####################################### Verlet List functions ####################################//


__global__ void gpuGenerateVerletList(FLOAT_ARRAY_TYPE *dPos, int *dVlists, size_t pitchDVlists, int *dVcount, float *dmax, float rListSq, int N, float boxSize, bool *doUpdate)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE d;
        dVcount[idx]=0;
        int* dVlist;
        for (int i=0; i<N; ++i)
        {
            if ( i != idx )
            {
                gpuDistPBC(&d, &dPos[idx], &dPos[i], &boxSize);

                if ( rListSq > (d.x * d.x + d.y * d.y + d.z *d.z) )
                {
                    dVlist = (int*)((char*)dVlists + dVcount[idx] * pitchDVlists);
                    dVlist[idx] = i;
                    dVcount[idx] = dVcount[idx] + 1;
                }
            }
        }
        dmax[idx]=0.0f;
    }
    //reset update flag
    (*doUpdate)=false;
}


__global__ void gpuGenerateVerletListVar1(FLOAT_ARRAY_TYPE *dPos, int *dVlists, size_t pitchDVlists, int *dVcount, float *dmax, float rListSq, int N, float boxSize, bool *doUpdate)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx]; //save reference particle position in register
        FLOAT_ARRAY_TYPE d;
        dVcount[idx]=0;
        int* dVlist;
        for (int i=0; i<N; ++i)
        {
            if ( i != idx )
            {
                gpuDistPBC(&d, &tPos, &dPos[i], &boxSize);
                if ( rListSq > (d.x * d.x + d.y * d.y + d.z *d.z) )
                {
                    dVlist = (int*)((char*)dVlists + dVcount[idx] * pitchDVlists);
                    dVlist[idx] = i;
                    dVcount[idx] = dVcount[idx] + 1;
                }
            }
        }
        dmax[idx]=0.0f;
    }
    //reset update flag
    (*doUpdate)=false;
}


__global__ void gpuGenerateVerletListSmem(FLOAT_ARRAY_TYPE *dPos, int *dVlists, size_t pitchDVlists, int *dVcount, float *dmax, float rListSq, int N, float boxSize, bool *doUpdate)
{
    //CONSTRAINT: numBlocks * numThreadsPerBlock == N, otherwise __synchthread() called can not work!
    extern __shared__ FLOAT_ARRAY_TYPE sPos[]; //allocate shared memory
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + tid;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE tPos = dPos[idx];//save reference particle position in register
        FLOAT_ARRAY_TYPE d;
        dVcount[idx]=0;
        int* dVlist;

        for(int j=0; j<gridDim.x; ++j)
        {
            int bOffset = j * blockDim.x;
            sPos[tid] = dPos[bOffset+tid]; //fetch block of particle positions and save in shared memory
            __syncthreads();
            for (int i=0; i<blockDim.x; ++i)
            {
                int id = bOffset + i;
                if ( id != idx && id < N)
                {
                    gpuDistPBC(&d, &tPos, &sPos[i], &boxSize); //use broadcast from shared memory -> full coalescence
                    if ( rListSq > (d.x * d.x + d.y * d.y + d.z *d.z) )
                    {
                        dVlist = (int*)((char*)dVlists + dVcount[idx] * pitchDVlists);
                        dVlist[idx] = id;
                        dVcount[idx] = dVcount[idx] + 1;
                    }
                }
            }
            __syncthreads();
        }
        dmax[idx]=0.0f;
    }
    //reset update flag
    (*doUpdate)=false;
}



__global__ void gpuGenerateVerletListSmemVar1(FLOAT_ARRAY_TYPE *dPos, int *dVlists, size_t pitchDVlists, int *dVcount, float *dmax, float rListSq, int N, float boxSize, bool *doUpdate)
{
    extern __shared__ FLOAT_ARRAY_TYPE sPos[]; //allocate shared memory
    int tid = threadIdx.x;
    int offset = blockDim.x * blockIdx.x;
    int idx = offset + tid;
    if (idx < N)
    {
        sPos[tid] = dPos[idx];//save reference particle position in smem
        __syncthreads();
        FLOAT_ARRAY_TYPE d;
        dVcount[idx]=0;
        int* dVlist;

        //block before thread block in device memory
        for(int j=0; j< offset; ++j)
        {
            gpuDistPBC(&d, &sPos[tid], &dPos[j], &boxSize);
            if ( rListSq > (d.x * d.x + d.y * d.y + d.z *d.z) )
            {
                dVlist = (int*)((char*)dVlists + dVcount[idx] * pitchDVlists);
                dVlist[idx] = j;
                dVcount[idx] = dVcount[idx] + 1;
            }
        }
        //block in shared memory
        for(int j=0; j<blockDim.x; ++j)
        {
            if(j != tid)//exclude self interaction
            {
                gpuDistPBC(&d, &sPos[tid], &sPos[j], &boxSize);
                if ( rListSq > (d.x * d.x + d.y * d.y + d.z *d.z) )
                {
                    dVlist = (int*)((char*)dVlists + dVcount[idx] * pitchDVlists);
                    dVlist[idx] = offset + j;
                    dVcount[idx] = dVcount[idx] + 1;
                }
            }
        }
        //block after threadblock in device memory
        for(int j=offset+blockDim.x; j<N; ++j)
        {
            gpuDistPBC(&d, &sPos[tid], &dPos[j], &boxSize);
            if ( rListSq > (d.x * d.x + d.y * d.y + d.z *d.z) )
            {
                dVlist = (int*)((char*)dVlists + dVcount[idx] * pitchDVlists);
                dVlist[idx] = j;
                dVcount[idx] = dVcount[idx] + 1;
            }
        }
        dmax[idx]=0.0f;
    }
    //reset update flag
    (*doUpdate)=false;
}


//####################################### NVT ensemble methods ####################################//


void rescaleT(int N, int dim, int numInterData, int numThreads, float T, int threads, int blocks, FLOAT_ARRAY_TYPE *dVel, float *sumSq)
{
    //get sum of squares of velocities
    SqVecSumReduction(dVel, sumSq, numInterData, numThreads, N);

    //rescale velocities to desired temperatures
    gpuRescaleTKernel<<<blocks, threads>>>(N, dim, T, dVel, sumSq);
    checkCudaError("Kernel error: gpuRescaleTKernel");
}

__global__ void gpuRescaleTKernel(int N, int dim, float T, FLOAT_ARRAY_TYPE *dVel, float *vSumSq)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    float factor = __fsqrt_rn(dim*N*T/vSumSq[0]);
    if (idx < N)
    {
        dVel[idx].x = dVel[idx].x * factor;
        dVel[idx].y = dVel[idx].y * factor;
        dVel[idx].z = dVel[idx].z * factor;
    }
}


//####################################### general MD helper functions ####################################//

void removeDrift(int N, int numInterData, int numThreads, int threads, int blocks, FLOAT_ARRAY_TYPE *dVel, FLOAT_ARRAY_TYPE *compSum)
{
    //get component wise sum of positon
    componentSumReduction(dVel, compSum, numInterData, numThreads, N);

    //remove drift
    gpuRemoveDriftKernel<<<blocks, threads>>>(N, dVel, compSum);
    checkCudaError("Kernel error: gpuRemoveDriftKernel");
}

__global__ void gpuRemoveDriftKernel(int N, FLOAT_ARRAY_TYPE *dVel, FLOAT_ARRAY_TYPE *compSum)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    FLOAT_ARRAY_TYPE factor;
    factor.x = -compSum[0].x/N;
    factor.y = -compSum[0].y/N;
    factor.z = -compSum[0].z/N;
    if (id < N)
    {
        dVel[id].x += factor.x;
        dVel[id].y += factor.y;
        dVel[id].z += factor.z;
    }
}

__global__ void gpuGetVelOnTime(FLOAT_ARRAY_TYPE *dVel, FLOAT_ARRAY_TYPE *dVt, FLOAT_ARRAY_TYPE *dForce, int N, float dt)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N)
    {
        dVt[idx].x = dVel[idx].x + 0.5f * dt * dForce[idx].x;
        dVt[idx].y = dVel[idx].y + 0.5f * dt * dForce[idx].y;
        dVt[idx].z = dVel[idx].z + 0.5f * dt * dForce[idx].z;
    }
}



//####################################### Observable functions ####################################//



float getKineticEnergy(FLOAT_ARRAY_TYPE *dVel, int &N, float *sumSq, int &maxThreads, int &maxBlocks)
{
    //get sum of squares
    SqVecSumReduction(dVel, sumSq, 64,64, N);
    float e;
    cudaMemcpy(&e, sumSq, sizeof(float), cudaMemcpyDeviceToHost);
    e *= 0.5f;
    return e;
}





float getHarmonicPotentialEnergy(float *dEpot, FLOAT_ARRAY_TYPE *dPos, float *sum, int *cnn, size_t &pitch, int &N, int &numNN, float &k, float &boxSize, dim3 &dimGrid, dim3 &dimBlock, int &maxThreads, int &maxBlocks)
{
    float e;
    gpuGetHarmonicPotentialEnergyKernel<<< dimGrid, dimBlock>>>(dEpot, dPos, cnn, pitch, N, numNN, k, boxSize);
    checkCudaError("Kernel error: gpu_get_harmonic_potential");
    sumReduction(dEpot, sum, 64, 64, N);
    cudaMemcpy(&e, sum, sizeof(float), cudaMemcpyDeviceToHost);
    return e;
}


__global__ void gpuGetHarmonicPotentialEnergyKernel(float *dEpot, FLOAT_ARRAY_TYPE *dPos, int *cnn, size_t pitch, int N, int numNN, float k, float boxSize)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    FLOAT_ARRAY_TYPE d;
    if (idx < N)
    {
        dEpot[idx] = 0.0f;
        int pid;
        int *row = (int*)((char*)cnn + idx * pitch);
        for (int c = 0; c < numNN; ++c)
        {
            pid = row[c];
            if (pid != -1)
            {
                gpuDistPBC(&d, &dPos[idx], &dPos[pid], &boxSize);
                dEpot[idx] += d.x * d.x + d.y * d.y + d.z * d.z;
            }
        }
        dEpot[idx] *= 0.25f*k;
    }
}


float getLJCutPotentialEnergy(FLOAT_ARRAY_TYPE *dPos, float *dEpot, float *sum, int &N, float &boxSize, float &rCutSq, dim3 &dimGrid, dim3 &dimBlock, int &maxThreads, int &maxBlocks)
{
    float e;
    gpuLJCutPotentialEnergyKernel<<< dimGrid, dimBlock>>>(dPos, dEpot, N, boxSize, rCutSq);
    checkCudaError("Kernel error: gpuLJCutPotentialEnergyKernel");
    sumReduction(dEpot, sum, 64, 64, N);
    cudaMemcpy(&e, sum, sizeof(float), cudaMemcpyDeviceToHost);
    return e;
}


__global__ void gpuLJCutPotentialEnergyKernel(FLOAT_ARRAY_TYPE *dPos, float *dEpot, int N, float boxSize, float rCutSq)
{
    //Reduced units are used -> epsilon=sigma=1
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < N)
    {
        FLOAT_ARRAY_TYPE d;
        float invR6;
        dEpot[idx]=0.0f;//clear
        for (int i=0; i<N; ++i)
        {
            if ( i != idx) //excude self interaction
            {
                gpuDistPBC( &d, &dPos[idx],  &dPos[i], &boxSize);
                invR6 = d.x * d.x + d.y * d.y + d.z * d.z;
                if ( rCutSq > invR6 )
                {
                    invR6 = 1.0f / invR6;
                    invR6 = powf( invR6, 3.0f );
                    dEpot[idx] += 2.0f * invR6 * (invR6 - 1.0f);
                }
            }
        }
    }
}

//####################################### CUDA helper functions ####################################//


inline static unsigned long inKB(unsigned long bytes)
{
    return bytes/1024;
}

inline static unsigned long inMB(unsigned long bytes)
{
    return bytes/(1048576);
}


bool memInfo()
{
    bool success=true;
    unsigned int free, total;
    CUresult res;
    res = cuMemGetInfo(&free, &total);
    if (res != CUDA_SUCCESS)
    {
        std::cout<<"uMemGetInfo failed: "<<res<<"\n";
        success=false;
    }
    else
    {
        std::cout<<"________________ Memory usage ________________\n";
        std::cout<<"Free : "<<free<<" bytes ("<<inKB(free)<<" KB) ("<<inMB(free)<<" MB)\n";
        std::cout<<"Used : "<<total-free<<" bytes ("<<inKB(total-free)<<" KB) ("<<inMB(total-free)<<" MB)\n";
        std::cout<<"Total: "<<total<<" bytes ("<<inKB(total)<<" KB) ("<<inMB(total)<<" MB)\n";
        std::cout<<"-----> "<<100.0*free/(double)total<<"% free, "<<100.0*(total - free)/(double)total<<"% used\n";
        std::cout<<"______________________________________________\n";
    }
    return success;
}


void checkCudaError(const char *identifier)
{
    cudaThreadSynchronize();
    cudaError_t err = cudaGetLastError();
    if ( cudaSuccess != err)
    {
        std::cout<<"Cuda error: "<<identifier<<" : "<< cudaGetErrorString( err) <<"\n";
        exit(EXIT_FAILURE);
    }
}



__global__ void copyDevToDev(FLOAT_ARRAY_TYPE *dA, FLOAT_ARRAY_TYPE *dB, int N)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < N)
    {
        dB[idx].x = dA[idx].x;
        dB[idx].y = dA[idx].y;
        dB[idx].z = dA[idx].z;
    }
}


__global__ void gpuSimpleSqSum(FLOAT_ARRAY_TYPE *dVal, float*sqSum, int N)
{
    int idx=blockIdx.x*blockDim.x+threadIdx.x;
    if (idx==0)
    {
        (*sqSum)=0.0f;
        for (int i=0; i<N; ++i)
        {
            (*sqSum) += dVal[i].x * dVal[i].x + dVal[i].y * dVal[i].y + dVal[i].z * dVal[i].z;
        }
    }
}

__global__ void gpuSimpleComponentSum(FLOAT_ARRAY_TYPE *dVal, FLOAT_ARRAY_TYPE *sum, int N)
{
    int idx=blockIdx.x*blockDim.x+threadIdx.x;
    if (idx==0)
    {
        (*sum).x=0.0f;
        (*sum).y=0.0f;
        (*sum).z=0.0f;
        for (int i=0; i<N; ++i)
        {
            (*sum).x += dVal[i].x;
            (*sum).y += dVal[i].y;
            (*sum).z += dVal[i].z;
        }
    }
}



//####################################### Screen output functions ####################################//


void showPositions(FLOAT_ARRAY_TYPE *hPos, int N)
{
    for (int i=0; i<N; ++i)
    {
        std::cout<<"hPos["<<i<<"] = "<<hPos[i].x<<" , "<<hPos[i].y<<" , "<<hPos[i].z<<"\n";
    }
}

void showVelocities(FLOAT_ARRAY_TYPE *hVel, int &N)
{
    for (int i=0; i<N; ++i)
    {
        std::cout<<"hVel["<<i<<"] = "<<hVel[i].x<<" , "<<hVel[i].y<<" , "<<hVel[i].z<<"\n";
    }
}

void showForces(FLOAT_ARRAY_TYPE *hForce, int &N)
{
    for (int i=0; i<N; ++i)
    {
        std::cout<<"hForce["<<i<<"] = "<<hForce[i].x<<" , "<<hForce[i].y<<" , "<<hForce[i].z<<"\n";
    }
}

void showDeviceArray(FLOAT_ARRAY_TYPE *dArray, int N)
{
    FLOAT_ARRAY_TYPE *hArray = new FLOAT_ARRAY_TYPE[N];
    cudaMemcpy(hArray, dArray, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyDeviceToHost);
    std::cout<<"---- Device array ----\n\n";
    for (int i=0; i<N; ++i)
    {
        std::cout<<i<<": "<<hArray[i].x<<" "<<hArray[i].y<<" "<<hArray[i].z<<"\n";
    }
    std::cout<<"----------------------\n\n";
    delete[] hArray;
}

void showNeareastNeighborList(int *dNNListArray, int &N, int &numNN)
{
    for (int i=0; i<N; ++i)
    {
        int id=numNN*i;
        std::cout<<"---------------------\nparticle: "<<i<<"\nneighbors:\n";
        for (int j=0; j<numNN; ++j)
        {
            std::cout<<dNNListArray[id+j]<<"\n";
        }
    }
    std::cout<<"------------------------\n";
    std::cout<<"------------------------\n";
}

void showVerletList(int *dVlists, size_t pitchVlists, int *dVcount, int N, int NList)
{
    int *hVerletList=new int[N*NList];
    int *hCount = new int[N];
    cudaMemcpy(hCount, dVcount, N * sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy2D(hVerletList, N*sizeof(int), dVlists, pitchVlists, N * sizeof(int), NList, cudaMemcpyDeviceToHost);
    for (int i=0; i<10; ++i)
    {
        std::cout<<"i: "<<i<<" ; vcount: "<<hCount[i]<<" -->";
        for (int j=0; j<hCount[i]; ++j)
        {
            std::cout<<" "<<hVerletList[i + j*N];
        }
        std::cout<<"\n_________________\n";
    }
    int maxentry=0;
    int indexmax;
    for (int i=0; i<N; ++i)
    {
        if (maxentry < hCount[i])
        {
            maxentry = hCount[i];
            indexmax=i;
        }
    }
    std::cout<<"max entry: "<<maxentry<<" @ index: "<<indexmax<<"\n";
    delete[] hVerletList;
    delete[] hCount;
}

//####################################### Simulation parameter functions ####################################//

bool setCudaDevice(int gpuDevNo)
{
    //check if desired device ID is there
    int deviceCount;
    cudaGetDeviceCount(&deviceCount);
    if(deviceCount <= gpuDevNo)
    {
        std::cout<<"Device ID '"<<gpuDevNo<<"' not present!\n";
        return false;
    }

    cudaSetDevice( gpuDevNo );
    cudaDeviceProp prop;         // get device properties
    cudaGetDevice(&gpuDevNo);
    cudaGetDeviceProperties(&prop, gpuDevNo);
    std::cout<<"_______________________________________________________\n";
    std::cout<<"Number of GPU devices: "<<deviceCount<<"\n";
    std::cout<<"Running on GPU device #: "<<gpuDevNo<<"\n";
    std::cout<<"GPU processor #: "<<8*prop.multiProcessorCount<<"\n";
    std::cout<<"Computing Capability: ["<<prop.major<<"."<< prop.minor<<"]\n";
    std::cout<<"_______________________________________________________\n";
    return true;
}



bool generateStartScript(std::string exeName)
{
    bool ok;
    std::string scriptname = "example.sh";
    std::string inputfile = "N16_T0.5_g2.0_k0.1_init.dat";
    std::string basename = "N16_T0.5_g2.0_k0.1_";
    int dim = 3;
    int N = 4096;
    int numNN = 6;
    float T = 0.5;
    float k = 0.1;
    float boxSize = 32;
    float dt = 0.004;
    int NOffset = 0;
    int NSim = 10000;
    int NOut = 1000;
    int NAdjust = 100;
    int gpuDevNo=0;
    ok =        generateScript(scriptname, exeName, inputfile, basename, dim, N, numNN, T, k, boxSize, dt, NOffset, NSim, NOut, NAdjust, gpuDevNo);
    return ok;
}


bool generateRestartScript(std::string exeName, std::string basename, int dim, int N, int numNN, float T, float k, float boxSize, float dt, int NOffset, int NSim, int NOut, int NAdjust, int gpuDevNo)
{
    std::ostringstream ssFile;
    std::string scriptname;
    ssFile<<basename<<"_continue.sh";
    scriptname = ssFile.str();
    ssFile.str("");
    ssFile.clear();
    ssFile.fill('0');
    ssFile<<basename<<"";
    ssFile.width(10);
    ssFile<<NOffset;
    ssFile.width(1);
    ssFile<<".dat";
    std::string inputfile = ssFile.str();
    bool ok;
    ok =        generateScript(scriptname, exeName, inputfile, basename, dim, N, numNN, T, k, boxSize, dt, NOffset, NSim, NOut, NAdjust, gpuDevNo);
    return ok;
}

bool generateScript(std::string scriptname, std::string exeName, std::string inputfile, std::string basename, int dim, int N, int numNN, float T, float k, float boxSize, float dt, int NOffset, int NSim, int NOut, int NAdjust, int gpuDevNo)
{
    std::ofstream out(scriptname.c_str());
    if (!out.good())
    {
        std::cout<<"'"<<scriptname<<"' couldn't be created!\n";
        return 0;
    }
    out<<"#!/bin/bash\n\n";
    out<<"# input x-v-data file name\n";
    out<<"input="<<inputfile<<" \n\n";
    out<<"# output base filename\n";
    out<<"basename="<<basename<<"\n\n";
    out<<"# Dimension of the simulation \n";
    out<<"dim="<<dim<<"\n\n";
    out<<"# particle number \n";
    out<<"N="<<N<<"\n\n";
    out<<"# maximum number of coupled nearest neighbors\n";
    out<<"num_nn="<<numNN<<"\n\n";
    out<<"# temperature\n";
    out<<"temperature="<<T<<"\n\n";
    out<<"# spring constant:\n";
    out<<"k="<<k<<"\n\n";
    out<<"# edge length of the simulation box\n";
    out<<"box="<<boxSize<<"\n\n";
    out<<"# timestep\n";
    out<<"timestep="<<dt<<"\n\n";
    out<<"# offset for simulation steps\n";
    out<<"Noffset="<<NOffset<<"\n\n";
    out<<"# simulate till NSim steps are done\n";
    out<<"Nsim="<<NSim<<"\n\n";
    out<<"# output intervall (CONSTRAINT: must be a multiple of 'Nadjust')\n";
    out<<"Nout="<<NOut<<"\n\n";
    out<<"# NVT velocity scaling and drift removal intervall\n";
    out<<"Nadjust="<<NAdjust<<"\n\n";
    out<<"# number of GPU device to run the job on\n";
    out<<"gpuDevNo="<<gpuDevNo<<"\n\n##################\n\n";
    out<<exeName<<" $input $basename $dim $N $num_nn $temperature $k $box $timestep $Noffset $Nsim $Nout $Nadjust $gpuDevNo\n";
    return 0;
}


//####################################### Parallel Primitives ####################################//



__global__ void gpuReduce1(float *dInData, float *dInterData, int numInterData, int numThreads, int N)
{
    extern __shared__ float smem[];
    int idx=blockIdx.x*blockDim.x+threadIdx.x;
    float s=0.0f;
    for (int j=idx; j<N; j+=numInterData*numThreads)
    {
        s += dInData[j];
    }
    smem[threadIdx.x]=s;
    __syncthreads();
    if (threadIdx.x==0)
    {
        float s=0.0f;
        for (int i=0; i<numThreads; i++)
        {
            s+=smem[i];
        }
        dInterData[blockIdx.x]=s;
    }
}


__global__ void gpuReduce2(float *dInterData, int numInterData, int numThreads)
{
    extern __shared__ float smem[];
    int idx=threadIdx.x;
    float s=0.0f;
    for (int j=idx; j<numInterData; j+=numThreads)
    {
        s+=dInterData[j];
    }
    smem[threadIdx.x]=s;
    __syncthreads();
    if (idx==0)
    {
        float s=0.0f;
        for (int i=0; i<numThreads; i++)
        {
            s+=smem[i];
        }
        dInterData[0]=s;
    }
}



__global__ void gpuSqSumReduce1(FLOAT_ARRAY_TYPE *dInData, float *dInterData, int numInterData, int numThreads, int N)
{
    extern __shared__ float smem[];
    int idx=blockIdx.x*blockDim.x+threadIdx.x;
    float s=0.0f;
    for (int j=idx; j<N; j+=numInterData*numThreads)
    {
        s += dInData[j].x * dInData[j].x + dInData[j].y * dInData[j].y + dInData[j].z * dInData[j].z;
    }
    smem[threadIdx.x]=s;
    __syncthreads();
    if (threadIdx.x==0)
    {
        float s=0.0f;
        for (int i=0; i<numThreads; i++)
        {
            s+=smem[i];
        }
        dInterData[blockIdx.x]=s;
    }
}


__global__ void gpuComponentSumReduce1(FLOAT_ARRAY_TYPE *dInData, FLOAT_ARRAY_TYPE *dInterData, int numInterData, int numThreads, int N)
{
    extern __shared__ FLOAT_ARRAY_TYPE smema[];
    int tid = threadIdx.x;
    int bid = blockIdx.x;
    int idx=bid*blockDim.x+tid;
    FLOAT_ARRAY_TYPE s;
    s.x = 0.0f;
    s.y = 0.0f;
    s.z = 0.0f;
    for (int j=idx; j<N; j+=numInterData*numThreads)
    {
        s.x += dInData[j].x;
        s.y += dInData[j].y;
        s.z += dInData[j].z;
    }
    smema[tid].x = s.x;
    smema[tid].y = s.y;
    smema[tid].z = s.z;
    __syncthreads();
    if (tid==0)
    {
        s.x = 0.0f;
        s.y = 0.0f;
        s.z = 0.0f;
        for (int i=0; i<numThreads; i++)
        {
            s.x += smema[i].x;
            s.y += smema[i].y;
            s.z += smema[i].z;
        }
        dInterData[bid].x = s.x;
        dInterData[bid].y = s.y;
        dInterData[bid].z = s.z;
    }
}


__global__ void gpuComponentSumReduce2(FLOAT_ARRAY_TYPE *dInterData, int numInterData, int numThreads)
{
    extern __shared__ FLOAT_ARRAY_TYPE smema[];
    int tid=threadIdx.x;
    FLOAT_ARRAY_TYPE s;
    s.x = 0.0f;
    s.y = 0.0f;
    s.z = 0.0f;

    for (int j=tid; j<numInterData; j+=numThreads)
    {
        s.x += dInterData[j].x;
        s.y += dInterData[j].y;
        s.z += dInterData[j].z;
    }

    smema[tid].x = s.x;
    smema[tid].y = s.y;
    smema[tid].z = s.z;
    __syncthreads();
    if (tid==0)
    {
        s.x = 0.0f;
        s.y = 0.0f;
        s.z = 0.0f;
        for (int i=0; i<numThreads; i++)
        {
            s.x += smema[i].x;
            s.y += smema[i].y;
            s.z += smema[i].z;
        }
        dInterData[0].x = s.x;
        dInterData[0].y = s.y;
        dInterData[0].z = s.z;
    }
}


void sumReduction(float *dInData, float *dInterData, int numInterData, int numThreads, int N)
{
    int smem = numThreads * sizeof(float);
    gpuReduce1<<<numInterData, numThreads, smem>>>(dInData, dInterData, numInterData, numThreads, N);
    checkCudaError("Kernel error: reduce 1");
    gpuReduce2<<<1, numThreads, smem>>>(dInterData, numInterData, numThreads);
    checkCudaError("Kernel error: reduce 2");
}



void SqVecSumReduction(FLOAT_ARRAY_TYPE *dInData, float *dInterData, int numInterData, int numThreads, int N)
{
    int smem = numThreads * sizeof(float);
    gpuSqSumReduce1<<<numInterData,numThreads, smem>>>(dInData, dInterData, numInterData, numThreads, N);
    checkCudaError("Kernel error: gpuSqSumReductionKernel");
    gpuReduce2<<<1,numThreads, smem>>>(dInterData, numInterData, numThreads);
    checkCudaError("Kernel error: reduce2");
}


void componentSumReduction(FLOAT_ARRAY_TYPE *dInData, FLOAT_ARRAY_TYPE *dInterData, int numInterData, int numThreads, int N)
{
    int smem = numThreads * sizeof(FLOAT_ARRAY_TYPE);
    gpuComponentSumReduce1<<<numInterData, numThreads, smem>>>(dInData, dInterData, numInterData, numThreads, N);
    checkCudaError("Kernel error: component reduction kernel");
    gpuComponentSumReduce2<<<1, numThreads, smem>>>(dInterData, numInterData, numThreads);
    checkCudaError("Kernel error: component reduce 2");
}

---