an alpha kernel for step integration of HH neuronal model

AlphaKernel.cl

    // OpenCL Kernel Function for Hodgkin Huxley integration step
    kernel void IntegrateHHStep(const float maxG_K,
    							const float maxG_Na,
    							const float maxG_Leak,
    							const float E_K,
    							const float E_Na,
    							const float E_Leak,
    							const float I_ext,
    							const float dt,
    							global const float* V_in, 
    							global const float* x_n_in,
    							global const float* x_m_in,
    							global const float* x_h_in,
    							global float* V_out,
    							global float* x_n_out,
    							global float* x_m_out,
    							global float* x_h_out, 
    							int numElements) {
        // get index into global data array
        int iGID = get_global_id(0);

        // bound check (equivalent to the limit on a 'for' loop for standard/serial C code
        if (iGID >= numElements)  {
            return;
        }
            
		// logic for step integration
        // alpha functions
		float alpha[3];
		alpha[0] = (10 - V_in[iGID]) / (100 * (exp((10 - V_in[iGID]) / 10) - 1));
		alpha[1] = (25 - V_in[iGID]) / (10 * (exp((25 - V_in[iGID]) / 10) - 1));
		alpha[2] = 0.07 * exp(-V_in[iGID] / 20);
		// beta functions
		float beta[3];
		beta[0] = 0.125 * exp(-V_in[iGID] / 80);
		beta[1] = 4 * exp(-V_in[iGID] / 18);
		beta[2] = 1 / (exp((30 - V_in[iGID]) / 10) + 1);

		// calculate tau and x0 with alpha and beta
		float tau[3];
		for (int i = 0; i < 3; i++) {
			tau[i] = 1 / (alpha[i] + beta[i]);
		}

		float x0[3];
		for (int i = 0; i < 3; i++) {
			x0[i] = alpha[i] * tau[i];
		}

		// leaky integration for Xs with eurler's method
		x_n_out[iGID] = (1 - dt / tau[0]) * x_n_in[iGID] + dt / tau[0] * x0[0];
		x_m_out[iGID] = (1 - dt / tau[1]) * x_m_in[iGID] + dt / tau[1] * x0[1];
		x_h_out[iGID] = (1 - dt / tau[2]) * x_h_in[iGID] + dt / tau[2] * x0[2];

		// calculate conductances for n, m, h
		float gnmh[3];
		gnmh[0] = maxG_K * pow(x_n_out[iGID], 4);
		gnmh[1] = maxG_Na * pow(x_m_out[iGID], 3) * x_h_out[iGID];
		gnmh[2] = maxG_Leak;

		// calculate current with Ohm's law
		float I[3];
		I[0] = gnmh[0] * (V_in[iGID] - E_K);
		I[1] = gnmh[1] * (V_in[iGID] - E_Na);
		I[2] = gnmh[2] * (V_in[iGID] - E_Leak);

		// given all the currents, update voltage membrane
		V_out[iGID] = V_in[iGID] + dt * (I_ext - (I[0] + I[1] + I[2]));
    }