采用CUDA Thrust实现的空间坐标系变换

采用CUDA Thrust实现的空间坐标系变换

系列文章目录

基于CUDA的空间目标轨道并行计算技术，共四节，其中第一、二节目录如下

1 基于CUDA的空间目标轨道计算需求与任务分析
2 基于CUDA的空间坐标系变换矩阵计算
2.1 基于SOFA的空间坐标系变换
2.2 采用CUDA Thrust实现的空间坐标系变换
2.3 采用CUDA Runtime实现的空间坐标系变换
————————————————

---

原理

CUDA支持基于Thrust进行开发，此时仅需实现对应的仿函数即可，但是不同于runtime模式，此时是通过序号来区分线程的。

基于CUDA的空间坐标系变换，实际上是并行计算各个时间点的变换矩阵。每个线程输入时间，输出所需变换矩阵。

由于Thrust的输入是序号（效率最高），因此需根据序号计算时间，为此在仿函数中需提供根据序号计算时间所需的成员变量。

同时，极移等数据是预先在CPU端加载到内存的，而不能GPU端从内存中读取数据，因此需要提供访问机制。

因此仿函数类定义为：

struct c2tTransFunctor	{
	timeOfSpace _bt;	//开始时间
	int		_stepMilliSecond;//以毫秒为单位的步长
	double	*_xps, *_yps, *_dut1s, *_ddp80s, *_dde80s;//极移数据地址
	int  _eopNum;	//极移数据个数
	__device__	mat3x3 operator()(int& t) {
		int off = (long long int)_stepMilliSecond * t / (86400 * 1000);
		if (off >= _eopNum) off = _eopNum - 1;
		double xp = _xps[off] * DAS2R;
		double yp = _yps[off] * DAS2R;
		double dut1 = _dut1s[off];
		timeOfSpace t0 = _bt;
		t0.addMilliSecond((long long int)_stepMilliSecond*t);
		iauCal2jd(t0._year, t0._month, t0._day, &djmjd0, &date);
		……
		return mat;	};

仿函数类中，定义了预推的开始时间_bt和步长_stepMiliSeond（ms），在函数体中，根据上述信息和序号，即可计算得到时间（t0）。

仿函数中定义了极移数据指针*_xps, *_yps, *_dut1s, *_ddp80s, *_dde80s，根据时间可以计算对应极移数据在极移数组中的偏移。

省略的代码即是上一节中基于SOFA库中实现所对应代码。

可以推算，仿函数中定义的成员变量均将映射到GPU中的常量内存。在Thrust模式下，上述成员变量的赋值也非常简单。

void  InitC2TFunctorPolarData(const timeOfSpace& t0, const timeOfSpace& t1, c2tTransFunctor& functor){
	thrust::host_vector<double> xps;
	timeOfSpace t = t0;
	while (!(t > t1))	{
		double xp, yp, dut1, ddp80, dde80;
		find(t._year, t._month, t._day, xp, yp, dut1, ddp80, dde80);
		xps.push_back(xp);
		t.addDays(1);	}
	functor._eopNum = xps.size();
	thrust::device_vector<double> d_xps = xps;
	functor._xps = thrust::raw_pointer_cast(d_xps.data());	}

按照Thrust规范，进行函数调用即可完成基于空间坐标系变换的并行计算。

propagateTransMat_Thrust(const timeOfSpace& t0, const timeOfSpace& t1, int stepMilliSecond, thrust::device_vector<mat3x3>& d_mats)	{
	struct c2tTransFunctor functor;
	functor._bt = t0;
	functor._stepMilliSecond = stepMilliSecond;
	InitC2TFunctorPolarData(t0, t1, functor);
	int num = getSamplesNums(t0,t1,stepMilliSecond) ;
	thrust::device_vector<int> d_tms(num);
	d_mats.resize(num);
	thrust::sequence(d_tms.begin(), d_tms.end(), 0, 1);
	thrust::transform(d_tms.begin(), d_tms.end(), d_mats.begin(), functor);	}

效率分析

在不同采样点的情况下，分别对CPU版本和Thrust版本进行计时，结果如表所示。

采样点数	CPU版本耗时（s）	Thrust版本耗时（s）
14400	0.155	0.079
86400	0.976	0.132
864000	5.611	0.432
864000*5	26.798	1.909

可以看出，当采样点数比较少时，并行计算的优化并不明显，而随着采样点数的增加，其优化可以达到十几倍，效率提升比较明显。

用nvpf对其进行分析（14400采样点），时间线如图

可以看出，cudaMalloc和driverAPI的调用占用的时间相当之高，远高于kernel所耗时间。同时随着采样点数的增加，cudaMalloc耗时并非线性增长，而是比较缓慢增长；另外测试表明，cudaMalloc的调用位置并不影响其耗时。因此，可以合理推断，其实cudaMalloc和driverAPI的调用所消耗时间应该包括了GPU初始化所需耗时。因此，在后续分析中，效率分析将主要针对kernel进行。

采用Thrust进行开发具有简单易用的优点，很多CUDA调用得以隐藏，在很多情况下效率与Runtime模式也差不多。

但是Runtime方式还是有其优点：一是采用Runtime方式可控性更强，可以对线程数目、访存等进行设置；二是Thrust通过仿函数返回值方式生成结果，对于生成多个计算结果的情况，只能采用结构体的方式，这样形成的结果数据是结构体数组，但是后续计算需要的往往是数组的结构，其效率会更高。因此后续重点对RunTime方式进行实现和优化。

完整代码链接

完整的源代码（结合到二体模型实现）已上传，链接为源代码