Cuda

Что такое WSL?

CUDA в WSL

NVML

GPU-контейнеры в WSL

Запуск контейнера N-body

user@PCName:/mnt/c$ curl https://get.docker.com | sh

user@PCName:/mnt/c$ distribution=$(. /etc/os-release;echo $ID$VERSION_ID)user@PCName:/mnt/c$ curl -s -L https://nvidia.github.io/nvidia-docker/gpgkey | sudo apt-key add -user@PCName:/mnt/c$ curl -s -L https://nvidia.github.io/nvidia-docker/$distribution/nvidia-docker.list | sudo tee /etc/apt/sources.list.d/nvidia-docker.listuser@PCName:/mnt/c$ curl -s -L https://nvidia.github.io/libnvidia-container/experimental/$distribution/libnvidia-container-experimental.list | sudo tee /etc/apt/sources.list.d/libnvidia-container-experimental.list

user@PCName:/mnt/c$ sudo apt-get updateuser@PCName:/mnt/c$ sudo apt-get install -y nvidia-docker2

user@PCName:/mnt/c$ sudo dockerd

user@PCName:/mnt/c$ docker run --gpus all nvcr.io/nvidia/k8s/cuda-sample:nbody nbody -gpu -benchmark

Запуск контейнера TensorFlow

user@PCName:/mnt/c$ docker pull tensorflow/tensorflow:latest-gpu-py3

user@PCName:/mnt/c$ vi ./matmul.pyimport sysimport numpy as npimport tensorflow as tffrom datetime import datetimedevice_name = sys.argv[1] # Choose device from cmd line. Options: gpu or cpushape = (int(sys.argv[2]), int(sys.argv[2]))if device_name == "gpu":device_name = "/gpu:0"else:device_name = "/cpu:0"tf.compat.v1.disable_eager_execution()with tf.device(device_name):random_matrix = tf.random.uniform(shape=shape, minval=0, maxval=1)dot_operation = tf.matmul(random_matrix, tf.transpose(random_matrix))sum_operation = tf.reduce_sum(dot_operation)startTime = datetime.now()with tf.compat.v1.Session(config=tf.compat.v1.ConfigProto(log_device_placement=True)) as session:result = session.run(sum_operation)print(result)# Вывод результатовprint("Shape:", shape, "Device:", device_name)print("Time taken:", datetime.now() - startTime)

user@PCName:/mnt/c$ docker run --runtime=nvidia --rm -ti -v "${PWD}:/mnt/c" tensorflow/tensorflow:latest-gpu-jupyter python /mnt/c/matmul.py gpu 20000

user@PCName:/mnt/c$ docker run -it --gpus all -p 8888:8888 tensorflow/tensorflow:latest-gpu-py3-jupyter

Обзор WSL

Linux-ядро WSL 2

Поддержка GPU в WSL

Основные сведения о GPU-PV

• Windows Sandbox
• Microsoft Defender Application Guard
• Эмулятор Microsoft HoloLens 2

Поддержка пользовательского режима

Попробуйте новые возможности WSL 2

• Для того чтобы получить доступ к Preview-версии WSL 2 с поддержкой GPU-ускорения, вы можете присоединиться к Windows Insider Program.
• Загрузите свежие драйверы NVIDIA, установите их и попробуйте запустить в WSL 2 какой-нибудь CUDA-контейнер.

В последнее время коллеги по "цеху" независимо друг от друга стали спрашивать меня: как получить c одного SDR-приемника одновременно все каналы Bluetooth? Полоса ведь позволяет, есть SDR с выходной полосой 80 МГц и более. Можно, конечно, сделать это на ПЛИС, но время такой разработки будет довольно большим. Мне давно было известно, что сделать такое на GPU довольно просто, но чтобы так!

Стандарт Bluetooth определяет физический уровень в двух версиях: Classic и Low Energy. Спецификация есть здесь. Документ ужасно большой, читать его целиком опасно для мозга. К счастью, большие компании, производящие измерительную технику, имеют средства для создания наглядных документов по теме. Tektronix и National Instruments, например. У меня совершенно нет шансов в конкуренции с ними по качеству представления материала. Интересующихся прошу изучать по ссылкам.

Все, что мне нужно знать о физическом уровне для создания многоканального фильтра, это шаг сетки частот и скорость модуляции. Они сведены в таблицу в одном из указанных документов:

Таким образом, нам нужно нарезать полосу 80 МГц на 79 фильтров с шагом настройки 1 МГц и, одновременно с этим, на 40 фильтров с шагом настройки 2 МГц. Частоты следования отсчетов с выходов фильтров должны быть 1 МГц и 2 МГц, соответственно.

Таким образом, нам необходимо две гребенки фильтров.

Для начала выберем параметры этих фильтров исходя из полос сигналов Bluetooth Classic и Bluetooth Low Energy. Нам нужны их импульсные характеристики, чтобы рассчитать нагрузку на вычислительное устройство фильтра. Здесь сразу стоит оговориться, что длины импульсных характеристик мы выбрали исходя из требований "быстрого" алгоритма фильтрации. Суть от этого не меняется. И число коэффициентов импульсной характеристики не должно быть слишком большим, чтобы фильтр был реализуем на вменяемой вычислительной аппаратуре.

Для фильтров с шагом 1 МГц выберем полосу пропускания ФНЧ (половина полосы пропускания полосового фильтра) 500 кГц, длину импульсной характеристики подгоним к 480 отводам. Для фильтров с шагом 2 МГц эти параметры выберем 1 МГц и 240 отводов, соответственно. Тип окна выбираем Кайзера. Рассчитаем импульсные характеристики в filterDesigner и выгрузим их в формате С-header:

Можно решать задачу лобовым способом: строить соответствующий по числу фильтров массив DDC (Digital Down Converter). Такой подход хорош для ПЛИС, где возможна экономия за счет уменьшения разрядности вычислителей первых ступеней. Также ПЛИС это наиболее энергетически-эффективный способ реализации. Но трудозатраты при этом способе наиболее высоки.

При выполнении гребенки фильтров на популярных нынче GPU появляется возможность реализации более хитрого алгоритма: полифазной гребенки фильтров на основе БПФ, который на CUDA доступен из библиотеки. В заграничной литературе алгоритм называется Polyphase or WOLA (Weight, Overlap and Add) FFT Filterbank. Лень к рисованию не дает мне возможности выполнить наглядное объяснение самостоятельно. В Сети есть много материалов по теме, особенно наглядный график выполнен здесь на странице 11 (огромное спасибо уважаемым авторам), вот он:

Железо выберем то, которое есть под рукой. Это плата ввода компании "Инструментальные Системы" FMC126P. О ней я уже писал в одной предыдущей статье. В разъем FMC платы вставлен субмодуль той же компании с трансивером AD9371 с полосой 100 МГц. Весь поток с трансивера может быть непрерывно передан в компьютер для обработки.

Выберем видео-карту с GPU GTX 1050. (Соврал, это она нас выбрала: это все что было под рукой, была выдрана из вычислителя для расчета антенн, но тем удивительней было увидеть рабочую гребенку). Перейдем к программной части.

К сожалению, из-за лицензий мы не можем опубликовать полный код. Можем показать только ядра GPU. Впрочем, остальной код и не особо интересен.

Вот ядро, которое выполняет перемножение сигнала на окно и сложение, и обертка для его вызова:

__global__ void cuComplexMultiplyWindowKernel(const cuComplex *data, const float *window, size_t windowSize, cuComplex *result) {    __shared__ cuComplex multiplicationResult[480];    multiplicationResult[threadIdx.x] = cuComplexMultiplyFloat(data[threadIdx.x + windowSize / 4 * blockIdx.x], window[threadIdx.x]);    __syncthreads();    cuComplex sum;    sum.x = sum.y = 0;    if (threadIdx.x < windowSize / 4) {        for(int i = 0; i < 4; i++) {            sum = cuComplexAdd(sum, multiplicationResult[threadIdx.x + i * windowSize / 4]);        }        result[threadIdx.x + windowSize / 4 * blockIdx.x] = sum;    }}cudaError_t cuComplexMultiplyWindow(const cuComplex *data, const float *window, size_t windowSize, cuComplex *result, size_t dataSize, cudaStream_t stream) {    size_t windowStep = windowSize / 4;    cuComplexMultiplyWindowKernel<<<dataSize / windowStep - 3, windowSize, 1024, stream>>>(data, window, windowSize, result);    return cudaGetLastError();}

Код обработки сигнала, который вызывает это ядро, в точности повторяет схему алгоритма, изображенную на рисунке выше, поэтому не вижу смысла приводить его здесь.

Гребенки были проверены на выходном спектре каналов в реальном времени. На вход AD9371 подавался сигнал генератора сигналов 2450 МГц, селективность фильтров соответствовала расчетной.

В планах: адаптация софта на плату XRTX и реализация поиска пакетов, если это кому-нибудь окажется нужно или будет свободное время.

Всю работу по софту выполнил gaudima, слава ему!

void ewald_rec(){    int mmin = 0;    int nmin = 1;    // массивы где хранятся iexp(x[i] * kx[l]),    double** elc;    double** els;    //... iexp(y[i] * ky[m]) и    double** emc;    double** ems;    //... iexp(z[i] * kz[n]),    double** enc;    double** ens;    // временные массивы для произведений iexp(x*kx)*iexp(y*ky)    double* lmc;    double* lms;    // и для q[i] * iexp(x*kx)*iexp(y*ky)*iexp(z*kz)    double* ckc;    double* cks;    // ПРЕДВАРИТЕЛЬНОЕ ЗАПОЛНЕНИЕ МАССИВОВ    eng = 0.0;    for (i = 0; i < Nat; i++)   // цикл по атомам    {        // emc/s[i][0] и enc/s[i][0] уже заполнены на этапе инициализации        // в массив elc/s нужно обновить, см. далее        elc[i][0] = 1.0;        els[i][0] = 0.0;        // iexp(kr)        sincos(twopi * xs[i] * ra, els[i][1], elc[i][1]);        sincos(twopi * ys[i] * rb, ems[i][1], emc[i][1]);        sincos(twopi * zs[i] * rc, ens[i][1], enc[i][1]);    }    // заполняем следующие элементы массива emc/s[i][l] = iexp(y[i]*ky[l]) итеративно, используя комплексное умножение    for (l = 2; l < ky; l++)        for (i = 0; i < Nat; i++)        {            emc[i][l] = emc[i][l - 1] * emc[i][1] - ems[i][l - 1] * ems[i][1];            ems[i][l] = ems[i][l - 1] * emc[i][1] + emc[i][l - 1] * ems[i][1];        }    // заполняем следующие элементы массива enc/s[i][l] = iexp(z[i]*kz[l]) итеративно, используя комплексное умножение    for (l = 2; l < kz; l++)        for (i = 0; i < Nat; i++)        {            enc[i][l] = enc[i][l - 1] * enc[i][1] - ens[i][l - 1] * ens[i][1];            ens[i][l] = ens[i][l - 1] * enc[i][1] + enc[i][l - 1] * ens[i][1];        }    // ГЛАВНЙ ЦИКЛ ПО ВСЕМ K-ВЕКТОРАМ:    for (l = 0; l < kx; l++)    {        rkx = l * twopi * ra;         // записываем exp(ikx[l]) в ikx[0] для сохранения памяти        if (l == 1)            for (i = 0; i < Nat; i++)            {                elc[i][0] = elc[i][1];                els[i][0] = els[i][1];            }        else if (l > 1)            for (i = 0; i < Nat; i++)            {                // iexp(kx[0]) = iexp(kx[0]) * iexp(kx[1])                x = elc[i][0];                elc[i][0] = x * elc[i][1] - els[i][0] * els[i][1];                els[i][0] = els[i][0] * elc[i][1] + x * els[i][1];            }        for (m = mmin; m < ky; m++)        {            rky = m * twopi * rb;            // заполняем временный массив произведением iexp(kx*x[i]) * iexp(ky*y[i])            if (m >= 0)                for (i = 0; i < Nat; i++)                {                    lmc[i] = elc[i][0] * emc[i][m] - els[i][0] * ems[i][m];                    lms[i] = els[i][0] * emc[i][m] + ems[i][m] * elc[i][0];                }            else // для отрицательных значений m используем комплексное сопряжение:                for (i = 0; i < Nat; i++)                {                    lmc[i] = elc[i][0] * emc[i][-m] + els[i][0] * ems[i][-m];                    lms[i] = els[i][0] * emc[i][-m] - ems[i][-m] * elc[i][0];                }            for (n = nmin; n < kz; n++)            {                rkz = n * twopi * rc;                rk2 = rkx * rkx + rky * rky + rkz * rkz;                if (rk2 < rkcut2) // используем радиус обрезания                {                    // вычисляем сумму (q[i]*iexp(kr[k]*r[i])) - зарядовую плотность                    sumC = 0; sumS = 0;                    if (n >= 0)                        for (i = 0; i < Nat; i++)                        {                            //считываем заряд частицы                            ch = charges[types[i]].charge;                            ckc[i] = ch * (lmc[i] * enc[i][n] - lms[i] * ens[i][n]);                            cks[i] = ch * (lms[i] * enc[i][n] + lmc[i] * ens[i][n]);                            sumC += ckc[i];                            sumS += cks[i];                        }                    else // для отрицательных индексов используем комплексное сопряжение:                        for (i = 0; i < Nat; i++)                        {                            //считываем заряд частицы                            ch = charges[types[i]].charge;                            ckc[i] = ch * (lmc[i] * enc[i][-n] + lms[i] * ens[i][-n]);                            cks[i] = ch * (lms[i] * enc[i][-n] - lmc[i] * ens[i][-n]);                            sumC += ckc[i];                            sumS += cks[i];                        }                    //наконец вычисляем энергию и силы                    akk = exp(rk2 * elec->mr4a2) / rk2;                    eng += akk * (sumC * sumC + sumS * sumS);                    for (i = 0; i < Nat; i++)                    {                        x = akk * (cks[i] * sumC - ckc[i] * sumS) * C * twopi * 2 * rvol;                        fxs[i] += rkx * x;                        fys[i] += rky * x;                        fzs[i] += rkz * x;                    }                }            } // end n-loop (over kz-vectors)            nmin = 1 - kz;        } // end m-loop (over ky-vectors)        mmin = 1 - ky;    }  // end l-loop (over kx-vectors)   engElec2 = eng * С * twopi * rvol;}

__global__ void recip_ewald(int atPerBlock, int atPerThread, cudaMD* md)// calculate reciprocal part of Ewald summ// the first part : summ (qiexp(kr)) evaluation{    int i;      // for atom loop    int ik;     // index of k-vector    int l, m, n;    int mmin = 0;    int nmin = 1;    float tmp, ch;    float rkx, rky, rkz, rk2;   // component of rk-vectors    int nkx = md->nk.x;    int nky = md->nk.y;    int nkz = md->nk.z;        // arrays for keeping iexp(k*r) Re and Im part    float2 el[2];    float2 em[NKVEC_MX];    float2 en[NKVEC_MX];    float2 sums[NTOTKVEC];          // summ (q iexp (k*r)) for each k    extern __shared__ float2 sh_sums[];     // the same in shared memory    float2 lm;     // temp var for keeping el*em    float2 ck;     // temp var for keeping q * el * em * en (q iexp (kr))    // invert length of box cell    float ra = md->revLeng.x;    float rb = md->revLeng.y;    float rc = md->revLeng.z;    if (threadIdx.x == 0)        for (i = 0; i < md->nKvec; i++)            sh_sums[i] = make_float2(0.0f, 0.0f);    __syncthreads();    for (i = 0; i < md->nKvec; i++)        sums[i] = make_float2(0.0f, 0.0f);    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    ik = 0;    for (i = id0; i < N; i++)    {        // save charge        ch = md->specs[md->types[i]].charge;        el[0] = make_float2(1.0f, 0.0f);    // .x - real part (or cos) .y - imagine part (or sin)                em[0] = make_float2(1.0f, 0.0f);        en[0] = make_float2(1.0f, 0.0f);        // iexp (ikr)        sincos(d_2pi * md->xyz[i].x * ra, &(el[1].y), &(el[1].x));        sincos(d_2pi * md->xyz[i].y * rb, &(em[1].y), &(em[1].x));        sincos(d_2pi * md->xyz[i].z * rc, &(en[1].y), &(en[1].x));        // fil exp(iky) array by complex multiplication        for (l = 2; l < nky; l++)        {             em[l].x = em[l - 1].x * em[1].x - em[l - 1].y * em[1].y;             em[l].y = em[l - 1].y * em[1].x + em[l - 1].x * em[1].y;        }        // fil exp(ikz) array by complex multiplication        for (l = 2; l < nkz; l++)        {             en[l].x = en[l - 1].x * en[1].x - en[l - 1].y * en[1].y;             en[l].y = en[l - 1].y * en[1].x + en[l - 1].x * en[1].y;        }        // MAIN LOOP OVER K-VECTORS:        for (l = 0; l < nkx; l++)        {            rkx = l * d_2pi * ra;                         // move exp(ikx[l]) to ikx[0] for memory saving (ikx[i>1] are not used)            if (l == 1)                el[0] = el[1];            else if (l > 1)                {                    // exp(ikx[0]) = exp(ikx[0]) * exp(ikx[1])                    tmp = el[0].x;                    el[0].x = tmp * el[1].x - el[0].y * el[1].y;                    el[0].y = el[0].y * el[1].x + tmp * el[1].y;                }            //ky - loop:            for (m = mmin; m < nky; m++)            {                rky = m * d_2pi * rb;                //set temporary variable lm = e^ikx * e^iky                if (m >= 0)                {                        lm.x = el[0].x * em[m].x - el[0].y * em[m].y;                         lm.y = el[0].y * em[m].x + em[m].y * el[0].x;                }                else // for negative ky give complex adjustment to positive ky:                {                        lm.x = el[0].x * em[-m].x + el[0].y * em[-m].y;                        lm.y = el[0].y * em[-m].x - em[-m].x * el[0].x;                }                //kz - loop:                for (n = nmin; n < nkz; n++)                {                    rkz = n * d_2pi * rc;                    rk2 = rkx * rkx + rky * rky + rkz * rkz;                    if (rk2 < md->rKcut2) // cutoff                    {                        // calculate summ[q iexp(kr)]   (local part)                        if (n >= 0)                         {                                ck.x = ch * (lm.x * en[n].x - lm.y * en[n].y);                                ck.y = ch * (lm.y * en[n].x + lm.x * en[n].y);                         }                        else // for negative kz give complex adjustment to positive kz:                         {                                ck.x = ch * (lm.x * en[-n].x + lm.y * en[-n].y);                                ck.y = ch * (lm.y * en[-n].x - lm.x * en[-n].y);                        }                        sums[ik].x += ck.x;                        sums[ik].y += ck.y;                                                // save qiexp(kr) for each k for each atom:                        md->qiexp[i][ik] = ck;                        ik++;                    }                } // end n-loop (over kz-vectors)                nmin = 1 - nkz;            } // end m-loop (over ky-vectors)            mmin = 1 - nky;        }  // end l-loop (over kx-vectors)    } // end loop by atoms    // save sum into shared memory    for (i = 0; i < md->nKvec; i++)    {        atomicAdd(&(sh_sums[i].x), sums[i].x);        atomicAdd(&(sh_sums[i].y), sums[i].y);    }    __syncthreads();    //...and to global    int step = ceil((double)md->nKvec / (double)blockDim.x);    id0 = threadIdx.x * step;    N = min(id0 + step, md->nKvec);    for (i = id0; i < N; i++)    {        atomicAdd(&(md->qDens[i].x), sh_sums[i].x);        atomicAdd(&(md->qDens[i].y), sh_sums[i].y);    }}// end 'ewald_rec' function__global__ void ewald_force(int atPerBlock, int atPerThread, cudaMD* md)// calculate reciprocal part of Ewald summ// the second part : enegy and forces{    int i;      // for atom loop    int ik;     // index of k-vector    float tmp;    // accumulator for force components    float3 force;    // constant factors for energy and force    float eScale = md->ewEscale;    float fScale = md->ewFscale;    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    for (i = id0; i < N; i++)    {        force = make_float3(0.0f, 0.0f, 0.0f);        // summ by k-vectors        for (ik = 0; ik < md->nKvec; ik++)        {            tmp = fScale * md->exprk2[ik] * (md->qiexp[i][ik].y * md->qDens[ik].x - md->qiexp[i][ik].x * md->qDens[ik].y);            force.x += tmp * md->rk[ik].x;            force.y += tmp * md->rk[ik].y;            force.z += tmp * md->rk[ik].z;        }        md->frs[i].x += force.x;        md->frs[i].y += force.y;        md->frs[i].z += force.z;    } // end loop by atoms    // one block calculate energy    if (blockIdx.x == 0)        if (threadIdx.x == 0)        {            for (ik = 0; ik < md->nKvec; ik++)                md->engCoul2 += eScale * md->exprk2[ik] * (md->qDens[ik].x * md->qDens[ik].x + md->qDens[ik].y * md->qDens[ik].y);        }}// end 'ewald_force' function

Вместо заключения

    //bonds:    int nBond;      // number of bonds    int mxBond;          // maximal number of bonds    int4* bonds;    // array of bonds     int* nbonds;    // count of bond for a given atom    int* neighToBind;   // a neighbor of a given atom for binding    int* canBind;       // flags that atom[iat] can be bind    int* r2Min;         // distances for the nearest neighbor (used for binding)    int* parents;       // indexes of one of the atom bonded with a given    cudaBond* bondTypes;     int** def_bonds;    // array[nSpec][nSpec] of default bond types    int** bindBonds;    // array[nSpec][nSpec] bond types created by binding    float** bindR2;        // square of binding distance [nSpec][nSpec]    //angles:    int nAngle;    // number of angles    int mxAngle;    int4* angles;   // array of angles      int* nangles;        // number of angles for given atom    int* oldTypes;          cudaAngle* angleTypes;     int* specAngles;    // [nSp] angle type formed by given species

struct cudaBond{    int type;  // potential type    int spec1, spec2; // type of atoms that connected by this bond type    int new_type[2];      // bond type after mutation    int new_spec1[2], new_spec2[2];    int mxEx, mnEx;     // flags: maximum or minimum of bond length exists    float p0, p1, p2, p3, p4;    // potential parameters    float r2min, r2max;          // square of minimal and maximal bond length    float (*force_eng)(float r2, float r, float &eng, cudaBond *bond); // return energy     int count;      // quantity of such bonds    float rSumm;        // summ of lentghs (for mean length calculation)    int rCount;          // number of measured lengths (for mean length calculation)    int ltSumm, ltCount;    // for calculation of lifetime};struct cudaAngle{    int type; // potential type    float p0, p1, p2;    // potential parameters    void (*force_eng)(int4* angle, cudaAngle* type, cudaMD* md, float& eng);};

__global__ void apply_bonds(int iStep, int bndPerBlock, int bndPerThread, cudaMD* md){    int def;    int id1, id2;       // atom indexes    int old, old_spec2, spec1, spec2, new_spec1, new_spec2;     // atom types    int new_bond_type;        int save_lt, need_r, loop;    // flags to save lifetime, to need to calculate r^2 and to be in while loop    int mnmx;   // flag minimum or maximum    int action; // flag: 0 - do nothing, 1 - delete bond, 2 - transform bond    cudaBond *old_bnd, *cur_bnd;// old bond type, current bond type    float dx, dy, dz, r2, r;    float f, eng = 0.0f;    __shared__ float shEng;#ifdef DEBUG_MODE    int cnt;    // count of change spec2 loops#endif    if (threadIdx.x == 0)    {        shEng = 0.0f;    }    __syncthreads();    int id0 = blockIdx.x * bndPerBlock + threadIdx.x * bndPerThread;    int N = min(id0 + bndPerThread, md->nBond);    int iBnd;    for (iBnd = id0; iBnd < N; iBnd++)      if (md->bonds[iBnd].z)  // the bond is not broken      {          // atom indexes          id1 = md->bonds[iBnd].x;          id2 = md->bonds[iBnd].y;          // atom types          spec1 = md->types[id1];          spec2 = md->types[id2];          old_bnd = &(md->bondTypes[md->bonds[iBnd].z]);          cur_bnd = old_bnd;          save_lt = 0;          need_r = 1;          loop = 1;#ifdef DEBUG_MODE          cnt = 0;#endif                    if ((cur_bnd->spec1 == spec1)&&(cur_bnd->spec2 == spec2))          {              //ok          }          else              if ((cur_bnd->spec1 == spec2) && (cur_bnd->spec2 == spec1))              {                  invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));                  //... then ok              }              else // atom types do not correspond to bond types              {                  save_lt = 1;              }          // end initial stage          while (loop)          {             if (save_lt)                    {                  def = md->def_bonds[spec1][spec2];                  if (def == 0)     // these atom types do not form a bond                  {#ifdef DEBUG_MODE                      printf("probably, something goes wrong\n");#endif                      action = 1;   // delete                      break;                  }                  else                  {                      //! change bond type and go on                      if (def < 0)                        {                          invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));                          def = -def;                      }                      md->bonds[iBnd].z = def;                      cur_bnd = &(md->bondTypes[def]);                  }             }  // end if (save_lt)             // calculate distance (only once)             if (need_r)             {                dx = md->xyz[id1].x - md->xyz[id2].x;                dy = md->xyz[id1].y - md->xyz[id2].y;                dz = md->xyz[id1].z - md->xyz[id2].z;                delta_periodic(dx, dy, dz, md);                r2 = dx * dx + dy * dy + dz * dz;                need_r = 0;             }             action = 0;   // 0 - just cultivate bond 1 - delete bond 2 - transform bond             if ((cur_bnd->mxEx) && (r2 > cur_bnd->r2max))             {                 mnmx = 1;                 if (cur_bnd->new_type[mnmx] == 0)  // delete bond                   action = 1;                else                   action = 2;   // modify bond             }             else if ((cur_bnd->mnEx) && (r2 < cur_bnd->r2min))             {                 mnmx = 0;                 action = 2;   // at minimum only bond modification possible             }             // end select action             // try to change atom types (if needed)             if (action)             {                 save_lt = 1;                 new_spec1 = cur_bnd->new_spec1[mnmx];                 new_spec2 = cur_bnd->new_spec2[mnmx];                 //the first atom                 old = atomicCAS(&(md->types[id1]), spec1, new_spec1);                 if (old != spec1)                 {                     spec1 = old;                     spec2 = md->types[id2];   // refresh type of the 2nd atom                     // return to begin of the while loop                 }                 else      // types[id1] have been changed                 {#ifdef USE_NEWANG   // save changes in atom type                     atomicCAS(&(md->oldTypes[id1]), -1, spec1);#endif                     old_spec2 = spec2;                     while ((old = atomicCAS(&(md->types[id2]), old_spec2, new_spec2)) != old_spec2)                     {                         //! the worst variant: this thread changes atom 1, other thread changes atom 2                         // imagine that we had A-old bond with the same behavior                         def = md->def_bonds[spec1][old];#ifdef DEBUG_MODE                         if (def == 0)                         {                             printf("UBEH[001]: in apply_bonds, change atom types. There are no bond types between Species[%d] and [%d]\n", spec1, old);                             break;                         }#endif                         if (def < 0)  // spec1 -> new_spec2 spec2 -> newSpec1                         {                             cur_bnd = &(md->bondTypes[-def]);                             new_spec2 = cur_bnd->new_spec1[mnmx];                         }                         else // direct order                         {                             cur_bnd = &(md->bondTypes[def]);                             new_spec2 = cur_bnd->new_spec2[mnmx];                         }                         old_spec2 = old;#ifdef DEBUG_MODE                         cnt++;                         if (cnt > 10)                         {                             printf("UBEH[002]: too many atempst to change spec2 = %d\n", spec2);                             break;                         }#endif                     }#ifdef USE_NEWANG   // save changes in atom type                     atomicCAS(&(md->oldTypes[id2]), -1, spec2);#endif                     loop = 0;                 }                 //end change types             } // end if (action)             else               loop = 0;    // action == 0, out of cycle          }  // end while(loop)          if (action == 2)          {              new_bond_type = cur_bnd->new_type[mnmx];              if (new_bond_type < 0)              {                  md->bonds[iBnd].x = id2;                  md->bonds[iBnd].y = id1;                  new_bond_type = -new_bond_type;              }              md->bonds[iBnd].z = new_bond_type;              cur_bnd = &(md->bondTypes[new_bond_type]);          }          // perform calculations and save mean bond length          if (action != 1)  // not delete          {              r = sqrt(r2);              f = cur_bnd->force_eng(r2, r, eng, cur_bnd);              atomicAdd(&(md->frs[id1].x), f * dx);              atomicAdd(&(md->frs[id2].x), -f * dx);              atomicAdd(&(md->frs[id1].y), f * dy);              atomicAdd(&(md->frs[id2].y), -f * dy);              atomicAdd(&(md->frs[id1].z), f * dz);              atomicAdd(&(md->frs[id2].z), -f * dz);                            atomicAdd(&(cur_bnd->rSumm), r);              atomicAdd(&(cur_bnd->rCount), 1);          }          else      //delete bond          {// decrease the number of bonds for atoms              atomicSub(&(md->nbonds[id1]), 1);              atomicSub(&(md->nbonds[id2]), 1);              md->bonds[iBnd].z = 0;              // change parents              exclude_parents(id1, id2, md);          }          if (save_lt)          {              keep_bndlifetime(iStep, &(md->bonds[iBnd]), old_bnd);              if (action != 1) // not delete                atomicAdd(&(cur_bnd->count), 1);              atomicSub(&(old_bnd->count), 1);          }      } // end main loop      // split energy to shared and then to global memory      atomicAdd(&shEng, eng);      __syncthreads();      if (threadIdx.x == 0)          atomicAdd(&(md->engBond), shEng);}

__device__ void exclude_parents(int id1, int id2, cudaMD* md)// exclude id1 and id2 from parents of each other (if they are)//  and seek other parents if able{    // flags to clear parent    int clear_1 = 0;        int clear_2 = 0;    int i, flag;        if (md->parents[id1] == id2)         clear_1 = 1;    if (md->parents[id2] == id1)        clear_2 = 1;    i = 0;    while ((i < md->nBond) && (clear_1 || clear_2))    {        if (md->bonds[i].z != 0)        {            flag = 0;            if (clear_1)            {                if (md->bonds[i].x == id1)                {                    md->parents[id1] = md->bonds[i].y;                    flag = 1;                }                else if (md->bonds[i].y == id1)                {                    md->parents[id1] = md->bonds[i].y;                    flag = 1;                }                if (flag)                {                    clear_1 = 0;                    i++;                    continue;                }            }            if (clear_2)            {                if (md->bonds[i].x == id2)                {                    md->parents[id2] = md->bonds[i].y;                    flag = 1;                }                else if (md->bonds[i].y == id2)                {                    md->parents[id2] = md->bonds[i].y;                    flag = 1;                }                if (flag)                {                    clear_2 = 0;                    i++;                    continue;                }            }        }        i++;    }    // be on the safe side    if (clear_1)            md->parents[id1] = -1;    if (clear_2)        md->parents[id2] = -1;}

__device__ void try_to_bind(float r2, int id1, int id2, int spec1, int spec2, cudaMD *md){    int r2Int;      //  (int)r2 * const    // save parents to exclude re-linking    if (md->parents[id1] == id2)        return;    if (md->parents[id2] == id1)        return;    if (md->bindBonds[spec1][spec2] != 0)    {        if (r2 < md->bindR2[spec1][spec2])        {            r2Int = (int)(r2 * 100);            if (atomicMin(&(md->r2Min[id1]), r2Int) > r2Int)    // replace was sucessfull            {                md->neighToBind[id1] = id2 + 1;     // as 0 is reserved for no neighbour                md->canBind[id1] = 1;            }            // similar for the second atom            if (atomicMin(&(md->r2Min[id2]), r2Int) > r2Int)    // replace was sucessfull            {                md->neighToBind[id2] = id1 + 1;     // as 0 is reserved for no bind                md->canBind[id2] = 1;            }        }    }}

__global__ void create_bonds(int iStep, int atPerBlock, int atPerThread, cudaMD* md)// connect atoms which are selected to form bonds{    int id1, id2, nei;    // neighbour index    int btype, bind;    // bond type index and bond index    cudaBond* bnd;    int spec1, spec2;   // species indexes        int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    int iat;    for (iat = id0; iat < N; iat++)    {        nei = md->neighToBind[iat];        if (nei)    // neighbour exists        {            nei--;  // (nei = spec_index + 1)            if (iat < nei)            {                id1 = iat;                id2 = nei;            }            else            {                id1 = nei;                id2 = iat;            }                        // try to lock the first atom            if (atomicCAS(&(md->canBind[id1]), 1, 0) == 0)  // the atom is already used                continue;            // try to lock the second atom            if (atomicCAS(&(md->canBind[id2]), 1, 0) == 0)  // the atom is already used            {                // unlock the first one back                atomicExch(&(md->canBind[id1]), 1);                continue;            }            // create bond iat-nei            bind = atomicAdd(&(md->nBond), 1);#ifdef DEBUG_MODE            if (bind >= md->mxBond)            {                printf("UBEH[003]: Exceed maximal number of bonds, %d\n", md->mxBond);            }#endif            spec1 = md->types[id1];            spec2 = md->types[id2];#ifdef USE_NEWANG   // save that we have changed atom type            atomicCAS(&(md->oldTypes[id1]), -1, spec1);            atomicCAS(&(md->oldTypes[id2]), -1, spec2);#endif            btype = md->bindBonds[spec1][spec2];                        if (btype < 0)            {                // invert atoms order                md->bonds[bind].x = id2;                md->bonds[bind].y = id1;                md->bonds[bind].z = -btype;                bnd = &(md->bondTypes[-btype]);                // change atom types according the formed bond                md->types[id1] = bnd->spec2;                md->types[id2] = bnd->spec1;            }            else             {                md->bonds[bind].x = id1;                md->bonds[bind].y = id2;                md->bonds[bind].z = btype;                bnd = &(md->bondTypes[btype]);                // change atom types according the formed bond                md->types[id1] = bnd->spec1;                md->types[id2] = bnd->spec2;            }                        atomicAdd((&bnd->count), 1);            md->bonds[bind].w = iStep;  // keep time of the bond creation for lifetime calculation            atomicAdd(&(md->nbonds[id1]), 1);            atomicAdd(&(md->nbonds[id2]), 1);            // replace parents if none:            atomicCAS(&(md->parents[id1]), -1, id2);            atomicCAS(&(md->parents[id2]), -1, id1);        }    }    // end loop by atoms}

__global__ void clear_bonds(cudaMD* md)// clear bonds with .z == 0{    int i = 0;    int j = md->nBond - 1;    while (i < j)    {        while ((md->bonds[j].z == 0) && (j > i))            j--;        while ((md->bonds[i].z != 0) && (i < j))            i++;        if (i < j)        {            md->bonds[i] = md->bonds[j];            md->bonds[j].z = 0;            i++;            j--;        }    }    if ((i == j) && (md->bonds[i].z == 0))        md->nBond = j;    else        md->nBond = j + 1;}

Валентные углы

__global__ void refresh_angles(int iStep, int atPerBlock, int atPerThread, cudaMD *md)// delete old angles and create new ones for atoms which change their type{int i, j, n, t, ang;int nei[8];// bonded neighbors of given atomint cnt;int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;int N = min(id0 + atPerThread, md->nAt);int iat;for (iat = id0; iat < N; iat++)if (md->oldTypes[iat] != -1){i = 0;n = md->nangles[iat];while (n && (i < md->nAngle)){if (md->angles[i].w)if (md->angles[i].x == iat){n--;md->angles[i].w = 0;}i++;}// create new anglest = md->specAngles[md->types[iat]];// get type of angle, which formed by current atom typeif (t && (md->nbonds[iat] > 1))// atom type supports angle creating and number of bonds is enough{// search of neighbors by bondsi = 0; cnt = 0;n = md->nbonds[iat];while (n && (i < md->nBond)){if (md->bonds[i].z)// if bond isn't deleted{if (md->bonds[i].x == iat){nei[cnt] = md->bonds[i].y;cnt++;n--;}else if (md->bonds[i].y == iat){nei[cnt] = md->bonds[i].x;cnt++;n--;}}i++;}// add new angles based on found neighbors:for (i = 0; i < cnt-1; i++)for (j = i + 1; j < cnt; j++){ang = atomicAdd(&(md->nAngle), 1);md->angles[ang].x = iat;md->angles[ang].y = nei[i];md->angles[ang].z = nei[j];md->angles[ang].w = t;}n = (cnt * (cnt - 1)) / 2;}md->nangles[iat] = n;// reset flagmd->oldTypes[iat] = -1;}}

__global__ void apply_angles(int iStep, int angPerBlock, int angPerThread, cudaMD* md)// apply valence angle potentials{cudaAngle* ang;// energies of angle potentialfloat eng;__shared__ float shEng;if (threadIdx.x == 0)shEng = 0.0f;__syncthreads();int id0 = blockIdx.x * angPerBlock + threadIdx.x * angPerThread;int N = min(id0 + angPerThread, md->nAngle);int i;for (i = id0; i < N; i++)if (md->angles[i].w){ang = &(md->angleTypes[md->angles[i].w]);ang->force_eng(&(md->angles[i]), ang, md, eng);}// split energy to shared and then to global memoryatomicAdd(&shEng, eng);__syncthreads();if (threadIdx.x == 0)atomicAdd(&(md->engAngl), shEng);}

__device__ void angle_hcos(int4* angle, cudaAngle* type, cudaMD* md, float& eng)// harmonic cosine valent angle potential:// U = k / 2 * (cos(th)-cos(th0))^{float k = type->p0;float cos0 = type->p1;// indexes of central atom and ligands:int c = angle->x;int l1 = angle->y;int l2 = angle->z;// vector ijfloat xij = md->xyz[l1].x - md->xyz[c].x;float yij = md->xyz[l1].y - md->xyz[c].y;float zij = md->xyz[l1].z - md->xyz[c].z;delta_periodic(xij, yij, zij, md);float r2ij = xij * xij + yij * yij + zij * zij;float rij = sqrt(r2ij);// vector ikfloat xik = md->xyz[l2].x - md->xyz[c].x;float yik = md->xyz[l2].y - md->xyz[c].y;float zik = md->xyz[l2].z - md->xyz[c].z;delta_periodic(xik, yik, zik, md);float r2ik = xik * xik + yik * yik + zik * zik;float rik = sqrt(r2ik);float cos_th = (xij * xik + yij * yik + zij * zik) / rij / rik;float dCos = cos_th - cos0; // delta cosinusfloat c1 = -k * dCos;float c2 = 1.0 / rij / rik;atomicAdd(&(md->frs[c].x), -c1 * (xik * c2 + xij * c2 - cos_th * (xij / r2ij + xik / r2ik)));atomicAdd(&(md->frs[c].y), -c1 * (yik * c2 + yij * c2 - cos_th * (yij / r2ij + yik / r2ik)));atomicAdd(&(md->frs[c].z), -c1 * (zik * c2 + zij * c2 - cos_th * (zij / r2ij + zik / r2ik)));atomicAdd(&(md->frs[l1].x), c1 * (xik * c2 - cos_th * xij / r2ij));atomicAdd(&(md->frs[l1].y), c1 * (yik * c2 - cos_th * yij / r2ij));atomicAdd(&(md->frs[l1].z), c1 * (zik * c2 - cos_th * zij / r2ij));atomicAdd(&(md->frs[l2].x), c1 * (xij * c2 - cos_th * xik / r2ik));atomicAdd(&(md->frs[l2].y), c1 * (yij * c2 - cos_th * yik / r2ik));atomicAdd(&(md->frs[l2].z), c1 * (zij * c2 - cos_th * zik / r2ik));eng += 0.5 * k * dCos * dCos;}

Примеры

Заключительные комментарии

• Космики;
• Жизнь на частицах.

Использование CUDA Runtime API для вычислений. Сравнение CPU и GPU вычислений.

intdegree_of=2;intdegree_of_max=Number_degree_of_max[0];//Массивхранящийзначениемаксимальнойстепеничислаintx=thread;//номервыполняемогопотокаintmax_number=INPUT[0];//Массивхранящийчисло,котороенеобходимополучитьintNumber=1;intDegree;boolBREAK=false;//Переменнаядлязавершенияwhilewhile(degree_of<=degree_of_max&&!BREAK){Number=1;for(inti=0;i<degree_of;i++){Number*=x;Degree=degree_of;}if(Number==max_number){OUT_NUMBER[thread]=X;//OUT_NUMBERМассивхранящийчислакоторыенеобходимовозвестивстепеньDegreeдляполученияисходногочислаOUT_DEGREE[thread]=Degree;//OUT_DEGREEМассивхранящийстепеньвкоторуюнужновозвестичислоXдляполученияисходногочисла}degree_of++;//Вслучаевыходазапредел:if(degree_of>degree_of_max||Number>max_number){BREAK=true;}}

#include <iostream>#include<vector>#include<string>//необходимо для getline#include<thread>#include<fstream>using namespace std;int Running_thread_counter = 0;void Upload_to_CPU(unsigned long long  *Number, unsigned long long  *Stepn, bool *Stop,unsigned long long  *INPUT, unsigned long long  *max, int THREAD);void Upload_to_CPU(unsigned long long  *Number, unsigned long long  *Stepn, bool *Stop,unsigned long long  *INPUT, unsigned long long  *max, int THREAD) {int thread = THREAD;Running_thread_counter++;unsigned long long  MAX_DEGREE_OF = max[0];int X = thread;unsigned long long  Calculated_number = 1;unsigned long long  DEGREE_OF = 2;unsigned long long   INP = INPUT[0];Stop[thread] = false;bool BREAK = false;if (X != 0 && X != 1) {while (!BREAK) {if (DEGREE_OF <= MAX_DEGREE_OF) {Calculated_number = 1;for (int counter = 0; counter < DEGREE_OF; counter++) {Calculated_number *= X;}if (Calculated_number == INP) {Stepn[thread] = DEGREE_OF;Number[thread] = X;Stop[thread] = true;BREAK = true;}DEGREE_OF++;}else { BREAK = true; }}}}void Parallelize_to_threads(unsigned long long  *Number, unsigned long long  *Stepn, bool *Stop,unsigned long long  *INPUT, unsigned long long  *max, int size);int main(){int size = 1000;unsigned long long  *Number = new unsigned long long[size], *Degree_of = new unsigned long long[size];unsigned long long  *Max_Degree_of = new unsigned long long[1];unsigned long long  *INPUT_NUMBER = new unsigned long long[1];Max_Degree_of[0] = 7900;INPUT_NUMBER[0] = 216 * 216 * 216;ifstream inp("input.txt");if (inp.is_open()) {string t;vector<unsigned long long>IN;while (getline(inp, t)) {IN.push_back(stol(t));}INPUT_NUMBER[0] = IN[0];//исходное числоMax_Degree_of[0] = IN[1];//значение максимальной степени}else {ofstream error("error.txt");if (error.is_open()) {error << "No file " << '"' << "input.txt" << '"' << endl;error << "Please , create a file" << '"' << "input.txt" << '"' << endl;error << "One read:input number" << endl;error << "Two read:input max stepen" << endl;error << "." << endl;error.close();INPUT_NUMBER[0] = 1;Max_Degree_of[0] = 1;}}//расскометрируйте следующий код , если хотите видеть исходные значения в окне консоли //cout << INPUT[0] << endl;bool *Elements_that_need_to_stop = new bool[size];Parallelize_to_threads(Number, Degree_of, Elements_that_need_to_stop, INPUT_NUMBER, Max_Degree_of, size);vector<unsigned long long>NUMBER, DEGREEOF;for (int i = 0; i < size; i++) {if (Elements_that_need_to_stop[i]) {if (Degree_of[i] < INPUT_NUMBER[0] && Number[i] < INPUT_NUMBER[0]) {//проверка на ошибки NUMBER.push_back(Number[i]);DEGREEOF.push_back(Degree_of[i]);}}}//расскометрируйте следующий код , если хотите вывести результаты в консоль//это может замедлить программу /*for (int f = 0; f < NUMBER.size(); f++) {cout << NUMBER[f] << "^" << DEGREEOF[f] << "=" << INPUT_NUMBER[0] << endl;}*/ofstream out("out.txt");if (out.is_open()) {for (int f = 0; f < NUMBER.size(); f++) {out << NUMBER[f] << "^" << DEGREEOF[f] << "=" << INPUT_NUMBER[0] << endl;}out.close();}}void Parallelize_to_threads(unsigned long long  *Number, unsigned long long  *Stepn, bool *Stop,unsigned long long  *INPUT, unsigned long long  *max, int size) {thread *T = new thread[size];Running_thread_counter = 0;for (int i = 0; i < size; i++) {T[i] = thread(Upload_to_CPU, Number, Stepn, Stop, INPUT, max, i);T[i].detach();}while (Running_thread_counter < size - 1);//дождаться завершения выполнения всех потоков }

//библиотеки cuda_runtime.h и device_launch_parameters.h//для работы с cyda#include "cuda_runtime.h"#include "device_launch_parameters.h"#include<vector>#include<string>//для getline#include <stdio.h>#include<fstream>using namespace std;__global__ void Upload_to_GPU(unsigned long long  *Number,unsigned long long  *Stepn, bool *Stop,unsigned long long  *INPUT,unsigned long long  *max) {int thread = threadIdx.x;unsigned long long  MAX_DEGREE_OF = max[0];    int X = thread;unsigned long long  Calculated_number = 1;unsigned long long  Current_degree_of_number = 2;    unsigned long long   Original_numberP = INPUT[0];Stop[thread] = false;bool BREAK = false;if (X!=0&&X!=1) {while (!BREAK) {if (Current_degree_of_number <= MAX_DEGREE_OF) {Calculated_number = 1;for (int counter = 0; counter < Current_degree_of_number; counter++) { Calculated_number*=X;}if (Calculated_number == Original_numberP) {Stepn[thread] = Current_degree_of_number;Number[thread] = X;Stop[thread] = true;BREAK = true;}Current_degree_of_number++;}else { BREAK = true; }}}}cudaError_t Configure_cuda(unsigned long long *Number, unsigned long long  *Stepn, bool *Stop,unsigned long long  *INPUT, unsigned long long  *max,unsigned int size);int main(){int size = 1000;    unsigned long long  *Number=new unsigned long long [size], *Degree_of=new unsigned long long [size];unsigned long long  *Max_degree_of = new unsigned long long [1];unsigned long long  *INPUT_NUMBER = new unsigned long long [1];   Max_degree_of[0] = 7900;ifstream inp("input.txt");if (inp.is_open()) {string text;vector<unsigned long long>IN;while (getline(inp, text)) {IN.push_back( stol(text));}INPUT_NUMBER[0] = IN[0];Max_degree_of[0] = IN[1];}else {ofstream error("error.txt");if (error.is_open()) {error<<"No file "<<'"'<<"input.txt"<<'"'<<endl;error<<"Please , create a file" << '"' << "input.txt" << '"' << endl;error << "One read:input number" << endl;error << "Two read:input max stepen" << endl;error << "." << endl;error.close();INPUT_NUMBER[0] = 1;Max_degree_of[0] = 1;}}bool *Elements_that_need_to_stop = new bool[size];    // Загрузка массивов в cudacudaError_t cudaStatus =  Configure_cuda(Number, Degree_of, Elements_that_need_to_stop, INPUT_NUMBER, Max_degree_of, size);    if (cudaStatus != cudaSuccess) {        fprintf(stderr, "addWithCuda failed!");        return 1;    }vector<unsigned long long>NUMBER, DEGREEOF;for (int i = 0; i < size; i++) {if (Elements_that_need_to_stop[i]) {NUMBER.push_back(Number[i]);//занести в вектор числоDEGREEOF.push_back(Degree_of[i]);//занести в вектор степень числа}}//раскоментируйте следующий код , чтобы вывести результаты в консоль/*for (int f = 0; f < NUMBER.size(); f++) {cout << NUMBER[f] << "^" << DEGREEOF[f] << "=" << INPUT_NUMBER[0] << endl;}*/ofstream out("out.txt");if (out.is_open()) {for (int f = 0; f < NUMBER.size(); f++) {out << NUMBER[f] << "^" << DEGREEOF[f] << "=" << INPUT_NUMBER[0] << endl;}out.close();}    //Очистить ресурсы связанные с устройством    cudaStatus = cudaDeviceReset();    if (cudaStatus != cudaSuccess) {        fprintf(stderr, "cudaDeviceReset failed!");        return 1;    }    return 0;}cudaError_t  Configure_cuda(unsigned long long  *Number, unsigned long long *Degree_of, bool *Stop,unsigned long long *INPUT, unsigned long long *max,unsigned int size) {unsigned long long *dev_Number = 0;unsigned long long *dev_Degree_of = 0;unsigned long long *dev_INPUT = 0;unsigned long long *dev_Max = 0;bool *dev_Elements_that_need_to_stop;cudaError_t cudaStatus;// УСТАНОВКА ИСПОЛЬЗУЕМОГО GPU cudaStatus = cudaSetDevice(0);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");goto Error;}// РЕЗЕРВИРОВАНИЕ МЕСТА В ПАМЯТИ ПОД ДАННЕcudaStatus = cudaMalloc((void**)&dev_Number, size * sizeof(unsigned long long));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!dev_Number");goto Error;}cudaStatus = cudaMalloc((void**)&dev_Degree_of, size * sizeof(unsigned long long));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!dev_Degree_of");goto Error;}cudaStatus = cudaMalloc((void**)&dev_Max, size * sizeof(unsigned long long int));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!dev_Max");goto Error;}cudaStatus = cudaMalloc((void**)&dev_INPUT, size * sizeof(unsigned long long));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!dev_INPUT");goto Error;}cudaStatus = cudaMalloc((void**)&dev_Elements_that_need_to_stop, size * sizeof(bool));if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMalloc failed!dev_Stop");goto Error;}// ПЕРЕМЕЩЕНИЕ ДАННХ В ПАМЯТЬ GPUcudaStatus = cudaMemcpy(dev_Max, max, size * sizeof(unsigned long long), cudaMemcpyHostToDevice);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!");goto Error;}cudaStatus = cudaMemcpy(dev_INPUT, INPUT, size * sizeof(unsigned long long), cudaMemcpyHostToDevice);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!");goto Error;}Upload_to_GPU<<<1, size>>>(dev_Number, dev_Degree_of, dev_Elements_that_need_to_stop, dev_INPUT, dev_Max);// Проверка сбоев ядраcudaStatus = cudaGetLastError();if (cudaStatus != cudaSuccess) {fprintf(stderr, "addKernel launch failed: %s\n", cudaGetErrorString(cudaStatus));goto Error;}// Ожидание завершения операций , выполняемых ядромcudaStatus = cudaDeviceSynchronize();if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching addKernel!\n", cudaStatus);goto Error;}// Перемещение данных из памяти GPU в системную памятьcudaStatus = cudaMemcpy(Number, dev_Number, size * sizeof(unsigned long long), cudaMemcpyDeviceToHost);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!");goto Error;}cudaStatus = cudaMemcpy(Degree_of, dev_Degree_of, size * sizeof(unsigned long long), cudaMemcpyDeviceToHost);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!");goto Error;}cudaStatus = cudaMemcpy(Stop, dev_Elements_that_need_to_stop, size * sizeof(bool), cudaMemcpyDeviceToHost);if (cudaStatus != cudaSuccess) {fprintf(stderr, "cudaMemcpy failed!");goto Error;}Error://Освобождение памяти GPU от данныхcudaFree(dev_INPUT);cudaFree(dev_Degree_of);cudaFree(dev_Max);cudaFree(dev_Elements_that_need_to_stop);cudaFree(dev_Number);return cudaStatus;}

__global__  

intRunning_thread_counter=0;

usingSystem;usingSystem.Collections.Generic;usingSystem.Linq;usingSystem.Text;usingSystem.Threading.Tasks;usingSystem.Diagnostics;usingSystem.IO;namespaceConsoleAppTESTSTEPEN_CPU_AND_GPU_{classProgram{staticstringUpload(Int64number,Int64degree_of){stringOUT="";string[]Chord_values=newstring[2];Int64Degree_of=degree_of;Int64Number=number;Chord_values[0]=Number.ToString();Chord_values[1]=Degree_of.ToString();File.WriteAllLines("input.txt",Chord_values);//файлвходныхданныхOUT+="inputnumber:"+Number.ToString()+"\n";OUT+="inputdegreeofnumber:"+Degree_of.ToString()+"\n";DateTimerunning_CPU_application=DateTime.Now;//записатьвремязапускапрограммыProcessproc=Process.Start("ConsoleApplication29.exe");//exeреализацияалгоритманаc++x64использующаяCPUдлявычисленийwhile(!proc.HasExited);//дождатсязавершенияпрограммыDateTimestop_CPU_application=DateTime.Now;//записатьвремяостановкипрограммыstring[]outs=File.ReadAllLines("out.txt");//получитьрезультатыFile.Delete("out.txt");OUT+="CPU:"+"\n";if(outs.Length>0){for(intj=0;j<outs.Length;j++){OUT+=outs[j]+"\n";}}else{OUT+="novalues"+"\n";}OUT+="running_CPU_application:"+running_CPU_application.ToString()+"\n";OUT+="stop_CPU_application:"+stop_CPU_application.ToString()+"\n";OUT+="GPU:"+"\n";//альтернативныедействиядляреализацииалгоритмаkorenXN.exex64использующегодлявычисленийGPUDateTimerunning_GPU_application=DateTime.Now;ProcessprocGPU=Process.Start("korenXN.exe");while(!procGPU.HasExited);DateTimestop_GPU_application=DateTime.Now;string[]outs2=File.ReadAllLines("out.txt");File.Delete("out.txt");if(outs2.Length>0){for(intj=0;j<outs2.Length;j++){OUT+=outs2[j]+"\n";}}else{OUT+="novalues"+"\n";}OUT+="running_GPU_application:"+running_GPU_application.ToString()+"\n";OUT+="stop_GPU_application:"+stop_GPU_application.ToString()+"\n";returnOUT;//возвратитьрезультат}staticvoidMain(){Int64start=36*36;//начальноезначениевходногочислаInt64degree_of_strat=500;//начальноезначениемаксимальнойстепениintsize=20-5;//количествоэлементоввмассивеInt64[]Number=newInt64[size];//массиввходныхчиселInt64[]Degree_of=newInt64[size];//массивмаксимальныхстепенейstring[]outs=newstring[size];//масссиврезультатовfor(intn=0;n<size;n++){if(n%2==0){Number[n]=start*start;}else{Number[n]=start*degree_of_strat;Number[n]-=n+n;}start+=36*36;Degree_of[n]=degree_of_strat;degree_of_strat+=1000;}for(intn=0;n<size;n++){outs[n]=Upload(Number[n],Degree_of[n]);Console.WriteLine(outs[n]);}System.IO.File.WriteAllLines("result.txt",outs);//записатьрезультатывфайлresult.txt}}}

график

График