【发布时间】：2017-02-03 16:01:37
【问题描述】：
我实际上正在使用 CUDA，并且正在尝试使用该技术优化程序。所以我有一个大内核，我必须在 100k+ 到 100M+ 时间之间启动，也许数十亿？
所以我使用 dim3 变量读取允许启动该数量的线程（参见：https://devtalk.nvidia.com/default/topic/621867/size-limitation-for-1d-arrays-in-cuda-/?offset=7）
我有一个示例代码（在我的 gtx970 上）有时会运行，有时不会。
#ifndef PROPAGATORSAT_CUH_
# define PROPAGATORSAT_CUH
# define M_PI (3.14159265358979323846)
# define TWO_PI (2 * M_PI)
# define TOTAL_TIME (615359.772)
# define STEP (0.771)
# define NB_IT (TOTAL_TIME / (double)STEP)
# define NB_THREADS (1024)
# define NB_BLOCKS (int)((NB_IT + NB_THREADS - 1) / NB_THREADS)
# include <cmath>
# include <cfloat>
# include <stdio.h>
# include "../common/book.h"
# include "cuda_runtime.h"
# include "device_launch_parameters.h"
# define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true)
{
if (code != cudaSuccess)
{
fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
if (abort) exit(code);
}
}
class Global
{
public:
const double _ITURadEarth = 6378145.0;
const double _ITUGravCst = 3.986012E5;
const double _ITUJ2 = 0.001082636;
const double _J2000AngleDeg = 0;//-79.8058;
const double _J2000AngleRad = 0;//TO_RAD(_J2000AngleDeg);
const double _ITUAngleRateEarthRot = 4.1780745823E-3;
const double _ITUAngleRateEarthRotRad = degToRad(_ITUAngleRateEarthRot);
public:
__device__ double myAsin(double angle);
__host__ __device__ double myAcos(double angle);
__host__ __device__ double negPiToPi(double angle);
__host__ __device__ double degToRad(double angle);
__device__ double radToDeg(double angle);
};
class Cartesian
{
public:
double _X;
double _Y;
double _Z;
private:
double _m;
public:
__host__ __device__ Cartesian(double x, double y, double z) : _X(x), _Y(y), _Z(z), _m(-1) {}
};
class Propagator
{
public:
double _iDeg;
double _a;
double _omega_0;
double _OMEGA_0;
double _omega_r;
double _OMEGA_r;
double _rho;
double _SinI;
double _CosI;
double _p;
double _e;
double _ReKm;
double _n0;
double _n_bar;
double _M0;
double _sqrt_e;
int _orbitCase = -1;
double _WdeltaRad;
double _precessionRateRad;
double _artificialPrecessionRad = DBL_MIN;
double _simulationDuration = DBL_MIN;
double _incrementWdeltaRad;
void propagator(double smaKm,
double incDeg,
double e,
double raanDeg,
double aopDeg,
double trueAnomalyDeg,
bool stationKeeping,
double WdeltaDeg,
bool precessionMechanismSupplied,
double precessionRateDeg);
__device__ Cartesian evaluate(double timeSec,
double simulationDuration,
double artificialPrecessionRad,
bool ECImode);
__device__ double solveKepler(double M,
double e,
double epsilon);
__device__ Cartesian rotateOrbitalElements(Cartesian pq0,
double omega,
double OMEGA,
double CosI,
double SinI);
};
#endif /* !PROPAGATORSAT_CUH_ */
__host__ __device__ double Global::myAcos(double angle)
{
return (acos(((angle > 1) ? (1) : (angle < -1) ? (-1) : (angle))));
}
__device__ double Global::myAsin(double angle)
{
return (asin(((angle > 1) ? (1) : (angle < -1) ? (-1) : (angle))));
}
__host__ __device__ double Global::degToRad(double angle)
{
return (angle * M_PI / 180.0);
}
__device__ double Global::radToDeg(double angle)
{
return (angle * 180.0 / M_PI);
}
__host__ __device__ double Global::negPiToPi(double angle)
{
double output;
output = fmod(angle, TWO_PI);
output = fmod(angle + TWO_PI, TWO_PI);
return ((output > M_PI) ? (output - TWO_PI) : (output));
}
void Propagator::propagator(double smaKm, double incDeg, double e, double raanDeg, double aopDeg, double trueAnomalyDeg, bool stationKeeping, double WdeltaDeg, bool precessionMechanismSupplied, double precessionRateDeg)
{
double iRad, trueAnomalyRad, cosV, E, mu;
Global global;
_iDeg = incDeg;
iRad = global.degToRad(_iDeg);
_CosI = cos(iRad);
_SinI = sin(iRad);
_e = e;
_a = smaKm;
trueAnomalyRad = global.degToRad(trueAnomalyDeg);
if (e == 0)
_M0 = trueAnomalyRad;
else
{
cosV = cos(trueAnomalyRad);
E = global.myAcos((e + cosV) / (1 + e * cosV));
if (global.negPiToPi(trueAnomalyRad) < 0)
E = M_PI * 2 - E;
_M0 = E - e * sin(E);
}
_OMEGA_0 = global.degToRad(raanDeg);
_omega_0 = global.degToRad(aopDeg);
_p = _a * (1 - e * e);
_ReKm = global._ITURadEarth / 1000;
mu = global._ITUGravCst;
_n0 = sqrt(mu / pow(_a, 3));
_n_bar = _n0 * (1.0 + 1.5 * global._ITUJ2 * pow(_ReKm, 2) / pow(_p, 2) * (1.0 - 1.5 * pow(_SinI, 2)) * pow(1.0 - pow(e, 2), 0.5));
_OMEGA_r = -1.5 * global._ITUJ2 * pow(_ReKm, 2) / pow(_p, 2) * _n_bar * _CosI;
_omega_r = 1.5 * global._ITUJ2 * pow(_ReKm, 2) / pow(_p, 2) * _n_bar * (2.0 - 2.5 * pow(_SinI, 2));
_sqrt_e = sqrt((1 + e) / (1 - e));
_WdeltaRad = global.degToRad(WdeltaDeg);
_precessionRateRad = global.degToRad(precessionRateDeg);
if (stationKeeping == false)
_orbitCase = 1;
else if (precessionMechanismSupplied == false)
_orbitCase = 2;
else
_orbitCase = 3;
}
__device__ Cartesian Propagator::rotateOrbitalElements(Cartesian pq0, double omega, double OMEGA, double CosI, double SinI)
{
double CosOMEGA, SinOMEGA, CosOmega, SinOmega, R11, R12, R13, R21, R22, R23, R31, R32, R33, x, y, z;
CosOMEGA = cos(OMEGA);
SinOMEGA = sin(OMEGA);
CosOmega = cos(omega);
SinOmega = sin(omega);
R11 = CosOMEGA * CosOmega - SinOMEGA * SinOmega * CosI;
R12 = -CosOMEGA * SinOmega - SinOMEGA * CosOmega * CosI;
R13 = SinOMEGA * SinI;
R21 = SinOMEGA * CosOmega + CosOMEGA * SinOmega * CosI;
R22 = -SinOMEGA * SinOmega + CosOMEGA * CosOmega * CosI;
R23 = -CosOMEGA * SinI;
R31 = SinOmega * SinI;
R32 = CosOmega * SinI;
R33 = CosI;
x = R11 * pq0._X + R12 * pq0._Y + R13 * pq0._Z;
y = R21 * pq0._X + R22 * pq0._Y + R23 * pq0._Z;
z = R31 * pq0._X + R32 * pq0._Y + R33 * pq0._Z;
Cartesian cart = Cartesian(x, y, z);
return (cart);
}
__device__ Cartesian Propagator::evaluate(double timeSec, double simulationDuration, double artificialPrecessionRad, bool ECImode = true)
{
double M, E, v, cosV, sinV, rotationAngleECF, omega, OMEGA;
Global global;
if (_simulationDuration != simulationDuration || _artificialPrecessionRad != artificialPrecessionRad)
{
_simulationDuration = simulationDuration;
_artificialPrecessionRad = artificialPrecessionRad;
_incrementWdeltaRad = (_WdeltaRad * 2) / _simulationDuration;
}
M = _M0 + ((_orbitCase == 3) ? _n0 : _n_bar) * timeSec;
E = E = (_e == 0) ? M : solveKepler(M, _e, 1e-8);
v = 2.0 * atan(_sqrt_e * tan(E / 2));
cosV = cos(v);
sinV = sin(v);
_rho = _p / (1 + _e * cosV);
rotationAngleECF = (ECImode) ? 0 : -1 * (global._J2000AngleRad + timeSec * global._ITUAngleRateEarthRotRad);
omega = _omega_0 + ((_orbitCase == 3) ? 0 : _omega_r * timeSec);
OMEGA = _OMEGA_0 + rotationAngleECF + ((_orbitCase == 3) ? 0 : _OMEGA_r * timeSec);
if (_orbitCase == 1)
OMEGA += artificialPrecessionRad * timeSec;
else if (_orbitCase == 2)
OMEGA += _WdeltaRad * ((2.0 * timeSec / _simulationDuration) - 1);
else if (_orbitCase == 3)
OMEGA += _precessionRateRad * timeSec - _WdeltaRad + _incrementWdeltaRad * timeSec;
Cartesian pq0 = Cartesian(1000 * _rho * cosV, 1000 * _rho * sinV, 0);
Cartesian positionECI = Propagator::rotateOrbitalElements(pq0, omega, OMEGA, _CosI, _SinI);
return (positionECI);
}
__device__ double Propagator::solveKepler(double M, double e, double epsilon)
{
double En, Ens;
En = M;
Ens = En - (En - e * sin(En) - M) / (1 - e * cos(En));
while (abs(Ens - En) > epsilon)
{
En = Ens;
Ens = En - (En - e * sin(En) - M) / (1 - e * cos(En));
}
return (Ens);
}
__global__ void kernel(Propagator *CUDA_prop)
{
size_t tid;
tid = (blockIdx.x + blockIdx.y * gridDim.x + gridDim.x * gridDim.y * blockIdx.z) * blockDim.x + threadIdx.x;
//if (tid < NB_IT)
Cartesian positionNGSOsatECI = CUDA_prop[0].evaluate(STEP * tid, 615359.772, 0);
}
int main(void)
{
cudaEvent_t start, stop;
HANDLE_ERROR(cudaEventCreate(&start));
HANDLE_ERROR(cudaEventCreate(&stop));
HANDLE_ERROR(cudaEventRecord(start, 0));
Propagator prop[1], *CUDA_prop;
dim3 block(1000, 1, 1);
dim3 thread(1024, 1, 1);
prop[0].propagator(7847.3, 53, 0, 18, 0, 67.5, true, 5, true, 3.4000000596279278E-05);
HANDLE_ERROR(cudaMalloc((void **)&CUDA_prop, sizeof(Propagator)));
HANDLE_ERROR(cudaMemcpy(CUDA_prop, prop, sizeof(Propagator), cudaMemcpyHostToDevice));
kernel <<< block, thread >>> (CUDA_prop);
gpuErrchk(cudaPeekAtLastError());
gpuErrchk(cudaDeviceSynchronize());
HANDLE_ERROR(cudaFree(CUDA_prop));
HANDLE_ERROR(cudaEventRecord(stop, 0));
HANDLE_ERROR(cudaEventSynchronize(stop));
float elapsedTime;
HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
printf("time : %f ms\n", elapsedTime);
HANDLE_ERROR(cudaEventDestroy(start));
HANDLE_ERROR(cudaEventDestroy(stop));
return (0);
}
如果我启动那么多“线程”？它一直工作到大约 300k 块。但有时对于相同的数量它不起作用。我收到一条错误消息：来自该行的“未知错误”：
gpuErrchk(cudaDeviceSynchronize());
或 cudaFree，或内核调用后的某个函数。 在此处输入代码
如果我只使用 1k 块和 1k 线程启动并使用 cuda-memcheck，我会得到与以前相同的错误，但如果没有 cuda-memcheck，它运行得很好。
我不知道是什么原因导致这个问题以及如何解决它
NB : HANDLE_ERROR 宏可以通过 gpuErrchk 宏来更改，它是一个库中的定义，可以做完全相同的事情
我还想知道如何确定我可以使用硬件或任何规格的规格启动的最大线程数。
【问题讨论】：
-
听起来很像watchdog timer 的作用。您的内核运行了多长时间才会出现故障？
-
您的 GPU 是专用的计算设备还是您也将其用于显示？
-
我还使用我的 GPU 进行显示，内核运行 300k 块和 1024 个线程的时间约为 2100 毫秒，100k 块和 1024 个线程的运行时间约为 700 毫秒。对于看门狗，我看起来更早，它设置为 6 或 7 秒。当它失败时，它不会运行超过 3 秒
-
听起来更像看门狗定时器。你在 Windows 上吗？您是否连续启动多个内核？ Windows 驱动程序默认批处理内核以减少启动开销。您可以通过在两次启动之间放置
cudaStreamQuery(0)来强制单独启动。 -
与其增加更长内核的超时时间，不如尝试启动更多内核。这将允许 GUI 在两者之间进行更新，这就是拥有看门狗计时器的动机。
标签： cuda