Why this CUDA program slower than a CPU program? How can I improve it's speed?

Question

I partially ported my lib (https://github.com/NicholasShatokhin/OpenNNL) in CUDA (https://github.com/NicholasShatokhin/OpenNNLCuda). But CUDA version slower than CPU version. Can anybody look function trainingBP and tell me what did I do wrong and how can I fix it? I'm novice in CUDA so I made many stupid errors. I have many memory transfers. How can I reduce their numbers?

Best regards, Nick

This is a piece of code:

__global__ void weighting(double * outputs, double * inputs, double * weights, int start, int inputsCount, int neuronsCount)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx < inputsCount * neuronsCount)
    {
        outputs[idx] = inputs[idx%inputsCount] * weights[start+idx];
    }
}

__global__ void calculateOutputsAndDerivatives(double * outputs, double * derivatives, double * layerOutputs, double * inputs, double * biases, int inputsCount, int neuronsCount, int neuronsInPreviousLayers)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx < neuronsCount)
    {
        double temp = 0;

        for(int k=0;k<inputsCount;k++)
        {
            temp += inputs[idx*inputsCount+k];
        }

        int neuron = neuronsInPreviousLayers + idx;

        temp -= biases[neuron];

        outputs[neuron] = layerOutputs[idx] = activation(temp);
        derivatives[neuron] = activation_derivative(temp);
    }
}

void OpenNNL::calculateNeuronsOutputsAndDerivatives(double * inputs, double * deviceOutputs, double * deviceDerivatives)
{
    int inputsCount = _inputsCount;

    double * deviceTemp;
    double * deviceInputs;

    cudaCall(cudaMalloc ( (void**)&deviceInputs, inputsCount*sizeof(double) ));

    cudaCall(cudaMemcpy ( deviceInputs, inputs, inputsCount*sizeof(double), cudaMemcpyDeviceToDevice ));

    for(int i=0;i<_layersCount;i++)
    {
        cudaCall(cudaMalloc((void**)&deviceTemp, _neuronsPerLayerCount[i]*inputsCount*sizeof(double)));

        dim3 threadsMul = dim3(BLOCK_SIZE, 1);
        int blocksCount = floor((double) _neuronsPerLayerCount[i]*inputsCount / threadsMul.x) + 1;
        dim3 blocksMul = dim3(blocksCount, 1);

        weighting<<<blocksMul, threadsMul>>>(deviceTemp, deviceInputs, _neuronsInputsWeights, _inputsInPreviousLayers[i], inputsCount, _neuronsPerLayerCount[i]);

        cudaCall(cudaFree(deviceInputs));

        cudaCall(cudaMalloc((void**)&deviceInputs, _neuronsPerLayerCount[i]*sizeof(double)));

        dim3 threadsSum = dim3(BLOCK_SIZE, 1);
        blocksCount = floor((double) _neuronsPerLayerCount[i] / threadsSum.x) + 1;
        dim3 blocksSum = dim3(blocksCount, 1);

        calculateOutputsAndDerivatives<<<blocksSum, threadsSum>>>(deviceOutputs, deviceDerivatives, deviceInputs, deviceTemp, _neuronsBiases, inputsCount, _neuronsPerLayerCount[i], _neuronsInPreviousLayers[i]);

        inputsCount = _neuronsPerLayerCount[i];

        cudaCall(cudaFree(deviceTemp));
    }

    cudaCall(cudaFree(deviceInputs));
}

__global__ void calculateLocalGradientsForLastLayer(double * localGradients, double * error, double * outputs, double * derivatives, double * trainingOutputs, double sample_weight, int neuronsCount, int neuronsInLastLayers)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx < neuronsCount)
    {
        double current_error = trainingOutputs[idx] - outputs[neuronsInLastLayers + idx];
        localGradients[neuronsInLastLayers + idx] = current_error * sample_weight * derivatives[neuronsInLastLayers + idx];

        error[idx] = current_error * current_error;
    }
}

__global__ void calculateLocalGradientsForAnotherLayers(double * localGradients, double * neuronsInputsWeights, double * derivatives, int neuronsCount, int neuronsInPreviousLayers, int neuronsInPreviousLayersWithCurrent, int neuronsInNextLayer, int inputsInPreviousLayers, int inputsInCurrentLayer)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx < neuronsCount)
    {
        int neuron = neuronsInPreviousLayers + idx;

        localGradients[neuron] = 0;

        // this to Kernel, then reduce localGradients.
        for(int k=0;k<neuronsInNextLayer;k++)
        {
            localGradients[neuron] += neuronsInputsWeights[inputsInPreviousLayers + k*inputsInCurrentLayer + idx]
                                                            * localGradients[neuronsInPreviousLayersWithCurrent + k];
        }

        localGradients[neuron] *= derivatives[neuron];
    }
}

__global__ void changeWeightsForFirstLayer(double * neuronsInputsWeights, double * trainingInputs, double * localGradients, double speed, int inputsInLayer, int inputsCount)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx < inputsInLayer)
    {
        double deltaWeight = speed * localGradients[(int) (idx / inputsCount)] * trainingInputs[idx % inputsCount];
        double temp = neuronsInputsWeights[idx] + deltaWeight;
        neuronsInputsWeights[idx] = temp;
    }
}

__global__ void changeBiasesForFirstLayer(double * neuronsBiases, double * localGradients, double speed, int neuronsCount)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx < neuronsCount)
    {
        double deltaBias = speed * localGradients[idx];
        double temp = neuronsBiases[idx] - deltaBias;
        neuronsBiases[idx] = temp;
    }
}

__global__ void changeWeightsForAnotherLayers(double * neuronsInputsWeights, double * localGradients, double * outputs, int * neuronsInPreviousLayers, int * inputsInPreviousLayers, int * inputsInCurrentLayer, double speed, int inputsCount, int layersCount)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx >= inputsInPreviousLayers[1] && idx < inputsCount)
    {
        // searching layer's number. ugly method -- TODO: write better
        int layer = layersCount-1;
        for(int i=0;i<layersCount-1;i++)
        {
            if(idx >= inputsInPreviousLayers[i] && idx < inputsInPreviousLayers[i+1])
            {
                layer = i;
                break;
            }
        }

        int index1 = neuronsInPreviousLayers[layer] + ((int) ((idx - inputsInPreviousLayers[layer]) / inputsInCurrentLayer[layer]));
        int index2 = neuronsInPreviousLayers[layer-1] + ((idx - inputsInPreviousLayers[layer]) % inputsInCurrentLayer[layer]);
        neuronsInputsWeights[idx] += speed * localGradients[index1] * outputs[index2];
    }
}

__global__ void changeBiasesForAnotherLayers(double * neuronsBiases, double * localGradients, double speed, int neuronsCount, int neuronsInFirstLayer)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if(idx >= neuronsInFirstLayer && idx < neuronsCount)
    {
        neuronsBiases[idx] -= speed * localGradients[idx];
    }
}

__global__ void reduceSum ( double * inData, double * outData )
{
    __shared__ double data [BLOCK_SIZE];
    int tid = threadIdx.x;
    int i = 2 * blockIdx.x * blockDim.x + threadIdx.x;

    data [tid] = inData [i] + inData [i+blockDim.x]; // load into shared memeory

    __syncthreads ();

    for ( int s = blockDim.x / 2; s > 32; s >>= 1 )
    {
        if ( tid < s )
            data [tid] += data [tid + s];

        __syncthreads ();
    }

    if ( tid < 32 ) // unroll last iterations
    {
        data [tid] += data [tid + 32];
        data [tid] += data [tid + 16];
        data [tid] += data [tid + 8];
        data [tid] += data [tid + 4];
        data [tid] += data [tid + 2];
        data [tid] += data [tid + 1];
    }

    if ( tid == 0 ) // write result of block reduction
        outData [blockIdx.x] = data [0];
}

double OpenNNL::_changeWeightsByBP(double * deviceTrainingInputs, double * deviceTrainingOutputs, double speed, double sample_weight)
{
    dim3 threads = dim3(BLOCK_SIZE, 1);
    int blocksCount;
    dim3 blocks;

    double error = 0, current_error;
    int neuronsInLastLayer = _neuronsPerLayerCount[_layersCount-1];
    int index = 0, n = neuronsInLastLayer;
    double * errors = new double[neuronsInLastLayer];

    double * deviceOutputs;
    double * deviceDerivatives;

    double * deviceLocalGradients;
    double * deviceErrors[2] = {NULL, NULL};

    cudaCall(cudaMalloc ( (void**)&deviceOutputs, _neuronsCount*sizeof(double) ));
    cudaCall(cudaMalloc ( (void**)&deviceDerivatives, _neuronsCount*sizeof(double) ));

    cudaCall(cudaMalloc ( (void**)&deviceLocalGradients, _neuronsCount*sizeof(double) ));
    cudaCall(cudaMalloc ( (void**)&deviceErrors[0], neuronsInLastLayer*sizeof(double) ));
    cudaCall(cudaMalloc ( (void**)&deviceErrors[1], neuronsInLastLayer*sizeof(double) ));

    calculateNeuronsOutputsAndDerivatives(deviceTrainingInputs, deviceOutputs, deviceDerivatives);

    blocksCount = floor((double) neuronsInLastLayer / threads.x) + 1;
    blocks = dim3(blocksCount, 1);

    calculateLocalGradientsForLastLayer <<<blocks, threads>>>(deviceLocalGradients, deviceErrors[0], deviceOutputs, deviceDerivatives, deviceTrainingOutputs, sample_weight, neuronsInLastLayer, _neuronsInPreviousLayers[_layersCount-1]);

    for (index = 0; n >= BLOCK_SIZE; n /= (2*BLOCK_SIZE), index ^= 1 )
    {
                               // set kernel launch configuration
        dim3 dimBlock ( BLOCK_SIZE, 1, 1 );
        dim3 dimGrid ( n / (2*dimBlock.x), 1, 1 );

        reduceSum <<<dimGrid, dimBlock>>> ( deviceErrors [index], deviceErrors [index^1] );
    }

    cudaMemcpy ( errors, deviceErrors [index], neuronsInLastLayer*sizeof(double), cudaMemcpyDeviceToHost );

    for ( index = 1; index < n; index++ )
            errors [0] += errors [index];

    if(_layersCount > 1)
    {
        for(int i=_layersCount-2;i>=0;i--)
        {
            // calculateLocalGradientsForAnotherLayers
            blocksCount = floor((double) _neuronsPerLayerCount[i] / threads.x) + 1;
            blocks = dim3(blocksCount, 1);

            calculateLocalGradientsForAnotherLayers <<<blocks, threads>>> (deviceLocalGradients, _neuronsInputsWeights, deviceDerivatives, _neuronsPerLayerCount[i], _neuronsInPreviousLayers[i], _neuronsInPreviousLayers[i+1], _neuronsPerLayerCount[i+1], _inputsInPreviousLayers[i], _inputsInCurrentLayer[i]);
        }
    }

    blocksCount = floor((double) _neuronsPerLayerCount[0] * _inputsCount / threads.x) + 1;
    blocks = dim3(blocksCount, 1);

    changeWeightsForFirstLayer <<<blocks, threads>>> (_neuronsInputsWeights, deviceTrainingInputs, deviceLocalGradients, speed, _neuronsPerLayerCount[0] * _inputsCount, _inputsCount);

    blocksCount = floor((double) _neuronsPerLayerCount[0] / threads.x) + 1;
    blocks = dim3(blocksCount, 1);

    changeBiasesForFirstLayer <<<blocks, threads>>> (_neuronsBiases, deviceLocalGradients, speed, _neuronsPerLayerCount[0]);

    blocksCount = floor((double) _weightsCount / threads.x) + 1;
    blocks = dim3(blocksCount, 1);

    changeWeightsForAnotherLayers <<<blocks, threads>>> (_neuronsInputsWeights, deviceLocalGradients, deviceOutputs, _deviceNeuronsInPreviousLayers, _deviceInputsInPreviousLayers, _deviceInputsInCurrentLayer, speed, _weightsCount, _layersCount);

    blocksCount = floor((double) _neuronsCount / threads.x) + 1;
    blocks = dim3(blocksCount, 1);

    changeBiasesForAnotherLayers <<<blocks, threads>>> (_neuronsBiases, deviceLocalGradients, speed, _neuronsCount, _neuronsInPreviousLayers[1]);

    cudaCall(cudaFree(deviceLocalGradients));
    cudaCall(cudaFree(deviceOutputs));
    cudaCall(cudaFree(deviceDerivatives));

    cudaCall(cudaFree(deviceErrors[0]));
    cudaCall(cudaFree(deviceErrors[1]));

    error = errors[0]/2;
    return error;
}