ver2

#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <shrUtils.h>
#include <cutil_inline.h>
#include "cutil_math.h"
#include "cutil_inline.h"
#include <math_functions.h>
#include <iostream>

#define uint32 unsigned int
#define int32  int


using namespace std;

size_t cmem_free;
size_t cmem_total;

#define BYTES 4


#define uchar unsigned char
#define BLOCK_SIZE 512
const long double DRO = 4.8481368122023329446e-6;
//float DRO = 4.8481368122023329446e-6;

#ifndef GEOD_CONSTANTS
#define GEOD_CONSTANTS
#define GEOD_PI  3.14159265358979323846
#define GEOD_PI2 1.57079632679489661923
#define GEOD_2PI 6.28318530717958647693
#define R_Z_GK 6367558.4968
#define R_Z_WGS 6367449.1457662908
#define A_WGS 6378137.
#define A_PZ90 6378136.
#define A_CK42 6378245.
#define F_WGS 298.257223563
#define F_PZ90 298.257839303
#define F_CK42 298.3
// e = 1/f
#define E_WGS 0.00335281066474748
#define E_PZ90 0.0033528037351843
#define E_CK42 0.00335232986925913
#endif // GEOD_CONSTANTS


__global__ void gdcTogdc(float *b1,float *l1, float *h1, float *b2,float *l2, float *h2)
{
    float val1 = (float)A_CK42;
    float val2 = (float)A_PZ90;
    float val3 = (float)E_CK42;
    float val4 = (float)E_PZ90;
    float val5 = 25.;
    float val6 = -141.;
    float val7 = -80.;
    float val8 = 0.;
    float val9 = -0.35;
    float val10 = -0.66;
    float val11 = 0.0;

    //int idx = threadIdx.x + blockIdx.x * blockDim.x;
    int idx = threadIdx.x;

    h1[idx] = 8;
    //h1[1] = 8;
    //h1[2] = 9;

 /*   float deltaA = (val1 - val2);
    float deltaE2 = (val3 - val4)*0.5;
    float a = (val1 + val2) *0.5;
    float eps2 = (val3 + val4) * 0.5;
    float es1 = - eps2  + 1.;
    float sinb = __sinf(b1[idx]);
    float cosb = __cosf(b1[idx]);
    float sinl = __sinf(l1[idx]);
    float cosl = __cosf(l1[idx]);
    float cos2b = __cosf(b1[idx]*2.);
    float sinBcosB = sinb * cosb;
    float sinB2 = sinb * sinb;
    float dXcosLplusdYsinL = val5 * cosl + val6 * sinl;
    float eps2sinB21 = - eps2 * sinB2 + 1.;
    float inFeps2sinB21 = 1./eps2sinB21;
    float temp1 = sqrtf(eps2sinB21);
    float N = a / temp1; // радиус кривизны первого вертикала
    float M = N * es1 * inFeps2sinB21; // радиус кривизны меридианного сечения
    float wycosLsubwxsinL = val9 * cosl - val8 * sinl;
    float temp = (N *sinBcosB* (deltaA * eps2  / a + (inFeps2sinB21 + 1.)* deltaE2) - dXcosLplusdYsinL * sinb + val7 * cosb) / (M + h1[idx]);
    temp -= eps2 * (val11)*sinBcosB;
    temp += wycosLsubwxsinL * (eps2 * cos2b + 1.);
    b2[idx] = b1[idx] + temp;
    l2[idx] = l1[idx] + (val6 * cosl - val5 * sinl) /((N + h1[idx]) * cosb) + tan(b1[idx]) * es1 * (val8 * cosl + val9 * sinl) - val10;
    temp = dXcosLplusdYsinL * cosb + val7 * sinb + (a*a / N + h1[idx]) * (val11);
    temp += N * deltaE2 * sinB2 - a * deltaA / N;
    temp += wycosLsubwxsinL * N * eps2 * sinBcosB;

    h2[idx] = h1[idx] + temp;*/
}

void gdc2gdc(float *b1, float *l1, float *h1, float *b2, float *l2, float *h2)
{
    //TSystemConvert convertSystem;
   // getSystemConvert(flag, convertSystem);

    const unsigned int numItems= 100;

    float * b1Dev, * l1Dev, *h1Dev;
    float * b2Dev, * l2Dev, *h2Dev;
    //float * dataDev;

    //float  data[11] = {val1,val2,val3,val4,val5,val6,val7,val8,val9,val10,val11};

/*
        cudaMalloc( (float**)&b1Dev, numItems*sizeof(float));
    cudaMalloc( (float**)&l1Dev, numItems*sizeof(float));
    cudaMalloc( (float**)&h1Dev, numItems*sizeof(float));
        cudaMalloc( (float**)&b2Dev, numItems*sizeof(float));
        cudaMalloc( (float**)&l2Dev, numItems*sizeof(float));
        cudaMalloc( (float**)&h2Dev, numItems*sizeof(float));
*/
        cutilSafeCall(cudaMalloc( (void**)&b1Dev, numItems*sizeof(float)));
    cutilSafeCall(cudaMalloc( (void**)&l1Dev, numItems*sizeof(float)));
    cutilSafeCall(cudaMalloc( (void**)&h1Dev, numItems*sizeof(float)));
        cutilSafeCall(cudaMalloc( (void**)&b2Dev, numItems*sizeof(float)));
        cutilSafeCall(cudaMalloc( (void**)&l2Dev, numItems*sizeof(float)));
        cutilSafeCall(cudaMalloc( (void**)&h2Dev, numItems*sizeof(float)));


  //    cudaMalloc( (void**)&dataDev, sizeof(data)*sizeof(float));

    dim3 threads = 512;
//dim3(512,1);
    dim3 blocks = 1;
//dim3(numItems/threads.x/4 + ((numItems/4 % threads.x) ? 1 : 0) ,1);
    cutilSafeCall( cudaThreadSynchronize() );

    cutilSafeCall(cudaMemcpy(b1Dev, b1, numItems*sizeof(float), cudaMemcpyHostToDevice));
    cutilSafeCall(cudaMemcpy(l1Dev, l1, numItems*sizeof(float), cudaMemcpyHostToDevice));
    cutilSafeCall(cudaMemcpy(h1Dev, h1, numItems*sizeof(float), cudaMemcpyHostToDevice));

    gdcTogdc <<<blocks,threads>>> (b1Dev,l1Dev,h1Dev,b2Dev,l2Dev,h2Dev);


    cutilSafeCall(cudaMemcpy((void*)b2,b2Dev,numItems*sizeof(float), cudaMemcpyDeviceToHost));
    cutilSafeCall(cudaMemcpy((void*)l2,l2Dev,numItems*sizeof(float), cudaMemcpyDeviceToHost));
    cutilSafeCall(cudaMemcpy((void*)h2,h2Dev,numItems*sizeof(float), cudaMemcpyDeviceToHost));

    cutilSafeCall( cudaThreadSynchronize() );
    cutilSafeCall(cudaFree(b1Dev));
    cutilSafeCall(cudaFree(l1Dev));
    cutilSafeCall(cudaFree(h1Dev));
    cutilSafeCall(cudaFree(b2Dev));
    cutilSafeCall(cudaFree(l2Dev));
    cutilSafeCall(cudaFree(h2Dev));
}

int main(int argc, char** argv)
{
 /*   float tr = l1/GEOD_2PI;
    l1 = l1 - GEOD_2PI * floor(tr);
    if (l1 < 0.)
    {
        l1 += GEOD_2PI;
    }
    else
        if (l1 > GEOD_2PI)
        {
            l1 -= GEOD_2PI;
        }
*/
    int arraySize = 100;

    float * b1 = new float[arraySize];
    float * l1 = new float[arraySize];
    float * h1 = new float[arraySize];

    for(int i=0; i<arraySize; i++)
    {
        *(b1+i) = i;
        *(l1+i) = i;
        *(h1+i) = i;
    }

    float  *b2 = new float[arraySize];
    float  *l2 = new float[arraySize];
    float  *h2 = new float[arraySize];

    for(int i = 0; i<100; i++)
    {
        cout<<"b1["<<i<<"]="<<b1[i]<<"\t";
        cout<<"l1["<<i<<"]="<<l1[i]<<"\t";
        cout<<"h1["<<i<<"]="<<h1[i]<<endl;
    }
        //printf("b1[%d] = %d\t l1[%d] = %d\t h1[%d] = %d\n",i,b1+i,i,l1+i,i,h1+i);

    printf("=====================================================================\n");

    gdc2gdc( b1, l1, h1, b2, l2, h2);

    for(int idx = 0; idx<100; idx++)
    {
        cout<<"b2["<<idx<<"]="<<b2[idx]<<"\t";
        cout<<"l2["<<idx<<"]="<<l2[idx]<<"\t";
        cout<<"h2["<<idx<<"]="<<h2[idx]<<endl;
    }

    return 0;
}