基于异构计算的大规模线性规划问题求解方法与流程

1.本发明属于高性能计算技术领域，具体涉及一种基于异构计算的大规模线性规划问题求解方法。


背景技术：

2.随着以gpu为代表的各类硬件加速卡技术的不断成熟，异构计算模式已经成为求解大规模计算问题的重要手段之一。异构计算需要协调包括cpu在内的各类计算设备(如gpu、fpga、dsp等)，将计算任务按照其特点分配到不同的计算设备上，充分发挥各种计算设备的不同特点，协同完成问题的求解。
3.在异构计算领域中，开放计算语言(open computing language,opencl)是一种能够屏蔽各类计算设备的异构性，提供统一接口供用户方便地编写异构计算程序的语言。目前opencl已经得到主流厂商的广泛支持，在服务器环境、桌面环境和嵌入式环境中都已得到广泛应用。
4.opencl平台模型如图1所示。opencl将计算机定义为主机(host)+opencl设备，其中主机一般指cpu+内存，opencl设备指opencl对各类异构计算设备的统一抽象。在计算设备内部，opencl通过计算单元(compute unit)实现对各个计算内核(kernel)的管理，计算任务由计算单元中的内核完成，每个内核执行一个线程。执行一个计算任务时，主机首先编译opencl程序，然后将待计算的数据和编译好的opencl程序通过总线发送到设备上，由设备执行指定的计算任务。在此过程中，主机与设备以及设备的不同存储区域之间可能需要完成多次数据交互。计算任务完成后，主机从设备获取计算结果并返回给最终用户。
5.除了将计算设备进行抽象外，opencl还针对计算过程中数据的传输问题设计了统一的存储器模型，如图2所示。opencl将计算设备的存储器划分为4个逻辑空间：私有存储器、局部存储器、常数存储器和全局存储器：
6.·
私有存储器：设备上每个内核都拥有的寄存器，只能被当前线程访问；
7.·
局部集群器：由一组内核共享的片上存储器，能够被一个工作组中的所有线程访问；
8.·
常数存储器：设备片外存储器的一部分，其内容在初始化之后不可更改。一般情况下常数存储器能够自动获得设备内部专用的常数缓存支持，因此访问延迟较小；常数存储器能够被所有线程访问；
9.·
全局存储器：设备片外物理存储器，存储容量较大，但访问延迟较高，能够被所有线程访问；
10.如何优化数据布局，设计高效的数据传输策略，并在主机及各类并行计算设备上合理分配计算任务，充分发挥设备各自的性能优势，是使用opencl求解大规模计算问题时需要考虑的主要问题。
11.线性规划是一类重要的最优化问题，解决的是在线性约束条件下求解线性目标函数极值的问题，广泛应用于社会经济生活的各个方面。线性规划问题的标准模型为：
[0012][0013]
其中，c为目标函数系数向量，a为约束线性方程组的系数矩阵，b为约束线性方程组的右端向量，x为决策变量，m＜n，rank(a)＝m。
[0014]
线性规划问题的基本求解方法为单纯形法，其基本过程如下：
[0015]
s1、令a的前m列为基向量，这些基向量组成基矩阵b，其余列向量组成非基矩阵d；
[0016]
s2、令c
t
的前m个分量为基向量对应的分量，记为向量cb，其余分量组成向量cd；
[0017]
s3、由公式(1)所示线性规划问题标准模型，构建单纯形表：
[0018][0019]
s4、将公式(2)等号右边的矩阵左乘将单纯形表中的增广矩阵[a b]变换为规范型：
[0020][0021]
其中im为m阶单位阵，b-1
为基矩阵b的逆矩阵。如果b-1
不存在，则该线性规划问题无解；
[0022]
s5、将公式(3)等号右边的矩阵左乘将公式(3)等号右边矩阵最后一行中与基列向量对应的元素变为零：
[0023][0024]
公式(4)中等号右边的矩阵称为线性规划问题关于基b的标准单纯形表，其中：
[0025]
c)标准单纯形表的最后一列中，前m个元素b-1
b是关于b的基变量；
[0026]
d)标准单纯形表的最后一行中，称为检验数，记为ri,i＝m+1,
…
,n，为当前基本可行解下目标函数值的相反数；
[0027]
s6、如果所有检验数ri≥0,i＝m+1,
…
,n，则停止计算，当前基本可行解即为最优解，否则进入下一步；
[0028]
s7、分别按照公式(5)和公式(6)计算进基列向量与离基列向量的下标，其中q为进基列向量的下标，p为离基列向量的下标：
[0029]
q＝min{i:ri＜0}
ꢀꢀ
(5)
[0030][0031]
其中y
ij
表示单纯形表中(i,j)处的元素；y
j0
表示单纯形表最后一列中b-1
b的第j个基变量的值；
[0032]
s8、以(p,q)为枢轴元素，按照公式(7)对单纯形表进行枢轴变换：
[0033][0034]
其中yi'j为变换后单纯形表中(i,j)处的元素；转步骤s6。
[0035]
通过对单纯形法的处理步骤进行分析可知，单纯形法可分为两大步骤：一是构建标准单纯形表；二是在标准单纯形表上执行枢轴变换。随着决策变量和约束条件的增多，一方面单纯形表会变得越来越庞大，实际问题中其规模可达数亿甚至数十亿个元素；另一方面，由于枢轴变换需要更新整个单纯形表，且随着决策变量的增多，枢轴变换的次数也不断增加，因此枢轴变换是单纯形法的一个关键性能瓶颈点。
[0036]
本发明提供一种基于异构计算的大规模线性规划问题求解方法，通过分析单纯形法地求解过程，在主机和并行计算设备上合理分配计算任务，能够较快速地求解大规模线性规划问题。


技术实现要素：

[0037]
(一)要解决的技术问题
[0038]
本发明要解决的技术问题是如何提供一种基于异构计算的大规模线性规划问题求解方法，以解决单纯形法求解线性规划问题的瓶颈问题。
[0039]
(二)技术方案
[0040]
为了解决上述技术问题，本发明提出一种基于异构计算的大规模线性规划问题求解方法，该方法包括如下步骤：
[0041]
步骤1：将线性规划标准模型中的相关参数c,a,b输入计算机，同时，还需给定初始基列向量下标集合v；
[0042]
步骤2：主机初始化opencl异构计算环境，其中，主机读取并解析枢轴变换模块；
[0043]
步骤3：主机按照公式(2)构建线性规划问题的单纯形表；
[0044]
步骤4：主机检查基矩阵是否可逆，如果基矩阵可逆，转下一步；否则该线性规划问题无解，算法结束；
[0045]
步骤5：主机按照公式(3)和公式(4)将步骤3中所得的单纯形表转换为标准单纯形表：
[0046]
步骤6：主机向opencl设备申请内存空间；
[0047]
步骤6.1：opencl设备按照标准单纯形表的大小在全局存储器上创建缓冲区；
[0048]
步骤6.2：opencl设备按照标准单纯形表1行的大小在常数存储器上创建缓冲区；
[0049]
步骤7：主机根据枢轴变换模块，向opencl设备申请创建内核(kernel)对象并设置参数，内核对象的参数与枢轴变换模块的参数相同；
[0050]
步骤8：主机将标准单纯形表传输到opencl设备上，具体传输到由步骤6.1所创建的全局存储器缓冲区中；
[0051]
步骤9：主机检查检验数是否全部为非负数；如果检验数全部为非负数，说明当前
基列向量对应的解即最优解，转步骤18；否则转下一步；
[0052]
步骤10：主机按照公式(5)计算进基变量q；
[0053]
步骤11：opencl设备将更新后的单纯形表的第q列传输到主机；
[0054]
步骤12：opencl设备将更新后的单纯形表的最后1列传输到主机；
[0055]
步骤13：主机按照公式(6)计算离基变量p；
[0056]
步骤14：主机更新基列向量下标集合v，令v
p
＝q；
[0057]
步骤15：主机将(p,q)传输到opencl设备；
[0058]
步骤16：opencl设备以(p,q)位置的元素为枢轴元素，对单纯形表进行并行枢轴变换；
[0059]
步骤17：opencl设备将更新后的单纯形表的最后1行传输到主机；转步骤9；
[0060]
步骤18：opencl设备将更新后的单纯形表的最后1列传输到主机；
[0061]
步骤19：主机按照当前基列向量的位置从标准单纯形表的最后一列中获取最优解；
[0062]
步骤20：主机清理opencl计算环境，释放分配的内存和相关对象，计算结束；
[0063]
线性规划问题的标准模型为：
[0064][0065]
其中，c为目标函数系数向量，a为约束线性方程组的系数矩阵，b为约束线性方程组的右端向量，x为决策变量，m＜n，rank(a)＝m。
[0066]
采用单纯形法求解线性规划问题的基本过程如下：
[0067]
s1、令a的前m列为基向量，这些基向量组成基矩阵b，其余列向量组成非基矩阵d；
[0068]
s2、令c
t
的前m个分量为基向量对应的分量，记为向量cb，其余分量组成向量cd；
[0069]
s3、由公式(1)所示线性规划问题标准模型，构建单纯形表：
[0070][0071]
s4、将公式(2)等号右边的矩阵左乘将单纯形表中的增广矩阵[a b]变换为规范型：
[0072][0073]
其中im为m阶单位阵，b-1
为基矩阵b的逆矩阵。如果b-1
不存在，则该线性规划问题无解；
[0074]
s5、将公式(3)等号右边的矩阵左乘将公式(3)等号右边矩阵最后一行中与基列向量对应的元素变为零：
[0075][0076]
公式(4)中等号右边的矩阵称为线性规划问题关于基b的标准单纯形表，其中：
[0077]
e)标准单纯形表的最后一列中，前m个元素b-1
b是关于b的基变量；
[0078]
f)标准单纯形表的最后一行中，称为检验数，记为ri,i＝m+1,
…
,n，为当前基本可行解下目标函数值的相反数；
[0079]
s6、如果所有检验数ri≥0,i＝m+1,
…
,n，则停止计算，当前基本可行解即为最优解，否则进入下一步；
[0080]
s7、分别按照公式(5)和公式(6)计算进基列向量与离基列向量的下标，其中q为进基列向量的下标，p为离基列向量的下标：
[0081]
q＝min{i:ri＜0}
ꢀꢀ
(5)
[0082][0083]
其中y
ij
表示单纯形表中(i,j)处的元素；y
j0
表示单纯形表最后一列中b-1
b的第j个基变量的值；
[0084]
s8、以(p,q)为枢轴元素，按照公式(7)对单纯形表进行枢轴变换：
[0085][0086]
其中yi'j为变换后单纯形表中(i,j)处的元素；转步骤s6。
[0087]
(三)有益效果
[0088]
本发明提出一种基于异构计算的大规模线性规划问题求解方法，本发明的创新之处体现在：一是针对主机和设备的不同特点，对单纯形法整体计算任务进行划分：由主机完成检验数的判断以及进基、离基向量的求解，由计算设备完成枢轴变换。这种设计充分利用了主机善于逻辑判断和计算设备善于并行执行大量简单计算的特点；二是针对上述任务分解，设计了合理的主机与设备间的数据传输策略，在每轮次枢轴变换过程中尽力保持较小的数据传输量(单纯形表的1行+2列)，进一步提高并行计算效率；三是在计算设备上，设计实现枢轴变换的并行计算方法，在每轮次枢轴变换前将单纯形表的第p行复制到常数存储区，充分利用并行计算设备常数存储区的性能优势，同时将整个单纯形表上的枢轴变换分配到计算设备的多个工作项上，进一步提高算法的执行效率。
附图说明
[0089]
图1为opencl平台模型示意图；
[0090]
图2为opencl存储器模型示意图；
[0091]
图3为基于异构计算的大规模线性规划问题求解方法的基本过程示意图；
[0092]
图4为初始化opencl计算环境的基本过程；
[0093]
图5为并行枢轴变换模块处理流程图；
[0094]
图6为主机向opencl设备申请内存空间的过程；
[0095]
图7为主机向opencl设备申请创建内核对象并设置参数的过程；
[0096]
图8为主机将(p,q)传输到opencl设备的过程；
[0097]
图9为opencl设备进行枢轴变换的过程。
具体实施方式
[0098]
本发明提出一种基于异构计算的大规模线性规划问题求解方法。利用主机和并行计算设备的不同特点，将线性规划问题的标准求解算法单纯形法进行并行化。首先对单纯形的计算任务进行分解，为主机和计算设备分配不同的计算任务，在此基础上，在主机与设备间设计了合理的数据传输策略，提高计算性能；最后设计了枢轴变换的并行算法，利用计算设备的并行性克服了单纯形法枢轴变换轮次多、计算量大的缺点。
[0099]
本发明提供一种基于异构计算的大规模线性规划问题求解方法，其基本流程如图3所示，其详细步骤包括：
[0100]
步骤1：将线性规划标准模型中的相关参数c,a,b输入计算机，同时，还需给定初始基列向量下标集合v；
[0101]
步骤2：主机初始化opencl异构计算环境，，其中，主机读取并解析枢轴变换模块，如其步骤如图4所示：
[0102]
步骤2.1：获取当前计算机上的异构计算平台id；
[0103]
步骤2.2：获取当前异构计算平台上的设备id；
[0104]
步骤2.3：使用当前异构计算平台上的设备id创建异构计算环境上下文对象；
[0105]
步骤2.4：使用异构计算环境上下文对象创建命令队列，主机通过命令队列向opencl设备发送命令；
[0106]
步骤2.5：主机读取并解析枢轴变换模块，该模块包含6个输入参数，分别为mat、row_p、n_row、n_col、p、q：
[0107]
1)mat：存储标准单纯形表的数组，位于全局存储器；
[0108]
2)row_p：存储单纯形表第p行的数组，位于常数存储器；
[0109]
3)n_row：单纯形表的行数，正整数
[0110]
4)n_col：单纯形表的列数，正整数
[0111]
5)p：枢轴变换参数
[0112]
6)q：枢轴变换参数
[0113]
枢轴变换模块负责在设备上按照公式(7)对单纯形表进行并行枢轴变换，其详细处理步骤如图5所示。后续步骤16.3中，将调用根据该模块创建的内核对象进行并行枢轴变换：
[0114]
步骤2.5.1：获取当前工作项编号，记为i；
[0115]
步骤2.5.2：从mat中获取单纯形表第i行、第q的元素，记为y
iq
；
[0116]
步骤2.5.3：获取row_p的第q个元素，记为y
pq
；
[0117]
步骤2.5.4：计算r＝y
iq
/y
pq
；
[0118]
步骤2.5.5：令j＝1；
[0119]
步骤2.5.6：检查是否j≥1并且j≤n_col，如果是，转步骤2.5.7，否则模块调用结
束；
[0120]
步骤2.5.7：检查是否i＝p，是则转步骤2.5.8，否则转步骤2.5.9；
[0121]
步骤2.5.8：计算y
′
ij
＝y
ij
/y
pq
，转步骤2.5.11；
[0122]
步骤2.5.9：获取row_p的第j个元素，记为y
pj
；
[0123]
步骤2.5.10：计算y
′
ij
＝y
ij-y
pj
×
r，转步骤2.5.11；
[0124]
步骤2.5.11：j＝j+1；
[0125]
步骤3：主机按照公式(2构建线性规划问题的单纯形表；
[0126]
步骤4：主机检查基矩阵是否可逆，如果基矩阵可逆，转下一步；否则该线性规划问题无解，算法结束；
[0127]
步骤5：主机按照公式(3)和公式(4)将步骤3中所得的单纯形表转换为标准单纯形表：
[0128]
步骤6：主机向opencl设备申请内存空间，如图6所示：
[0129]
步骤6.1：opencl设备按照标准单纯形表的大小在全局存储器上创建缓冲区；
[0130]
步骤6.2：opencl设备按照标准单纯形表1行的大小在常数存储器上创建缓冲区；
[0131]
步骤7：主机根据枢轴变换模块，向opencl设备申请创建内核(kernel)对象并设置参数，内核对象的参数与枢轴变换模块的参数相同，如图7所示：
[0132]
步骤7.1：opencl设备创建内核对象kernel(mat,prow,m,n,p,q)，其中p、q为枢轴变换参数，将在后续步骤15中给定；
[0133]
步骤7.2：opencl设备将内核对象的第1个参数mat设置为步骤6.1中所创建的全局存储器缓冲区；
[0134]
步骤7.3：opencl设备将内核对象的第2个参数prow设置为步骤6.2中所创建的常数存储器缓冲区；
[0135]
步骤7.4：opencl设备将内核对象的第3个参数m设置为标准单纯形表的行数；
[0136]
步骤7.5：opencl设备将内核对象的第4个参数n设置为标准单纯形表的列数；
[0137]
步骤8：主机将标准单纯形表传输到opencl设备上，具体传输到由步骤6.1所创建的全局存储器缓冲区中；
[0138]
步骤9：主机检查检验数是否全部为非负数；如果检验数全部为非负数，说明当前基列向量对应的解即最优解，转步骤18；否则转下一步；
[0139]
步骤10：主机按照公式(5)计算进基变量q；
[0140]
步骤11：opencl设备将更新后的单纯形表的第q列传输到主机；
[0141]
步骤12：opencl设备将更新后的单纯形表的最后1列传输到主机；
[0142]
步骤13：主机按照公式(6)计算离基变量p；
[0143]
步骤14：主机更新基列向量下标集合v，令v
p
＝q；
[0144]
步骤15：主机将(p,q)传输到opencl设备，如图8所示：
[0145]
步骤15.1：opencl设备将内核对象的第5个参数p的值设置为p；
[0146]
步骤15.2：opencl设备将内核对象的第6个参数q的值设置为q；
[0147]
步骤16：opencl设备以(p,q)位置的元素为枢轴元素，对单纯形表进行并行枢轴变换，如图9所示：
[0148]
步骤16.1：opencl设备将当前单纯形表的第p行复制到由步骤6.2所创建的常数存
储器缓冲区中；
[0149]
步骤16.2：结合opencl设备的相关信息，确定执行并行计算的相关参数的值，如工作组个数、工作项个数等；
[0150]
步骤16.3：opencl设备创建相关工作组和工作项，每个工作项按照公式(7)对单纯形表的1行进行枢轴变换；其中所需的第p行的数据从步骤6.2中所创建的常数存储器缓冲区中读取；
[0151]
步骤17：opencl设备将更新后的单纯形表的最后1行传输到主机；转步骤9；
[0152]
步骤18：opencl设备将更新后的单纯形表的最后1列传输到主机；
[0153]
步骤19：主机按照当前基列向量的位置从标准单纯形表的最后一列中获取最优解；
[0154]
步骤20：主机清理opencl计算环境，释放分配的内存和相关对象，计算结束。
[0155]
单纯形法的求解过程中，首先构建标准单纯形表，在此基础上对单纯形表进行多轮次的枢轴变换，每次枢轴变换后检查检验数是否存在负值，如存在负值则须确定离基列向量和进基列向量，然后再次执行枢轴变换，直至检验数全部非负，得到最优解。在此过程中，多轮次的枢轴变换是整个算法的性能瓶颈点。
[0156]
实施例1：
[0157]
通过实际求解一个线性规划问题来说明本发明的具体实施方式。该问题将在cpu+gpu的异构计算环境中进行求解，其中cpu为intel-i711800h，gpu为nvidiageforce rtx 3050ti，计算过程中数据类型采用32位单精度浮点数。
[0158]
需要说明的是，限于篇幅，求解问题的规模并不庞大，实际情况中可应用本发明快速求解更大规模的线性规划问题。
[0159]
问题定义
[0160]
待求解的线性规划问题如下：
[0161][0162]
首先将该问题整理为如公式(1)所示的线性规划问题标准模型，其中：
[0163][0164]
观察到方程ax＝b的增广矩阵的前3列已经是一个基矩阵，因此当前对应的基变量
已经构成一个基本可行解，基列向量下标为v＝[1 2 3]，因此可以从该基本可行解出发，利用本发明提出的并行单纯形法求最优解。
[0165]
求解步骤
[0166]
步骤1：将c,a,b,v作为参数输入计算机，发起计算请求；
[0167]
步骤2：主机初始化opencl异构计算环境，具体步骤包括：
[0168]
步骤2.1：获取当前计算机上的异构计算平台id；
[0169]
步骤2.2：获取当前异构计算平台上的设备id；
[0170]
步骤2.3：使用当前异构计算平台上的设备id创建异构计算环境上下文对象；
[0171]
步骤2.4：使用异构计算环境上下文对象创建命令队列，主机通过命令队列向opencl设备发送命令；
[0172]
步骤2.5：主机读取并解析并行枢轴变换模块；
[0173]
步骤3：主机接收到参数c,a,b,v后，按照公式(2)构造单纯形表：
[0174][0175]
公式(10)所示的单纯形表中，基矩阵
[0176]
步骤4：基矩阵b可逆，转步骤5；
[0177]
步骤5：主机按照公式(3)和(4)将步骤3中得到的单纯形表转换为标准单纯形表(这里标准单纯形表恰巧与原单纯形表相同)：
[0178][0179]
步骤6：主机向opencl设备申请内存空间：
[0180]
步骤6.1：opencl设备在全局存储器上创建4
×8×
4＝128字节大小的缓冲区；
[0181]
步骤6.2：opencl设备在常数存储器上创建8
×
4＝32字节大小的缓冲区；
[0182]
步骤7：主机根据步骤2.5得到的并行枢轴变换模块，向opencl设备申请创建内核对象并为其设置参数：
[0183]
步骤7.1：opencl设备创建内核对象kernel(mat,prow,m,n,p,q)；
[0184]
步骤7.2：opencl设备将内核对象的第1个参数mat设置为步骤6.1中所创建的全局存储器缓冲区；
[0185]
步骤7.3：opencl设备将内核对象的第2个参数prow设置为步骤6.2中所创建的常数存储器缓冲区；
[0186]
步骤7.4：opencl设置将内核对象的第3个参数m设置为4，即标准单纯形表的行数；
[0187]
步骤7.5：opencl设置将内核对象的第4个参数n设置为8，即标准单纯形表的列数；
[0188]
步骤8：主机将标准单纯形表传输到opencl设备上，具体传输到由步骤6.1所创建的全局存储器缓冲区中；
[0189]
第1次变换
[0190]
步骤9：主机检查检验数是否全部非负，检验数因此需要对单纯形表进行枢轴变换：
[0191]
步骤10：主机由公式(5)得到进基变量q＝4；
[0192]
步骤11：opencl设备将单纯形表的第4列a4＝[1/4 1/2 0
ꢀ‑
3/4]
t
传输到主机；
[0193]
步骤12：opencl设备将单纯形表的最后1列[0 0 1 0]
t
传输到主机；
[0194]
步骤13：主机由公式(6)计算离基变量p＝1：
[0195][0196]
步骤14：主机更新基列向量下标集合，得到v＝[4 2 3]；
[0197]
步骤15：主机将p＝1,q＝4传输到opencl设备
[0198]
步骤15.1：opencl设备将内核对象的第5个参数p的值设置为1；
[0199]
步骤15.2：opencl设备将内核对象的第6个参数q的值设置为4；
[0200]
步骤16：opencl设备以y
14
为枢轴元素，对单纯形表进行枢轴变换：
[0201]
步骤16.1：opencl设备将当前单纯形表的第1行复制到步骤6.2所创建的常数存储器缓冲区中；
[0202]
步骤16.2：工作组个数设为4，每个工作组中1个工作项；每个工作项负责更新单纯形表中的1行；
[0203]
步骤16.3：每个工作项按照公式(7)对单纯形表中的1行进行枢轴变换，计算过程中y
1j
的值从步骤6.2中所创建的常数存储器缓冲区中读取；
[0204]
得到更新后的单纯形表，如下所示：
[0205][0206]
步骤17：opencl设备将更新后的单纯形表的最后1行[3 0 0 0
ꢀ‑4ꢀ‑
7/2 33 0]传输到主机；
[0207]
第2次变换
[0208]
步骤18：主机检查检验数是否全部非负，检验数r5＝-4＜0，因此需要对单纯形表进行枢轴变换；
[0209]
步骤19：主机由公式(5)得到进基变量q＝5；
[0210]
步骤20：opencl设备将更新后的标准单纯形表的第5列[-32 4 0
ꢀ‑
4]
t
传输到主机；
[0211]
步骤21：opencl设备将更新后的标准单纯形表的最后1列[0 0 1 0]
t
传输到主机；
[0212]
步骤22：主机由公式(6)计算离基变量p＝2：
[0213][0214]
步骤23：主机更新基列向量集合下标集合，得到v＝[4 5 3]；
[0215]
步骤24：主机将参数p＝2,q＝5传输到opencl设备
[0216]
步骤24.1：opencl设备将内核对象的第5个参数p的值设置为2；
[0217]
步骤24.2：opencl设备将内核对象的第6个参数q的值设置为5；
[0218]
步骤25：opencl设备以y
25
为枢轴元素，对单纯形表进行枢轴变换：
[0219]
步骤25.1：opencl设备将当前单纯形表的第2行复制到步骤6.2所创建的常数存储器缓冲区中；
[0220]
步骤25.2：工作组个数设为4，每个工作组中1个工作项；每个工作项负责更新单纯形表中的1行；
[0221]
步骤25.3：每个工作项按照公式(7)对单纯形表中的1行进行枢轴变换，计算过程中y
2j
的值从步骤6.2中所创建的常数存储器缓冲区中读取；
[0222]
得到更新后的单纯形表，如下所示：
[0223][0224]
步骤26：opencl设备将更新后的单纯形表的最后1行[1 1 0 0 0
ꢀ‑
2 18 0]传输到主机；
[0225]
第3次变换
[0226]
步骤27：主机检查检验数是否全部非负，r6＝-2＜0，因此需要对单纯形表进行枢轴变换；
[0227]
步骤28：主机由公式(5)得到进基变量q＝6；
[0228]
步骤29：opencl设备将单纯形表的第6列a6＝[8 3/8 1
ꢀ‑
2]
t
传输到主机；
[0229]
步骤30：opencl设备将单纯形表的最后1列[0 0 1 0]
t
传输到主机；
[0230]
步骤31：主机由公式(6)计算离基变量p＝1：
[0231][0232]
步骤32：主机更新基列向量下标集合，得到v＝[6 5 3]；
[0233]
步骤33：主机将p＝1,q＝6传输到opencl设备
[0234]
步骤33.1：opencl设备将内核对象的第5个参数p的值设置为1；
[0235]
步骤33.2：opencl设备将内核对象的第6个参数q的值设置为6；
[0236]
步骤34：opencl设备以y
16
为枢轴元素，对单纯形表进行枢轴变换：
[0237]
步骤34.1：opencl设备将当前单纯形表的第1行复制到步骤6.2所创建的常数存储
器缓冲区中；
[0238]
步骤34.2工作组个数设为4，每个工作组中1个工作项；每个工作项负责更新单纯形表中的1行；
[0239]
步骤34.3：每个工作项按照公式(7)对单纯形表中的1行进行枢轴变换，计算过程中y
1j
的值从步骤6.2中所创建的常数存储器缓冲区中读取；
[0240]
得到更新后的单纯形表，如下所示：
[0241][0242]
步骤35：opencl设备将更新后的单纯形表的最后1行[-2 3 0 1/4 0 0
ꢀ‑
3 0]传输到主机；
[0243]
第4次变换
[0244]
步骤36：主机检查检验数是否全部非负，r1＝-2＜0，因此需要对单纯形表进行枢轴变换：
[0245]
步骤37：主机由公式(5)得到进基变量q＝1；
[0246]
步骤38：opencl设备将更新后的单纯形表的第1列a1＝[-3/2 1/16 3/2
ꢀ‑
2]
t
传输到主机；
[0247]
步骤39：opencl设备将更新后的单纯形表的最后1列[0 0 1 0]
t
传输到主机；
[0248]
步骤40：主机由公式(6)计算得到离基变量p＝2：
[0249][0250]
步骤41：主机更新基列向量下标集合，得到v＝[6 1 3]；
[0251]
步骤42：主机将p＝2,q＝1传输到opencl设备
[0252]
步骤42.1：opencl设备将内核对象的第5个参数p的值设置为2；
[0253]
步骤42.2：opencl设备将内核对象的第6个参数q的值设置为1；
[0254]
步骤43：opencl设备以y
21
为枢轴元素，对单纯形表进行枢轴变换：
[0255]
步骤43.1：opencl设备将当前单纯形表的第2行复制到步骤6.2所创建的常数存储器缓冲区中；
[0256]
步骤43.2：工作组个数设为4，每个工作组中1个工作项；每个工作项负责更新单纯形表中的1行；
[0257]
步骤43.3：每个工作项按照公式(7)对单纯形表中的1行进行枢轴变换，计算过程中y
2j
的值从步骤6.2中所创建的常数存储器缓冲区中读取；
[0258]
得到更新后的单纯形表，如下所示：
[0259][0260]
第5次变换
[0261]
步骤44：主机检查检验数是否全部非负，r2＝-1＜0，因此需要对单纯形表进行枢轴变换：
[0262]
步骤45：主机由公式(5)得到进基变量q＝2；
[0263]
步骤46：opencl设备将单纯形表的第2列a2＝[-2
ꢀ‑
2 2
ꢀ‑
1]
t
传输到主机；
[0264]
步骤47：opencl设备将单纯形表的最后1列[0 0 1 0]
t
传输到主机；
[0265]
步骤48：主机由公式(6)计算离基变量p＝3：
[0266][0267]
步骤49：主机更新基列向量下标集合，得到v＝[6 1 2]；
[0268]
步骤50：主机将p＝3,q＝2传输到opencl设备
[0269]
步骤50.1：opencl设备将内核对象的第5个参数p的值设置为3；
[0270]
步骤50.2：opencl设备将内核对象的第6个参数q的值设置为2；
[0271]
步骤51：opencl设备以y
32
＝2为枢轴元素，对单纯形表进行枢轴变换：
[0272]
步骤51.1：opencl设备将当前单纯形表的第3行复制到步骤6.2所创建的常数存储器缓冲区中；
[0273]
步骤51.2：工作组个数设为4，每个工作组中1个工作项；每个工作项负责更新单纯形表中的1行；
[0274]
步骤51.3：每个工作项按照公式(7)对单纯形表中的1行进行枢轴变换，计算过程中y
3j
的值从步骤6.2中所创建的常数存储器缓冲区中读取；
[0275]
得到更新后的单纯形表，如下所示：
[0276][0277]
步骤52：opencl设备将单纯形表的最后1行[0 0 1/2
ꢀ‑
3/4 20 0 6 1/2]传输到主机；
[0278]
第6次变换
[0279]
步骤53：主机检查检验数是否全部非负，其中因此需要对单纯形表进行枢轴变换；
[0280]
步骤54：主机由公式(5)得到进基变量q＝4；
[0281]
步骤55：opencl设备将单纯形表的第4列a4＝[0 1/4 1/2
ꢀ‑
3/4]
t
传输到主机；
[0282]
步骤56：opencl设备将单纯形表的最后1列[1 1 1/2 1/2]
t
传输到主机；
[0283]
步骤57：主机由公式(6)计算离基变量p＝3：
[0284][0285]
步骤58：主机更新基列向量下标集合，得到v＝[6 1 4]；
[0286]
步骤59：主机将p＝3,q＝4传输到opencl设备
[0287]
步骤59.1：opencl设备将内核对象的第5个参数p的值设置为3；
[0288]
步骤59.2：opencl设备将内核对象的第6个参数q的值设置为4；
[0289]
步骤60：opencl设备以y
34
为枢轴元素，对单纯形表进行枢轴变换：
[0290]
步骤60.1：opencl设备将当前单纯形表的第3行复制到步骤6.2所创建的常数存储器缓冲区中；
[0291]
步骤60.2：工作组个数设为4，每个工作组中1个工作项；每个工作项负责更新单纯形表中的1行；
[0292]
步骤60.3：每个工作项按照公式(7)对单纯形表中的1行进行枢轴变换，计算过程中y
3j
的值从步骤6.2中所创建的常数存储器缓冲区中读取；
[0293]
得到更新后的单纯形表，如下所示：
[0294][0295]
步骤61：opencl设备将单纯形表的最后1行[0 3/2 5/4 0 2 0 21/2 5/4]传输到主机；
[0296]
得到最优解
[0297]
步骤62：主机检查检验数全部为非负值，说明当前对应的基本可行解即最优解；
[0298]
步骤63：opencl设备将单纯形表的最后1列[1 3/4 1 5/4]
t
传输到主机；
[0299]
步骤64：基列向量对应下标为v＝[6 1 4]，因此最优解为x
*
＝[3/4 0 0 1 0 1 0]
t
，目标函数值为
[0300]
步骤65：主机清理opencl计算环境，释放分配的内存和相关对象，计算结束。
[0301]
实施例2：
[0302]
问题定义：
[0303]
本发明的方法可以应用于多种线性规划问题所涉及的领域，例如，制定生产计划领域，每种产品都有不同的原材料需求，制定生产计划时，在满足原材料总量约束的条件下，如何安排每种商品的生产数量使得总收入最大。
[0304]
例如，某厂家生产4种不同的产品x1、x2、x3和x4，生产过程需要3种原材料m1、m2和m3，每种产品都有不同的原材料需求。制定生产计划时，厂家需要同时考虑3种原材料的总量，如表1所示，每种生产决策必须满足这些约束条件。假设产品x1、x2、x3和x4的售价分别为6元、4元、7元和5元，在满足原材料总量约束的条件下，厂家应如何安排每种商品的生产数
量x1、x2、x3和x4，以使得总收入最大？
[0305]
表1某厂家生产计划表
[0306][0307]
参考表1数据，得到线性约束方程组：
[0308][0309]
4种产品的总收入可表示为：f(x1,x2,x3,x4)＝6x1+4x2+7x3+5x4[0310]
因此，原问题即求解如下的最大化问题：
[0311][0312]
将f(x1,x2,x3,x4)中x1～x4之前的系数取相反数，将求最大值问题转换为求最小值问题；引入松弛变量x5、x6和x7，将公式(11)所示的最优化问题转换为如公式所示的最优化问题标准模型，相关参数为：
[0313][0314]
观察到方程ax＝b的增广矩阵的后3列已经是一个基矩阵，因此当前对应的基变量已经构成一个基本可行解，基列向量下标为v＝[5 6 7]，因此可以从该基本可行解出发，利用本发明提出的并行单纯形法求最优解。
[0315]
将参数c
t
、a、b和v输入算法(即发明内容步骤1)，其中c
t
为标准线性规划问题目标函数的系数向量；a为约束线性方程组的系数矩阵；b为约束线性方程组的右端向量；v为约束线性方程组的基列向量下标。构造标准单纯形表：
[0316][0317]
按照本发明所述步骤1～步骤20进行求解，求取最优解的过程共涉及4次并行枢轴变换，每次枢轴变换的相关参数及中间结果如表2所示。第4次枢轴变换完成后，检验数已全部非负，因此最优解为x＝[16 0 0 2]
t
，即在当前资源条件约束下，厂家生产16件x1产品以及2件x4产品可达到最大收入106元。
[0318]
表2具体实施方式2中最优化问题的计算过程和中间结果
[0319][0320][0321]
本发明的创新之处体现在：一是针对主机和设备的不同特点，对单纯形法整体计算任务进行划分：由主机完成检验数的判断以及进基、离基向量的求解，由计算设备完成枢轴变换。这种设计充分利用了主机善于逻辑判断和计算设备善于并行执行大量简单计算的特点；二是针对上述任务分解，设计了合理的主机与设备间的数据传输策略，在每轮次枢轴变换过程中尽力保持较小的数据传输量(单纯形表的1行+2列)，进一步提高并行计算效率；三是在计算设备上，设计实现枢轴变换的并行计算方法，在每轮次枢轴变换前将单纯形表的第p行复制到常数存储区，充分利用并行计算设备常数存储区的性能优势，同时将整个单纯形表上的枢轴变换分配到计算设备的多个工作项上，进一步提高算法的执行效率。
[0322]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

技术特征：
1.一种基于异构计算的大规模线性规划问题求解方法，其特征在于，该方法包括如下步骤：步骤1：将线性规划标准模型中的相关参数c,a,b输入计算机，同时，还需给定初始基列向量下标集合v；步骤2：主机初始化opencl异构计算环境，其中，主机读取并解析枢轴变换模块；步骤3：主机按照公式(2)构建线性规划问题的单纯形表；步骤4：主机检查基矩阵是否可逆，如果基矩阵可逆，转下一步；否则该线性规划问题无解，算法结束；步骤5：主机按照公式(3)和公式(4)将步骤3中所得的单纯形表转换为标准单纯形表：步骤6：主机向opencl设备申请内存空间；步骤6.1：opencl设备按照标准单纯形表的大小在全局存储器上创建缓冲区；步骤6.2：opencl设备按照标准单纯形表1行的大小在常数存储器上创建缓冲区；步骤7：主机根据枢轴变换模块，向opencl设备申请创建内核(kernel)对象并设置参数，内核对象的参数与枢轴变换模块的参数相同；步骤8：主机将标准单纯形表传输到opencl设备上，具体传输到由步骤6.1所创建的全局存储器缓冲区中；步骤9：主机检查检验数是否全部为非负数；如果检验数全部为非负数，说明当前基列向量对应的解即最优解，转步骤18；否则转下一步；步骤10：主机按照公式(5)计算进基变量q；步骤11：opencl设备将更新后的单纯形表的第q列传输到主机；步骤12：opencl设备将更新后的单纯形表的最后1列传输到主机；步骤13：主机按照公式(6)计算离基变量p；步骤14：主机更新基列向量下标集合v，令v
p
＝q；步骤15：主机将(p,q)传输到opencl设备；步骤16：opencl设备以(p,q)位置的元素为枢轴元素，对单纯形表进行并行枢轴变换；步骤17：opencl设备将更新后的单纯形表的最后1行传输到主机；转步骤9；步骤18：opencl设备将更新后的单纯形表的最后1列传输到主机；步骤19：主机按照当前基列向量的位置从标准单纯形表的最后一列中获取最优解；步骤20：主机清理opencl计算环境，释放分配的内存和相关对象，计算结束；其中，线性规划问题的标准模型为：其中，c为目标函数系数向量，a为约束线性方程组的系数矩阵，b为约束线性方程组的右端向量，x为决策变量，m＜n，rank(a)＝m；采用单纯形法求解线性规划问题的过程如下：s1、令a的前m列为基向量，这些基向量组成基矩阵b，其余列向量组成非基矩阵d；
s2、令的前m个分量为基向量对应的分量，记为向量c
b
，其余分量组成向量c
d
；s3、由公式(1)所示线性规划问题标准模型，构建单纯形表：s4、将公式(2)等号右边的矩阵左乘将单纯形表中的增广矩阵[a b]变换为规范型：其中i
m
为m阶单位阵，b-1
为基矩阵b的逆矩阵；如果b-1
不存在，则该线性规划问题无解；s5、将公式(3)等号右边的矩阵左乘将公式(3)等号右边矩阵最后一行中与基列向量对应的元素变为零：公式(4)中等号右边的矩阵称为线性规划问题关于基b的标准单纯形表，其中：a)标准单纯形表的最后一列中，前m个元素b-1
b是关于b的基变量；b)标准单纯形表的最后一行中，称为检验数，记为为检验数，记为为当前基本可行解下目标函数值的相反数；s6、如果所有检验数r
i
≥0,i＝m+1,
…
,n，则停止计算，当前基本可行解即为最优解，否则进入下一步；s7、分别按照公式(5)和公式(6)计算进基列向量与离基列向量的下标，其中q为进基列向量的下标，p为离基列向量的下标：q＝min{i:r
i
＜0}
ꢀꢀꢀꢀꢀꢀ
(5)其中y
ij
表示单纯形表中(i,j)处的元素；y
j0
表示单纯形表最后一列中b-1
b的第j个基变量的值；s8、以(p,q)为枢轴元素，按照公式(7)对单纯形表进行枢轴变换：其中y
′
ij
为变换后单纯形表中(i,j)处的元素；转步骤s6。2.如权利要求1所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，所述步骤2具体包括：
步骤2.1：获取当前计算机上的异构计算平台id；步骤2.2：获取当前异构计算平台上的设备id；步骤2.3：使用当前异构计算平台上的设备id创建异构计算环境上下文对象；步骤2.4：使用异构计算环境上下文对象创建命令队列，主机通过命令队列向opencl设备发送命令；步骤2.5：主机读取并解析枢轴变换模块。3.如权利要求2所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，枢轴变换模块包含6个输入参数，分别为mat、row_p、n_row、n_col、p和q，mat：存储标准单纯形表的数组，位于全局存储器；row_p：存储单纯形表第p行的数组，位于常数存储器；n_row：单纯形表的行数，正整数，n_col：单纯形表的列数，正整数，p：枢轴变换参数，q：枢轴变换参数。4.如权利要求3所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，枢轴变换模块负责在设备上按照公式(7)对单纯形表进行并行枢轴变换。5.如权利要求3所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，步骤2.5具体包括如下步骤：步骤2.5.1：获取当前工作项编号，记为i；步骤2.5.2：从mat中获取单纯形表第i行、第q的元素，记为y
iq
；步骤2.5.3：获取row_p的第q个元素，记为y
pq
；步骤2.5.4：计算r＝y
iq
/y
pq
；步骤2.5.5：令j＝1；步骤2.5.6：检查是否j≥1并且j≤n_col，如果是，转步骤2.5.7，否则模块调用结束；步骤2.5.7：检查是否i＝p，是则转步骤2.5.8，否则转步骤2.5.9；步骤2.5.8：计算y
′
ij
＝y
ij
/y
pq
，转步骤2.5.11；步骤2.5.9：获取row_p的第j个元素，记为y
pj
；步骤2.5.10：计算y
′
ij
＝y
ij-y
pj
×
r，转步骤2.5.11；步骤2.5.11：j＝j+1。6.如权利要求5所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，所述步骤7具体包括如下步骤：步骤7.1：opencl设备创建内核对象kernel(mat,prow,m,n,p,q)，其中p、q为枢轴变换参数，将在后续步骤15中给定；步骤7.2：opencl设备将内核对象的第1个参数mat设置为步骤6.1中所创建的全局存储器缓冲区；步骤7.3：opencl设备将内核对象的第2个参数prow设置为步骤6.2中所创建的常数存储器缓冲区；步骤7.4：opencl设备将内核对象的第3个参数m设置为标准单纯形表的行数；步骤7.5：opencl设备将内核对象的第4个参数n设置为标准单纯形表的列数。7.如权利要求6所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，所述步骤15具体包括如下步骤：步骤15.1：opencl设备将内核对象的第5个参数p的值设置为p；
步骤15.2：opencl设备将内核对象的第6个参数q的值设置为q。8.如权利要求7所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，所述步骤16具体包括如下步骤：步骤16.1：opencl设备将当前单纯形表的第p行复制到由步骤6.2所创建的常数存储器缓冲区中；步骤16.2：结合opencl设备的相关信息，确定执行并行计算的相关参数的值；步骤16.3：opencl设备创建相关工作组和工作项，每个工作项按照公式(7)对单纯形表的1行进行枢轴变换；其中所需的第p行的数据从步骤6.2中所创建的常数存储器缓冲区中读取。9.如权利要求8所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，步骤16.2中，相关参数包括：工作组个数和工作项个数。10.如权利要求1所述的基于异构计算的大规模线性规划问题求解方法，其特征在于，该方法用于制定生产计划领域。

技术总结
本发明涉及一种基于异构计算的大规模线性规划问题求解方法，属于高性能计算技术领域。本发明利用主机和并行计算设备的不同特点，将线性规划问题的标准求解算法单纯形法进行并行化。首先对单纯形的计算任务进行分解，为主机和计算设备分配不同的计算任务，在此基础上，在主机与设备间设计了合理的数据传输策略，提高计算性能；最后设计了枢轴变换的并行算法，利用计算设备的并行性克服了单纯形法枢轴变换轮次多、计算量大的缺点。计算量大的缺点。计算量大的缺点。


技术研发人员：李勤勇 濮约刚 孙大东 张明庆 姜有田 马帅 韩琼 吴磊
受保护的技术使用者：北京计算机技术及应用研究所
技术研发日：2023.07.19
技术公布日：2023/10/19