一种使用gpu加速lcs计算的方法
【专利摘要】一种使用GPU加速LCS计算的方法，步骤如下：(1)在CPU中，在给定的庞加莱界面上进行二维网格划分，网格划分个数为N×N，分配GPU的常量内存，设置GPU的线程块的数量为1，线程块的线程为二维线程，将每个网格点对应一个二维线程；(2)确定所有网格点的初始状态参数，将所有网格点的初始状态参数复制至上述GPU的常量内存；(3)设计变步长4阶龙格库塔积分器，利用计分器计算每个网格点的有限时间李雅普诺夫时间指数；(4)将所有网格点的有限时间李雅普诺夫时间指数传回CPU，利用回传的数据绘制二维云图获得有限时间李雅普洛夫指数域进而提取LCS。
【专利说明】
一种使用GPU加速LCS计算的方法
技术领域
[0001]本发明涉及一种计算拉格朗日拟序结构（LCS)的方法，更特别地说，是采用GPU通 用计算平台开发一种用于拉格朗日拟序结构计算的并行算法。
【背景技术】
[0002] 在平面圆型限制性三体问题(PCR3BP)中，在其共线平动点附近周期轨道的不变流 形，以渐进性质和运动分界面的特点被广泛用于小推力转移、地月低能轨道转移以及编队 飞行等空间任务。然而对于限制性四体问题等时间相关的动力系统，由于无法摆脱对时间 的显性依赖，相关研究难以推广。
[0003] 另一方面，拉格朗日拟序结构，定义为有限时间李雅普诺夫指数域中的脊线。有限 时间李雅普诺夫指数，指在相流作用下，邻近两点在有限时间内的平均最大分离速率，度量 了系统对于初值的敏感性。由动力系统理论可知，不变流形位于不同测的邻近两点会以指 数速率快速分离，即不变流形对应于较大的李雅普诺夫指数，从而拉格朗日拟序结构(LCS) 可以定义为有限时间李雅普诺夫指数域中的脊线。相流向后时，LCS蕴含稳定流形，向前时 蕴含不稳定流形。因此LCS可以作为时变流形场中的不变流形来研究动力系统的相空间，实 现在非自治系统（日-地-月-航天器系统）中对地月低能转移轨道的直接构建。然而，计算 LCS时，网格划分的越细，LCS越精确的同时计算也越耗时，用CPU设计串行的计算算法精确 绘制LCS计算量大，耗时长，这极大的限制了 LCS在实际轨道任务设计中的应用。北京航空航 天大学的祁瑞在论文《限制性四体问题中的时间相关不变流形》中，利用时间相关不变流形 的分界面性质，采用二分法对LCS进行提取，较粗的网格划分就能确定LCS的内部和外部，较 大的提高了计算效率。然后，其算法本质上还是串行算法，没有充分利用LCS "易并行化"的 计算特点，无法从本质上大幅减小计算时间，无法使LCS成为一个轨道任务设计的便利工 具。
【发明内容】
[0004] 本发明充分利用LCS计算易于并行化的特点，设计一种基于GPU的并行算法，解决 四体问题下的LCS计算耗时太长的问题。
[0005] 本发明的技术解决方案是:一种使用GPU加速LCS计算的方法，步骤如下：
[0006] (1)在CPU中，在给定的庞加莱界面上进行二维网格划分，网格划分个数为N X N，分 配GPU的常量内存，设置GPU的线程块的数量为1，线程块的线程为二维线程，将每个网格点 对应一个二维线程；
[0007] (2)确定所有网格点的初始状态参数，将所有网格点的初始状态参数复制至上述 GHJ的常量内存；
[0008] (3)设计变步长4阶龙格库塔积分器，利用计分器计算每个网格点的有限时间李雅 普诺夫时间指数；
[0009] (4)将所有网格点的有限时间李雅普诺夫时间指数传回CPU，利用回传的数据绘制 二维云图获得有限时间李雅普洛夫指数域进而提取LCS。
[0010] 所述设计变步长4阶龙格库塔积分器的步骤如下：
[0011] (3.1)建立在GPU上运行的矩阵类matrix，用于实现有限时间李雅普诺夫时间指数 计算过程中所有矩阵的运算；
[0012] (3.2)设置龙格库塔法计算过程中的中间变量matrix ki和状态变量matrix yn;
[0013] (3.3)分别以步长11 = 8丨6口和11 = 0.58丨6口为步长，计算丨+11时刻的状态变量;1:为当 前时刻；
[0014] (3.4)判断（3.3)中两个步长下的状态变量的矢量差模量是否小于预设的精度要 求，若小于，则将状态转移到t+h时刻，并将步长h加倍，重复步骤(3.3)直至积分时间结束； 否则，步长h减半，重复步骤(3.3)。
[0015]所述的矩阵类matrix满足如下条件：能在GPU上运行，使用宏定义的方式，不用或 少用动态分配内存;能完成LCS计算中使用的矩阵基本运算如加减乘、求逆、转置和求特征 值。
[0016] 所述步骤(4)中提取LCS的步骤如下：
[0017]已知LCS定义为李雅普诺夫域脊线，局部最大值的点即属于LCS;基于此，对LCS上 的点按行提取，假设LCS网格数为NXN，则每一行有两个点属于LCS，以提取第i行点为例，从 左边遍历得首个最大值为点ii，从右边遍历得首个最大值为点i2,则所有行点ii，i2分别相 连便可提取LCS。
[0018]本发明与现有技术相比有益效果为：
[0019] (1)由于目前国内外对LCS在航天器轨道设计领域处于起步阶段，LCS超大计算量 极大的限制了其应用，专门针对LCS优化计算的研究较少。结合LCS计算的易并行化特点，本 发明首次提出利用GPU完成LCS计算的算法。仿真试验结果表明，使用本发明算法，与使用传 统CPU相比，可获得近30倍的加速比。以400 X 400网格的LCS计算为例，传统CPU计算需要30 个小时左右才能完成，而在GPU设备Tesla C2075上只需要1个小时左右。（2)由于有限时间 李雅普诺夫指数的计算需要求解20X1维的状态转移矩阵微分方程，计算量相当大，而且涉 及较多的矩阵运算，矩阵类matrix的使用使得GPU上的LCS算法开发工作量大为减小，同时 也较好的完成了GPU上的内存管理。
[0020] (3)本发明的LCS提取算法摆脱了传统依赖有限时间李雅普诺夫指数域云图的方 法，使用云图必须依靠肉眼观察有限时间李雅普诺夫指数域的脊线或者使用图像处理获 得，而使用本发明可以提取LCS数据进而绘制LCS图。
【附图说明】
[0021 ]图1为本发明方法流程图；
[0022]图2为本发明实施例得到的有限时间李雅普诺夫指数域；
[0023]图3为本发明实施例在给定庞加莱界面上有限时间李雅普洛夫指数域中提取的 LCS0
【具体实施方式】
[0024]如下面将结合附图和实施例对本发明做进一步的详细说明。
[0025] 如图1所示，算法主题部分在GPU上完成，CPU负责GPU的线程配置，内存分配工作并 将GPU传回的有限时间李雅普诺夫指数进行精确提取从而获得拉格朗日拟序结构。
[0026] 开发基于GPU的变步长4阶龙格库塔积分器、在给定的庞加莱截面上进行二维网格 划分、根据二维网格特点，对应使用GPU并行计算架构CUDA的二维线程实现线程分配、对每 个二维线程上的二维网格点计算有限时间李雅普诺夫指数、精确提取LCS。
[0027]步骤一:CPU上，给定庞加莱截面
[0028]
[0029] 这里?为y轴方向速度大小，H为初始时刻系统能量。
[0030] 步骤二:确定每个格点的初始状态:假设网格划分个数为NXN，由于这里没有利用 GPU共享缓存必要，因此线程块数量DimGrid设置为1，而每个线程块的线程为二维线程，用 CUDA编程变量表示为：
[0031] dim3threads(N,N)
[0032] 确定好网格划分后，在给定的庞加莱截面山上，可以获得所有网格点的初始状态 参数：
[0033]
[0034] 其中i, j = l，2, . . .N，疗(AJ)为系统的等效势函数。
[0035]步骤三：由于把变量数据存入GPU常量内存后，将把变量的访问限制为只读，这样 与普通的全局内存中读数据相比，从常量内存中读取相同的数据可以节约内存的带宽约1/ 16,因此选择将所有网格点的初始状态参数复制至GPU上的常量内存。
[0036] 步骤四：计算每个网格点的李雅普诺夫时间指数：
[0037] 有限时间李雅普诺夫指数定义为
[0038] λ
/
[0039：
为to到t时刻的状态转移矩阵，可以通过状态转移矩阵微分方 程
[004C
[0041] 数值积分获得，其中
％系统的雅克比矩阵，数值积分采用 专门开发用于GHJ平台的变步长4阶龙格库塔积分器。
[0042] 所述设计变步长4阶龙格库塔积分器的步骤如下：
[0043] (3.1)建立在GPU上运行的矩阵类matrix，用于实现有限时间李雅普诺夫时间指数 计算过程中所有矩阵的运算；
[0044] (3.2)设置龙格库塔法计算过程中的中间变量matrix ki和状态变量matrix yn;
[0045] (3.3)分别以步长11 = 8丨6口和11 = 0.58丨6口为步长，计算丨+11时刻的状态变量;1:为当 前时刻；
[0046] (3.4)判断（3.3)中两个步长下的状态变量的矢量差模量是否小于预设的精度要 求，若小于，则将状态转移到t+h时刻，并将步长h加倍，重复步骤(3.3)直至积分时间结束； 否则，步长h减半，重复步骤(3.3)。
[0047]步骤五:将所有网格点的有限时间李雅普诺夫指数传回CPU，LCS上的点按行提取， 假设LCS网格数为NXN，则每一行有两个点属于LCS，以提取第i行点为例，从左边遍历得首 个最大值为点ii，从右边遍历得首个最大值为点i 2,则所有行点h，i2分别相连便可提取 LCS0
[0048] 实施例1
[0049] 在实施例中，采用本发明所提出的使用GPU加速LCS计算的方法进行LCS计算。
[0050] 步骤一:CPU上，给定庞加莱截面
[0051] i7I ^{(6)' i') I == -0.01215,// ===： -1.571722820493,0.2 < j < 0.6, -0.6 < y < O.l}. ?
[0052] 步骤二:确定每个格点的初始状态:这里网格划分个数为400 X 400，线程块大小设 置为1，每个线程块的二维线程大小设置为虹!113也^&(18(400,400)，其中第一维表示 7，第二 位表示>。
[0053]确定好网格划分后，在给定的庞加莱截面山上，可以获得所有网格点的初始状态 参数：
[0054]
[0063] 1、首先给出矩阵类和龙格库塔积分器的构建示例，其总的要求是能在GPU运行并 尽量少用内存动态分配，减小算法开发难度。
[0064]矩阵类的构建如下所示：
[0065] -deviee- class matrix public： -device-matrix(iiUrow，imcol);//，构造函数，创建矩昨，分?矩昨数 ^device - matrix (): /7，析构闲数，销毁矩阵，释放矩昨数辦内// -device____ matrix(const niatrix&); //?制构造函数，执行矩丨昨深度S制 __device_ matrix& opcrator=(const matrix&); /7这里及以丨、*几个两数市载运算符 -device- friend const mairix opcialor^iconst matrix&, const matrix&); device- friend consi mairix operator十(const malriXfc, consi tnalrix&); -device- !ricnd const matrix opcralor-(con^l rna1rix&, const mairix&); _ device_ friend const matrix 〇perator/(c〇nst matrix&!, const matrix&); __device_ ft lend const double norm(consi matrix&); //计界矩昨的 2 范数 protected: hu mRuvv; //矩阵行数 ini mCol; //矩昨列数 double* mAmiy; /7Y/M矩阵数据，按列存傭 }
[0066] 2:基于矩阵类的GPU上的龙格库塔积分算法开发，其原型如下，
[0067] -device--matrix rungekutta(matrix(*f)(double t,const matrix&),matrix t，matrix x0，double tol)；
[0068] 其中matrix(*pf)(double t，const matrix&)为所要求解的微分方程；
[0069] matrix t为积分时间；
[0070] matrix x0为积分初值；
[0071 ] double tol为要求最低积分精度；
[0072]结果返回积分状态参数。使用了矩阵类后的变步长龙格库塔积分器开发大为简 化，具体步骤如下所示，
[0073] 步骤1 · 1 设置初始步长double step = le_8;
[0074] 步骤1.2设置中间变量11^1^1^1，1^2，1^3，1^4和状态变量11^1^711;
[0075] 步骤1.3分别以步长step，和0.5step计算t = t+step后的状态值^二和）^
[0076] kl=f(t,yn)；
[0077] k2 = f(t+0.5*step,yn+0.5*step*kl);
[0078] k3 = f(t+0.5*step,yn+0.5*step*k2);
[0079] k4 = f(t+0.5*step，yn+step*k3);
[0080] ynl=yn+step/6*(kl+2*k2+2*k3+k4);
[0081 ]步骤1.4比较两个状态值的差值模量是否小于精度要求tol，
[0082]
[0083] 状态转移到丨=丨+8丨6口;711 = 7111;步长加倍8丨6口 = 2*8丨6口;如果积分时间到终点， 停止计算，否者转到步骤1.2.3继续求解；
[0084]
、步长减半8丨6口 = 8丨6口/2，重新回到步骤1.2.3计算。
[0085] 步骤五:将所有网格点的有限时间李雅普诺夫指数传回CPU并提取LCS，其结果如 图2、3所示，其中图2为有限时间李雅普诺夫指数域云图，图3为从图2提取的LCS，比较两图 可知，LCS把运动的分界面精确的绘制出来。计算总共耗时3677.408s。
[0086] 本发明未详细说明部分属于本领域技术人员公知常识。
【主权项】
1. 一种使用GHJ加速LCS计算的方法，其特征在于步骤如下： (1) 在CPU中，在给定的庞加莱界面上进行二维网格划分，网格划分个数为NXN，分配 GPU的常量内存，设置GPU的线程块的数量为1，线程块的线程为二维线程，将每个网格点对 应一个二维线程； (2) 确定所有网格点的初始状态参数，将所有网格点的初始状态参数复制至上述GPU的 常量内存； (3) 设计变步长4阶龙格库塔积分器，利用计分器计算每个网格点的有限时间李雅普诺 夫时间指数； (4) 将所有网格点的有限时间李雅普诺夫时间指数传回CPU，利用回传的数据绘制二维 云图获得有限时间李雅普洛夫指数域进而提取LCS。2. 根据权利要求1所述的一种使用GPU加速LCS计算的方法，其特征在于:所述设计变步 长4阶龙格库塔积分器的步骤如下： (3.1) 建立在GPU上运行的矩阵类matrix，用于实现有限时间李雅普诺夫时间指数计算 过程中所有矩阵的运算； (3.2) 设置龙格库塔法计算过程中的中间变量matrix ki和状态变量matrixyn; (3.3) 分别以步长h = step和h = 0.5step为步长，计算t+h时刻的状态变量;t为当前时 刻； (3.4) 判断(3.3)中两个步长下的状态变量的矢量差模量是否小于预设的精度要求，若 小于，则将状态转移到t+h时刻，并将步长h加倍，重复步骤(3.3)直至积分时间结束;否则， 步长h减半，重复步骤(3.3)。3. 根据权利要求2所述的一种使用GPU加速LCS计算的方法，其特征在于:所述的矩阵类 matrix满足如下条件:能在GPU上运行，使用宏定义的方式，不用或少用动态分配内存;能完 成LCS计算中使用的矩阵基本运算如加减乘、求逆、转置和求特征值。4. 根据权利要求1所述的一种使用GPU加速LCS计算的方法，其特征在于:所述步骤(4) 中提取LCS的步骤如下： 已知LCS定义为李雅普诺夫域脊线，局部最大值的点即属于LCS;基于此，对LCS上的点 按行提取，假设LCS网格数为NXN，则每一行有两个点属于LCS，以提取第i行点为例，从左边 遍历得首个最大值为点ii，从右边遍历得首个最大值为点i2,则所有行点ii，i2分别相连便 可提取LCS。
【文档编号】G06F19/00GK106055890SQ201610363856
【公开日】2016年10月26日
【申请日】2016年5月26日
【发明人】徐 明, 林名培, 贾向华, 王召辉, 马跃辰, 潘晓
【申请人】北京航空航天大学