Метки

, , , , , , , , ,

CUDA является расширением языка C++ для параллельных вычислений с помощью видеокард NVIDIA. Перенос вычислений на видеокарту может оказаться оптимальным для решения задач с применением одинаковых действий к большому числу элементов. Прирост скорости на порядки был получен в обработке томографических изображений, так же CUDA нашла применение в моделировании физических систем и программах восстановления паролей (см. тесты http://www.ixbt.com/video3/cuda-1.shtml и http://www.ixbt.com/video3/cuda-2.shtml).

Бешенный прирост скорости по сравнению с центральными процессорами персоналок и даже с большинством общедоступных для академического использования суперкомпьютерах и кластерах достигается благодаря нескольким фактам. Во-первых, даже на недорогой видеокарте вычислениями заняты десятки процессоров, так в уже довольно старых чипсетах GeForce GT220 48 ядер, а GeForce GT 280 аж 240. Так же существуют специальные видеокарты специально для параллельных вычислений (линейка чипсетов Tesla с производительностью порядка 500 ГФЛОПС на арифметике двойной точности).

Существуют и ограничения для использования. Так Tesla со стоимостью в несколько тысяч долларов для домашнего использования могут оказаться дороговаты, а видеокарты на чипсетах геймерского класса могут сдать на порядки при смене типов переменных с float на double. Оптимальным является как можно меньшее число ветвлений в коде, непосредственно исполняемом на видеокарте: видеопроцессоры группируются в так называемые мультипроцессоры, как правило объединяющие 8 ядер, которые могут выполнять только одинаковый код над разными участками памяти. В случае ветвления кода мультипроцессором сначала будет выполняться одна ветка, затем иная, что приводит к простаиванию ядер и увеличению времени выполнения кода. И как часто бывает, наиболее значимым является использование памяти, которой насчитывается несколько уровней (чем быстрее тип памяти, тем меньше его будет доступно на ядро).

Ниже представлен код интегрирования системы дифференциальных уравнений
\dot{\mathbf y} = f(t, \mathbf y),
где \mathbf y — вектор из N переменных методами Рунге-Кутты 4-го порядка и Рунге-Кутты-Фельдберга 4-5-го порядков (по описанию на страничке Фуллертоновского университета — тут надо отметить, что там же у них лежит пдф с описанием алгоритма и опечаткой в одном из коеффициентов). Лицензия As Is. Естественно, код не претендует на особую оригинальность и может быть значительно улучшен 🙂

Идея состоит в распараллеливании вычисления производной и операций над векторами. В качестве тестового используется никакой пример интегрирования уравнения движения обычного осциллятора.

Общие определения

Файл funcdef.cuh. Вначале определим тип TYPE для быстрого переключения между типами float и double при необходимости. Тип функции вычисления производных задается в обобщенном виде через typedef Deriv. Произвольные параметры в неё можно передать с помощью редактирования под свои нужды структуры DerivParams.

///////////////////////////////////////////////////////////////////////
//
// ODE system solving by CUDA
// Main types definition
// Source: https://engraver.wordpress.com/
//
#ifndef _FUNC_DEF_
#define _FUNC_DEF_

#include

////////////////////////////////////////////////////////////////////////
//
// Type definition for fast switching between float and double precision
//
typedef float TYPE;

////////////////////////////////////////////////////////////////////////
//
// Structure for sending all needed variables to function for derivative
// calculation
//
struct DerivParams
{
	TYPE omega;
};

////////////////////////////////////////////////////////////////////////
//
// Definition for function for derivative calculation in generic form
// (all parameters can be passed by DerivParams type editing)
//
typedef void (*Deriv)(int, TYPE, TYPE *, TYPE *, TYPE, DerivParams * params);

////////////////////////////////////////////////////////////////////////
//
// CUDA integration management (needed in runge_cuda.cu and funct.cu)
//
const unsigned int RUNGE_B_SIZE = 5; // 5
const unsigned int RUNGE_BLOCK_SIZE = 1 << RUNGE_B_SIZE;

#endif

Хеадер с прототипом функции интегрирования

Файл runge_cuda.cuh. В rungeCuda задается количество уравнений N в системе, указатель на функцию, возвращающую производную func, вектор начальных значений y0, вектор для хранения результата y, интервал аргументов [t0 t1], начальный шаг интегрирования h0 (как указатель, возвращается последнее значение), точность tol. Когда параметр verbose принимает значение true, в терминал печатается информация об изменении шага интегрирования и прочее. В переменной time возвращается время вычислений, а параметр method определяет алгоритм, которым будет выполняться интегрирование: Рунге-Кутты или Рунге-Кутты-Фельберга.

Все вызовы CUDA происходят через cudaSafeCall для диагностики ошибок, если такова происходит.

#ifndef _RUNGE_CUDA_H_
#define _RUNGE_CUDA_H_

/*
 * ********************************************************************
 * \dot y = f(t), t \in (t0, t1), y(t0) = y0
 * ********************************************************************
 */

#include "funcdef.cuh"

////////////////////////////////////////////////////////////////////////
//
// Constants for selecting integration method
//
enum RungeMethod
{
	RM_RK4,   // Runge-Kutta 4th order (RK4)
	RM_RKF45  // Runge-Kutta-Fehlbert 4-5th order (RKF45)
};

////////////////////////////////////////////////////////////////////////
//
// Constants for exit codes
//
enum RungeExitCodes
{
	RUNGE_SUCCESS, // All Ok
	RUNGE_FAIL     // Error
};

////////////////////////////////////////////////////////////////////////
//
// Integration process management for RKF45
//
const TYPE RUNGE_MIN_H = 1e-5;	// Minimal step
const TYPE RUNGE_MAX_H = 0.75;	// Maximal step
const TYPE RUNGE_MIN_S = 0.1;	// Minimal step multiplier
const TYPE RUNGE_MAX_S = 5.0;	// Maximal step multiplier

////////////////////////////////////////////////////////////////////////
//
// Integration function
//
RungeExitCodes rungeCuda(int N, // Problem size
								Deriv func,     // Function, that return derivative
								TYPE * y0,      // Initial values
								TYPE * y,       // Result
								TYPE t0,        // Start argument
								TYPE t1,        // End argument
								TYPE *h0,       // Step (contain last step on exit)
								TYPE tol,       // Tolerance
								bool	verbose,  // Make a noise if true 🙂
								DerivParams * params,// Parameters for func
								TYPE * time,    // Calculation time
								RungeMethod method   // Integration method: RK4 or RKF45
								);

////////////////////////////////////////////////////////////////////////
//
// For save call of CUDA functions: print error messages to terminal,
// if exists.
//
RungeExitCodes cudaSafeCall(cudaError_t err, TYPE num);

#endif

Функция интегрирования

Файл runge_cuda.cu. Распараллеливание непосредственно функции интегрирования (а не вычисления производной, которое параллелим в ином месте) касается операций с векторами: сложение массивов, хранящих коэффициенты Рунге-Кутты и поиск максимального элемента при вычислении ошибки метода на итерации.

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include "runge_cuda.cuh"

////////////////////////////////////////////////////////////////////////
//
// All constants for Runge-Kutta and Runge-Kutta-Fehlber algorithm of
// 4th and 5th order
//
const TYPE a1 = 1.0/4; // Integration constants here...
const TYPE a2 = 3.0/8; // Tabbed constants for RKF45, other for RK4 too.
const TYPE a3 = 12.0/13;
 const TYPE a4 = 1.0/2;

const TYPE b1 = 1.0/4;
const TYPE b2 = 3.0/32;
const TYPE b3 = 9.0/32;
const TYPE b4 = 1932.0/2197;
const TYPE b5 = -7200.0/2197;
const TYPE b6 = 7296.0/2197;
const TYPE b7 = 439.0/216;
const TYPE b8 = -8.0;
const TYPE b9 = 3680.0/513;
const TYPE b10 = -845.0/4104;
 const TYPE b11 = -8.0/27;
 const TYPE b12 = 2.0;
 const TYPE b13 = -3544.0/2565;
 const TYPE b14 = 1859.0/4104;
 const TYPE b15 = -11.0/40;

const TYPE c1 = 25.0/216.0;
const TYPE c2 = 1408.0/2565.0;
const TYPE c3 = 2197.0/4104.0;
const TYPE c4 = -1.0/5.0;
 const TYPE c5 = 16.0/135;
 const TYPE c6 = 6656.0/12825;
 const TYPE c7 = 28561.0/56430;
 const TYPE c8 = -9.0/50;
 const TYPE c9 = 2.0/55;

////////////////////////////////////////////////////////////////////////
//
// Global variables
//
int runge_sizeByte = 0;
static TYPE * runged_k1;
static TYPE * runged_k2;
static TYPE * runged_k3;
static TYPE * runged_k4;
static TYPE * runged_k5;
 static TYPE * runged_k6;
static TYPE * runged_y;
static TYPE * runged_y1;
static TYPE * runged_z;
static TYPE * runged_ytemp;

////////////////////////////////////////////////////////////////////////
//
// Function prototypes
//
RungeExitCodes initArrays(RungeMethod method);
RungeExitCodes freeArrays(RungeMethod method);
bool checkRKFError(int N, TYPE * h, TYPE * s, TYPE tol, RungeMethod method);
void vectorAdd(int N, TYPE * d_v1, TYPE * d_v2, TYPE multiplier);
void calcDiff(int N, TYPE * d_v1, TYPE * d_v2);
__global__ void addKernel(int N, TYPE * d_a, TYPE * d_b, TYPE multiplier);
__global__ void calcDiffKernel(int N, TYPE * d_a, TYPE * d_b);
template<unsigned int blockSize> __global__ void reduceKernel(int N, TYPE * d_xin, TYPE * d_xout);

////////////////////////////////////////////////////////////////////////
//
// Main function for integration
//
RungeExitCodes rungeCuda(int N,                // Problem size
 Deriv func,     // Function, that return derivative
 TYPE * y0,        // Initial values
 TYPE * y,        // Result
 TYPE t0,            // Start argument
 TYPE t1,            // End argument
 TYPE *h0,        // Step (contain last step on exit)
 TYPE tol,        // Tolerance
 bool    verbose,// Make a noise if true 🙂
 DerivParams * params,    // Parameters for func
 TYPE * time,    // Calculation time
 RungeMethod method // Integration method: RK4 or RKF45
 )
{
 runge_sizeByte = N * sizeof(TYPE);
 initArrays(method);
 if (method == RM_RKF45) if (cudaSafeCall(cudaMemcpy(runged_y, y0, runge_sizeByte, cudaMemcpyHostToDevice), 1) != RUNGE_SUCCESS) return RUNGE_FAIL;

 TYPE t = t0;
 TYPE h = 2*(*h0);
 TYPE h1; // For RK4
 TYPE s;

 DerivParams * sparams;
 if (cudaSafeCall(cudaMalloc((void**)&sparams, sizeof(DerivParams)), 1.3) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMemcpy(sparams, params, sizeof(DerivParams), cudaMemcpyHostToDevice), 1.6) != RUNGE_SUCCESS) return RUNGE_FAIL;

 cudaEvent_t start, stop;
 cudaEventCreate(&start);
 cudaEventCreate(&stop);
 *time = 0.0;
 cudaEventRecord(start, 0);

 ////////////////////////////////////////////////////////////////////////
 //
 // Integration by Runge-Kutta method of 4th order
 //
 if (method == RM_RK4)
 {
 do
 {
 h = h/2;
 h1 = h/2;
 if (verbose)
 {
 printf("RKF4h2>> step: %f %f\n", h, h1);
 }

 // integrate with step h
 t = t0;
 if (cudaSafeCall(cudaMemcpy(runged_y, y0, runge_sizeByte, cudaMemcpyHostToDevice), 1) != RUNGE_SUCCESS) return RUNGE_FAIL;
 while (t < t1 + 0.5*h)
 {
 // 1
 func(N, t, runged_y, runged_k1, h, sparams);
 // 2
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 2) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b1);
 func(N, t + a1*h, runged_y1, runged_k2, h, sparams);
 // 3
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 3) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b2);
 vectorAdd(N, runged_y1, runged_k2, b3);
 func(N, t + a2*h, runged_y1, runged_k3, h, sparams);
 // 4
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 4) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b4);
 vectorAdd(N, runged_y1, runged_k2, b5);
 vectorAdd(N, runged_y1, runged_k3, b6);
 func(N, t + a3*h, runged_y1, runged_k4, h, sparams);
 // 5
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 5) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b7);
 vectorAdd(N, runged_y1, runged_k2, b8);
 vectorAdd(N, runged_y1, runged_k3, b9);
 vectorAdd(N, runged_y1, runged_k4, b10);
 func(N, t + h, runged_y1, runged_k5, h, sparams);

 // y_temp
 if (cudaSafeCall(cudaMemcpy(runged_ytemp, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 7) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_ytemp, runged_k1, c1);
 vectorAdd(N, runged_ytemp, runged_k3, c2);
 vectorAdd(N, runged_ytemp, runged_k4, c3);
 vectorAdd(N, runged_ytemp, runged_k5, c4);
 if (cudaSafeCall(cudaMemcpy(runged_y, runged_ytemp, runge_sizeByte, cudaMemcpyDeviceToDevice), 1) != RUNGE_SUCCESS) return RUNGE_FAIL;

 t += h;
 }

 // integrate with step h/2
 if (cudaSafeCall(cudaMemcpy(runged_y, y0, runge_sizeByte, cudaMemcpyHostToDevice), 1) != RUNGE_SUCCESS) return RUNGE_FAIL;
 t = t0;
 while (t < t1 + 0.5*h1)
 {
 // 1
 func(N, t, runged_y, runged_k1, h1, sparams);
 // 2
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 2) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b1);
 func(N, t + a1*h1, runged_y1, runged_k2, h1, sparams);
 // 3
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 3) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b2);
 vectorAdd(N, runged_y1, runged_k2, b3);
 func(N, t + a2*h1, runged_y1, runged_k3, h1, sparams);
 // 4
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 4) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b4);
 vectorAdd(N, runged_y1, runged_k2, b5);
 vectorAdd(N, runged_y1, runged_k3, b6);
 func(N, t + a3*h1, runged_y1, runged_k4, h1, sparams);
 // 5
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 5) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b7);
 vectorAdd(N, runged_y1, runged_k2, b8);
 vectorAdd(N, runged_y1, runged_k3, b9);
 vectorAdd(N, runged_y1, runged_k4, b10);
 func(N, t + h1, runged_y1, runged_k5, h1, sparams);

 // z
 if (cudaSafeCall(cudaMemcpy(runged_z, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 7) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_z, runged_k1, c1);
 vectorAdd(N, runged_z, runged_k3, c2);
 vectorAdd(N, runged_z, runged_k4, c3);
 vectorAdd(N, runged_z, runged_k5, c4);
 if (cudaSafeCall(cudaMemcpy(runged_y, runged_z, runge_sizeByte, cudaMemcpyDeviceToDevice), 1) != RUNGE_SUCCESS) return RUNGE_FAIL;

 t += h1;
 }
 }
 while (checkRKFError(N, &h, &s, tol, method));
 if (cudaSafeCall(cudaMemcpy(y, runged_z, runge_sizeByte, cudaMemcpyDeviceToHost), 9) != RUNGE_SUCCESS) return RUNGE_FAIL;
 *h0 = h1;
 }

 ////////////////////////////////////////////////////////////////////////
 //
 // Integration by Runge-Kutta-Fehlber method of 4th-5th order
 //
 if (method == RM_RKF45)
 {
 while (t < t1 + 0.5*h)
 {
 s = 1;
 do
 {
 h *= s;
 if (h < RUNGE_MIN_H) h = RUNGE_MIN_H;
 if (h > RUNGE_MAX_H) h = RUNGE_MAX_H;
 // 1
 func(N, t, runged_y, runged_k1, h, sparams);
 // 2
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 2) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b1);
 func(N, t + a1*h, runged_y1, runged_k2, h, sparams);
 // 3
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 3) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b2);
 vectorAdd(N, runged_y1, runged_k2, b3);
 func(N, t + a2*h, runged_y1, runged_k3, h, sparams);
 // 4
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 4) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b4);
 vectorAdd(N, runged_y1, runged_k2, b5);
 vectorAdd(N, runged_y1, runged_k3, b6);
 func(N, t + a3*h, runged_y1, runged_k4, h, sparams);
 // 5
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 5) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b7);
 vectorAdd(N, runged_y1, runged_k2, b8);
 vectorAdd(N, runged_y1, runged_k3, b9);
 vectorAdd(N, runged_y1, runged_k4, b10);
 func(N, t + h, runged_y1, runged_k5, h, sparams);
 // 6
 if (cudaSafeCall(cudaMemcpy(runged_y1, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 6) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_y1, runged_k1, b11);
 vectorAdd(N, runged_y1, runged_k2, b12);
 vectorAdd(N, runged_y1, runged_k3, b13);
 vectorAdd(N, runged_y1, runged_k4, b14);
 vectorAdd(N, runged_y1, runged_k5, b15);
 func(N, t + a4*h, runged_y1, runged_k6, h, sparams);

 // y_temp
 if (cudaSafeCall(cudaMemcpy(runged_ytemp, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 7) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_ytemp, runged_k1, c1);
 vectorAdd(N, runged_ytemp, runged_k3, c2);
 vectorAdd(N, runged_ytemp, runged_k4, c3);
 vectorAdd(N, runged_ytemp, runged_k5, c4);

 // z
 if (cudaSafeCall(cudaMemcpy(runged_z, runged_y, runge_sizeByte, cudaMemcpyDeviceToDevice), 8) != RUNGE_SUCCESS) return RUNGE_FAIL;
 vectorAdd(N, runged_z, runged_k1, c5);
 vectorAdd(N, runged_z, runged_k3, c6);
 vectorAdd(N, runged_z, runged_k4, c7);
 vectorAdd(N, runged_z, runged_k5, c8);
 vectorAdd(N, runged_z, runged_k6, c9);
 }
 while (checkRKFError(N, &h, &s, tol, method));
 if (verbose)
 {
 printf("RKF45>> t: %10g; step: %10g\n", t, h);
 }
 if (cudaSafeCall(cudaMemcpy(runged_y, runged_z, runge_sizeByte, cudaMemcpyDeviceToDevice), 9) != RUNGE_SUCCESS) return RUNGE_FAIL;
 t += h;
 }
 if (cudaSafeCall(cudaMemcpy(y, runged_y, runge_sizeByte, cudaMemcpyDeviceToHost), 9) != RUNGE_SUCCESS) return RUNGE_FAIL;
 *h0 = h;
 }
 cudaEventRecord(stop, 0);
 cudaEventSynchronize(stop);
 cudaEventElapsedTime(time, start, stop);

 freeArrays(method);
 return RUNGE_SUCCESS;
}

////////////////////////////////////////////////////////////////////////
//
// Save call of CUDA functions
//
RungeExitCodes cudaSafeCall(cudaError_t err, TYPE num)
{
 if (err != cudaSuccess)
 {
 printf("CUDA message: %s\t\t%f\n", cudaGetErrorString(err), num);
 return RUNGE_FAIL;
 }
 return RUNGE_SUCCESS;
}

////////////////////////////////////////////////////////////////////////
//
// Initialize global arrays
//
RungeExitCodes initArrays(RungeMethod method)
{
 if (cudaSafeCall(cudaMalloc((void**)&runged_k1, runge_sizeByte), 10) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_k2, runge_sizeByte), 11) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_k3, runge_sizeByte), 12) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_k4, runge_sizeByte), 13) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_k5, runge_sizeByte), 14) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (method == RM_RKF45) if (cudaSafeCall(cudaMalloc((void**)&runged_k6, runge_sizeByte), 15) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_y, runge_sizeByte), 16) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_y1, runge_sizeByte), 17) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_z, runge_sizeByte), 18) != RUNGE_SUCCESS) return RUNGE_FAIL;
 if (cudaSafeCall(cudaMalloc((void**)&runged_ytemp, runge_sizeByte), 19) != RUNGE_SUCCESS) return RUNGE_FAIL;

 return RUNGE_SUCCESS;
}

////////////////////////////////////////////////////////////////////////
//
// Free global arrays
//
RungeExitCodes freeArrays(RungeMethod method)
{
 cudaFree(runged_k1);
 cudaFree(runged_k2);
 cudaFree(runged_k3);
 cudaFree(runged_k4);
 cudaFree(runged_k5);
 if (method == RM_RKF45) cudaFree(runged_k6);
 cudaFree(runged_y);
 cudaFree(runged_y1);
 cudaFree(runged_z);
 cudaFree(runged_ytemp);

 return RUNGE_SUCCESS;
}

////////////////////////////////////////////////////////////////////////
//
// Host function
// Addition of two N-dimensional vectors v1 and v2 placed in global
// memory (on videocard)
//     v1 := v1 + multiplier*v2
//
void vectorAdd(int N, TYPE * d_v1, TYPE * d_v2, TYPE multiplier)
{
 dim3    dimBlock(RUNGE_BLOCK_SIZE);
 dim3    dimGrid((N + RUNGE_BLOCK_SIZE - 1) / RUNGE_BLOCK_SIZE);
 addKernel<<<dimGrid, dimBlock>>>(N, d_v1, d_v2, multiplier);
}

////////////////////////////////////////////////////////////////////////
//
// Host function
// Calculate absolute values vector of differences between two
// N-dimensional vectors v1 and v2 placed in global memory (on
// videocard)
//     v1 := |v1 - v2|
//
void calcDiff(int N, TYPE * d_v1, TYPE * d_v2)
{
 dim3    dimBlock(RUNGE_BLOCK_SIZE);
 dim3    dimGrid((N + RUNGE_BLOCK_SIZE - 1) / RUNGE_BLOCK_SIZE);
 calcDiffKernel<<<dimGrid, dimBlock>>>(N, d_v1, d_v2);
}

////////////////////////////////////////////////////////////////////////
//
// Host function
// Calculate error of integration method and make conclusions for
// integration step h
//
bool checkRKFError(int N, TYPE * h, TYPE * s, TYPE tol, RungeMethod method)
{
 calcDiff(N, runged_ytemp, runged_z);
 dim3        dimBlock(RUNGE_BLOCK_SIZE);
 dim3        dimGrid((N + RUNGE_BLOCK_SIZE - 1) / RUNGE_BLOCK_SIZE);

 int N1 = N;
 while (N1 > 1)
 {
 reduceKernel<RUNGE_BLOCK_SIZE><<<dimGrid, dimBlock>>>(N1, runged_ytemp, runged_ytemp);
 N1 = (N1 + RUNGE_BLOCK_SIZE - 1) >> RUNGE_B_SIZE;
 }
 TYPE err;
 cudaMemcpy(&err, runged_ytemp, sizeof(TYPE), cudaMemcpyDeviceToHost);

 if (err < tol) return false;

 if (method == RM_RKF45)
 {
 (*s) = 0.840896 * exp(0.25 * log(tol*(*h)/err));
 (*h) *= (*s);
 if (*s < RUNGE_MIN_S) *s = RUNGE_MIN_S;
 if (*s > RUNGE_MAX_S) *s = RUNGE_MAX_S;
 }

 return true;
}

////////////////////////////////////////////////////////////////////////
//
// Kernel function
// Addition of two N-dimensional vectors
//
__global__ void addKernel(int N, TYPE * d_a, TYPE * d_b, TYPE multiplier)
{
 int i = blockDim.x * blockIdx.x + threadIdx.x;
 if (i < N) d_a[i] = d_a[i] + (d_b[i] * multiplier);
}

////////////////////////////////////////////////////////////////////////
//
// Kernel function
// Reduction: sum all elements of d_xin vector
//
template<unsigned int blockSize>
__global__ void reduceKernel(int N, TYPE * d_xin, TYPE * d_xout)
{
 __shared__ TYPE sdata[RUNGE_BLOCK_SIZE];

 size_t tid = threadIdx.x;
 size_t i = blockDim.x*2*blockIdx.x + threadIdx.x;
 size_t gridSize = blockSize*2*gridDim.x;
 sdata[tid] = 0;
 while (i < N)
 {
 if (i + blockDim.x < N) sdata[tid] += d_xin[i] + d_xin[i + blockDim.x];
 else if (i < N) sdata[tid] += d_xin[i];
 i += gridSize;
 }
 __syncthreads();

 if (blockSize >= 512)
 if (tid < 256) sdata[tid] += sdata[tid + 256]; __syncthreads();
 if (blockSize >= 256)
 if (tid < 128) sdata[tid] += sdata[tid + 128];  __syncthreads();
 if (blockSize >= 128)
 if (tid < 64) sdata[tid] += sdata[tid + 64];  __syncthreads();

 if (tid < 32)
 {
 if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
 if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
 if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
 if (blockSize >=  8 ) sdata[tid] += sdata[tid + 4];
 if (blockSize >=  4) sdata[tid] += sdata[tid + 2];
 if (blockSize >=  2) sdata[tid] += sdata[tid + 1];
 }

 if (tid == 0) d_xout[blockIdx.x] = sdata[0];
}

////////////////////////////////////////////////////////////////////////
//
// Kernel function
// Calculate absolute value of differences between vectors a and b
//
__global__ void calcDiffKernel(int N, TYPE * d_a, TYPE * d_b)
{
 int i = blockDim.x * blockIdx.x + threadIdx.x;
 if (i < N) d_a[i] = fabs(d_a[i] - d_b[i]);
}

Пример функции вычисления производной

Файл funct.cuh

#ifndef _FUNCT_H_
#define _FUNCT_H_

#include "funcdef.cuh"

void func(int N, TYPE t, TYPE * y, TYPE * deriv, TYPE h, DerivParams * params);

#endif

Файл funct.cu

#include "funct.cuh"

////////////////////////////////////////////////////////////////////////
//
// Kernel function for derivative calculation
//     \ddot y + \omega^2 y = 0
//
__global__ void getDeriv(TYPE t, TYPE * y, TYPE * deriv, TYPE h, DerivParams * params)
{
 int i = blockDim.x * blockIdx.x + threadIdx.x;
 __shared__ DerivParams sprm;
 sprm.omega = params->omega;

 if (i == 0) deriv[0] = y[1] * h;
 if (i == 1) deriv[1] = - sprm.omega * sprm.omega * y[0] * h;
}

////////////////////////////////////////////////////////////////////////
//
// Host function for derivative calculation
//
void func(int N, TYPE t, TYPE * y, TYPE * deriv, TYPE h, DerivParams * params)
{
 dim3    dimBlock(RUNGE_BLOCK_SIZE);
 dim3    dimGrid((N + RUNGE_BLOCK_SIZE - 1) / RUNGE_BLOCK_SIZE);
 getDeriv<<<dimGrid, dimBlock>>>(t, y, deriv, h, params);
}

main.cpp

#include <cstdio>
#include <cstdlib>

#include "funct.cuh"
#include "runge_cuda.cuh"

int main(int argc, char ** argv)
{
 int N = 2;
 TYPE y0[N]; y0[0] = 0; y0[1] = 1;
 TYPE y[N];
 TYPE t0 = 0;
 TYPE t1 = 3.1415926/2;
 TYPE h0 = 5;
 TYPE tol = 1e-5;
 bool verbose = true;
 DerivParams params;
 params.omega = 1;
 TYPE time;

 rungeCuda(N,
 func,
 y0,
 y,
 t0,
 t1,
 &h0,
 tol,
 verbose,
 &params,
 &time,
 RM_RKF45);

 printf("Result: %f Time: %f milliseconds\n", y[0], time);

 return 0;
}

Makefile

NVCXX := nvcc
CXX := g++
CXXFLAGS := -I/usr/include/cuda
NVFLAGS :=
LIBRARIES := -lm -lcuda
OUTPUT := rkf

$(OUTPUT): main.o funct.o runge_cuda.o
 $(NVCXX) main.o funct.o runge_cuda.o -o $(OUTPUT) $(LIBRARIES)

main.o: main.cpp
 $(CXX) $(CXXFLAGS) -c main.cpp -o main.o

funct.o: funct.cu
 $(NVCXX) $(CXXFLAGS) -c funct.cu -o funct.o

runge_cuda.o: runge_cuda.cu
 $(NVCXX) $(NVFLAGS) -c runge_cuda.cu -o runge_cuda.o

clean:
 rm -rf *.o *.linkinfo $(OUTPUT)
Реклама