专利名称：基于线程构造块的超大规模矩阵多核并行克劳特分解方法
技术领域：
本发明涉及一种超大规模矩阵分解方法，尤其是涉及基于线程构造块的超大规模矩阵多核并行克劳特(Crout)分解方法。
背景技术：
Intel 刚刚推广的线程构建模块(Threading Building Blocks, TBB)是基于 C++ 的多线程并行编程模型，用于支持多核处理器的并行计算，具有成熟的数据结构，支持可扩展的线程嵌套并行，支持可扩展内存分配、支持多种系统平台以及多种编译器支持。TBB的编程模式是使用模板作为并行迭代模型。这使得程序员很方便的处理同步、负载平衡、缓存优化等问题，而能轻松地实现自动调度的并行程序，充分利用CPU的多核运算能力及其他系统资源。TBB提升了程序的伸缩性及可扩展性，并且完全支持嵌套的并行编程。程序员可以创建自己的并行组件，以构建大型的并行程序。同时由于TBB是芯片制造商htel开发的，对于多核平台有着更好的支持。TBB是一个开源的C++模板库，它可将线程抽象成任务，然后创建可靠、可移植且可扩展的并行应用程序。使用TBB编写基于任务的并行应用程序，有助于提高跨多核平台上运行的可扩展软件的开发效率。通本地线程和线程封装器等其它线程化方法相比，TBB是执行并行应用程序最高效的方式，它能够充分利用多核平台的并行性能。TBB有以下三个显著的优点1)TBB对线程进行抽象，将其提高到任务的高度，简化了并行应用程序的开发工作；2)基于TBB开发的应用程序的性能可以随着处理器内核数量的增加而自动提升；3)TBB 能够减少死锁和资源竞争等常见线程错误，提供了一个跨平台、可扩展的并行解决方案。目前，在理论研究和实际应用中，矩阵求逆都有广泛的应用，如求解线性方程组、 动态规划、计算机图形处理、控制工程等等。目前，矩阵求逆的方法主要有以下几种(1)伴随矩阵法；(2)初等变换法；(3)高斯-约当消去法；(4)三角分解法。下面分别加以介绍。(1)伴随矩阵法。利用矩阵的伴随矩阵及矩阵的行列式来求得矩阵的逆矩阵。此方法一般适用于维数较小的矩阵(2)初等变换法。原矩阵和单位矩阵同时进行初等行(或列)变换，当原矩阵变成单位矩阵的时候，单位矩阵就变成了原矩阵的逆矩阵。(3)高斯-约当消去法又称高斯消去法、高斯消元法，实际上就是我们俗称的加减消元法。数学上，高斯-约当消去法，由高斯和约当得名，它是线性代数中的一个算法，用于求解线性方程组的解，求解矩阵的秩，以及求解可逆方矩阵的逆。当用于一个矩阵时，高斯消去产生“行消去梯形形式”。在消元的过程中，可能出现主元为零的情况，这时消元法无法进行；即使主元不为零但很小时，若用其作除数，会导致其它元素的数量级严重增长和舍入误差的扩散，最后也使得计算解不可靠。为使高斯-约当消去法具有较好的数值稳定性，可以使用完全主元素消去法，但完全主元素消去法在选主元时要花费较多机器时间，为此又考虑使用列主元消去法，但由于仍须选主元，计算量也较大。对于良态问题，高斯-约当消去法也可能给出很坏的结果，这说明高斯-约当消去法的算法很不稳定。事实上，一般的矩阵都是病态矩阵， 采用高斯-约当消去法不能得到满意的结果。(4)三角分解法三角分解法是将方阵分解成一个上三角矩阵和一个下三角矩阵，该方法又被称为 LU分解法。该分解方法的用途主要在简化大矩阵的行列式值的计算，矩阵求逆运算和求解联立方程组。需要注意的是，这种分解方法所得到的上下三角形矩阵不是唯一的，还可找到若干对不同的上下三角矩阵对，它们的乘积也会得到原矩阵。矩阵可以进行LU分解是有条件的，它要求其为方阵且它的所有顺序主子式均不等于零。LU分解的主要优点在于它能用公式表示费时的消去步骤，只对系数矩阵进行操作。LU分解的一个动机是它为矩阵求逆提供了一个有效的方式，也为评估方程组的状态提供了一个方法。在上述的矩阵求逆方法中，伴随矩阵法及初等变换法，易于理解，适合手动求解， 而对于大规模矩阵的求解，这两种方法不太适合；高斯-约当消去法及LU分解法便于程序的实现，在大多数条件下，能够适应大规模矩阵的求解。但是，对于大规模或超大规模矩阵， 高斯-约当消去法及传统的串行LU分解法的计算量往往过大，对计算机的性能要求较高， 算法的实时性不高，对于某些时间有严格要求的应用不太适合。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种能够大幅度提高矩阵LU分解的运行效率、且能跨平台、可扩展、应用广泛的基于线程构造块的超大规模矩阵多核并行克劳特分解方法。本发明的目的可以通过以下技术方案来实现基于线程构造块的超大规模矩阵多核并行克劳特分解方法，其特征在于，该分解方法包括以下步骤1)TBB并行计算平台的安装与环境设置；2)分析并抽取传统Crout方法中可以并行计算的部分，利用TBB并行模板类改写为符合TBB需要的规范类；幻设置问题初始值；4)进入Crout方法的行循环；5)调用parallelreduce并行模块类，将每一行的最大主元保存在临时向量中，同时计算并保存比例因子；6)进入Crout方法的列循环；7)确定新的主元及比例因子，修改TINY，每一行除以主元，判断列循环是否结束，如果判断为否则转到步骤6)，判断行循环是否结束，如果判断为否则转到步骤4) ；8)运行结束，完成超大规模矩阵分解，得到因部分主元法而改变了的行列排列次序。所述的步骤幻中的问题初始值包括矩阵的规模、矩阵的行数、列数，输入矩阵。与现有技术相比，本发明具有效率高、跨平台、可扩展、应用广泛等优点，本发明的并行Crout方法能够大幅度提高矩阵LU分解的运行效率；在此基础上，通过本发明能够大幅度提高矩阵求逆的效率，进而验证了本发明的有效性与高效性。并行Crout方法的适用条件为如果矩阵的所有子式都是非零的，即可使用本发明来进行LU分解，进而可以求得该矩阵的逆或线性方程组的解等。当矩阵规模达到百万级别的时候，本发明比串行Crout 方法约提高了 13.7%，并比高斯-约当消去法约提高了 77%。故本发明具有非常大的实用价值和广泛的应用前景。
图1为本发明的满足TBB的求最大主元类结构图；图2为本发明的并行Crout方法流程图；图3为利用本发明的并行Crout方法求矩阵逆阵的流程图；图4为利用本发明求解矩阵逆阵的结果示意图。
具体实施例方式下面结合附图和具体实施例对本发明进行详细说明。实施例TBB定义了任务的概念，在初始化TBB任务调度时，由任务调度器对象task_ scheduler_init实现多任务的分配和并行计算，支持多线程的划分。在调用并行计算的模板类时，由模板类参数指定循环处理的数值范围以及任务粒度参数。任务粒度参数决定了任务划分的粒度，如果粒度太大，不能充分提高运行效率；如果粒度太小，过度的并行化任务分配造成的开销反而降低了运行效率。在无法获得合适任务粒度情况下可以使用TBB提供的自动分配函数aut0_partiti0ner ()帮助用户设置合适的任务粒度参数。以parallelreduce设计的求解最大主元素模板类为例，详细说明TBB设计并行最大主元素算法的实现过程，以供Crout方法调用。首先，基于parallelreduce模板将此过程封装成类MaxElement，其作用是求解最大主元素。类MaxElement必须包括的接口为 operator和构造函数，具体形式如下class MaxElement {......void operator () (const blocked_range<size_t>&r){void operator () (const blocked_range<size_t>&r){for(size_t i = r. begin () ；i ！ = r. end() ；++i) {......}}}Il分出一个支线，如果要访问X，应该保证原子操作MaxE 1 ement (MaxE 1 ement&x, tbb: split)my_a(χ. my_a), value_of_max(FLT_MIN), index_of_max(_1){......}H合并支线void join(const MaxElement&y){......}}H构造函数MaxElement(DOUBLE^a)my_a(a)，value_of_max(FLT_MIN)，index_of_max(-1){}}；其中operator接口是并行处理的主要部分，主要功能是进行并行循环优化，将循环的参数修改成TBB定义的blockecLrange模板类，能够支持循环体内任务的并行划分。 在operator接口中，并行地执行各小块中的比较大小操作，并返回小块中的最大值。当TBB 决定分出一个 blocked_range 时，会调用 MaxElement (MaxElement&x，tbb sp 1 it)，这个tbb: :split只是一个点位符，用于和拷贝构造函数相区别。当TBB合并blockecLrange 时，会调用voidjoin (MaxElement&y)，本例中合并就意味着把多个结果进行运算，将最终的结果保存并返回。构造函数主要实现参数的初始化以及从整个任务空间中分离并构建子任务。在完成并行模板类的编写之后，初始化TBB任务调度器，调用上述编写的并行模板，返回计算结果，最后结束TBB任务调度。利用TBB设计并行Crout方法的基本步骤如下步骤一 TBB并行计算平台的安装与环境设置；步骤二分析并抽取传统Crout方法中可以并行计算的部分，利用TBB并行模板类改写为符合TBB需要的规范类；步骤三设置问题初始值，即矩阵的规模、矩阵的行数、列数，输入矩阵等；步骤四进入Crout方法的行循环；步骤五调用parallelreduce并行模块类，将每一行的最大主元保存在临时向量中，同时计算并保存比例因子，因此采用TBB进行并行设计对提高方法效率至关重要；步骤六进入Crout方法的列循环；步骤七确定新的主元及比例因子，修改TINY，每一行除以主元。若列循环未结束则转到步骤六，判断行循环是否结束，如果判断为否则转到步骤4)；否则转到步骤四；步骤八结束，得到因部分主元法而改变了的行列排列次序。为了验证并行Crout方法的正确性、有效性及评价本发明的实用性，下面通过本发明求解矩阵的逆来验证与评价，其基本步骤如下步骤一 TBB并行计算平台的安装与环境设置；步骤二 读入文件，其中包括输入矩阵及相关参数设置；步骤三设置问题初始值，即矩阵的行数、列数，输入矩阵；步骤四利用本发明的并行Crout方法对输入矩阵进行LU分解；步骤五进入矩阵的行循环；步骤六进入矩阵的列循环；步骤七通过回代函数将矩阵按列求逆；步骤八增加列计数，若列循环未结束则转到步骤六；增加行计数，若行循环未结束则转到步骤五；步骤九结束，得到输入矩阵的逆矩阵，写入到文件。如图1、2所示，本发明对传统Crout分解方法进行分析，把其中的关键步骤进行并行化处理，并利用Intel公司的线程构建块来实现并行化。本发明充分利用了线程构建块的灵活性、跨平台性及稳定性，实现基于线程构建块的并行Crout分解方法，并将该方法应用于矩阵求逆，经过大量仿真实验显示，本发明能够大幅度提高传统Crout分解方法的效率，进而大幅度提高使用了本发明的矩阵求逆方法的效率。本发明由一个符合线程构建块规范的求最大主元类、最大主元调用类和利用了最大主元调用类的Crout分解方法组成，符合线程构建块规范的求最大主元类提供了具体的并行处理的部分、任务的线程分支划分及合并等功能；最大主元类调用类提供了对求最大主元类的调用封装，且提供了对其进行线程划分的粒度控制参数；并行化的Crout分解方法利用了线程并行化的优点，充分利用了计算机的资源，提高了传统Crout分解方法的效率。本发明可以运行于多个系统平台，如Windows，Linux, Unix等系统。图1为本发明求最大主元类的结构图，以下对图中的各部分的详细描述在部分101中，实现operator接口，其是并行处理的主要部分，主要功能是进行并行循环优化，将循环的参数修改成TBB定义的blockecLrange模板类，能够支持循环体内任务的并行划分；在部分102中，实现带tbb:: split的分支创建函数接口，当TBB决定分出一个 blocked_range 时，会调用 MaxElement (MaxElement&x, tbb: split),这个 tbb: split 只是一个点位符，用于和拷贝构造函数相区别；在部分103中，实现join合并函数接口，当TBB合并blocked_range时，会调用 voidjoin (MaxElement&y)，本例中合并就意味着把多个结果进行运算，将最终的结果保存并返回；在部分104中，实现构造函数接口，构造函数主要实现参数的初始化以及从整个任务空间中分离并构建子任务。图2为本发明的并行Crout方法流程图，下面对图中的各步骤进行详细描述在步骤201中，读入需要分解的矩阵，初始化行值i，列值j，内含比例因子向量 w[i]，极小数TINY，临时变量等。然后执行步骤202 ；在步骤202中，按行循环，递增行，求内含比例因子。然后执行步骤203 ；在步骤203中，调用并行的求最大主元类算法求解各行的主元素和比例因子。然后执行步骤204 ；在步骤204中，按列循环，递增列。然后执行步骤205 ；在步骤205中，计算新的候选主元素及交换比例因子，替换原有的候选主元素及交换比例因子。然后执行步骤206;在步骤206中，判断矩阵是否为奇异矩阵，如果是则赋TINY为零，否则TINY保持不变。然后执行步骤207;在步骤207中，矩阵中的每一行除以主元素。如果按列循环未结束，则执行步骤 204 ；如果按列循环结束但按行循环未结束，则执行步骤202 ；如果按列循环与按行循环都已结束，则执行步骤208;在步骤208中，结束函数输出因部分主元法而改变了的行排列次序，可以根据此排列次序得到矩阵的LU分解矩阵。图3为利用并行Crout方法求矩阵的逆阵的流程图，下面对图中的各步骤进行详细描述在步骤301中，通过对话框选择需要的输入文件，选择的矩阵文件为TXT格式，数据结构格式为第一行为矩阵的行数(row)；第二行为矩阵的列数(col)，在本应用中行数等于列数(row = col = η)，即这里的矩阵为η维方阵；第三行开始为矩阵的值。然后执行步骤302 ；在步骤302中，根据读入的矩阵维数读取需要的数值并初始化相关变量及数据结构。然后执行步骤303;在步骤303中，利用本发明并行Crout方法对读入的矩阵进行分解。然后执行步
7骤 304 ；在步骤304中，控制矩阵的行循环i，然后执行步骤305 ；在步骤305中，控制矩阵的列循环j，然后执行步骤306 ；在步骤306中，通过回代函数对计算的矩阵进行处理。若i < η且j < n，则返回步骤305 ；若i < η且j > = n，则返回步骤304 ；若i > = η且j > = n，执行步骤307 ；在步骤307中，结束矩阵求逆计算，输出矩阵的逆及相关的计算参数方法的名称、运算时间、方法线程数，方法子矩阵大小等。执行步骤308 ；在步骤308中，将结果写入到TXT文件中，供使用者进一步分析。图4中详细比较了并行crout方法、串行crout方法、高斯-约当消去法在求解矩阵时的效率。由参考线可以看出，当矩阵维数大于等于256时，并行crout方法的运行时间优于串行crout方法，更明显优于高斯-约当消去法的运行时间；并且随着矩阵维数的增加，这种趋势越明显。当矩阵规模达到百万级别的时候，本发明比串行Crout方法约提高了 13.7%，并比高斯-约当消去法约提高了 77%。即对于大规模或超大规模矩阵而言，并行 crout方法明显优于高斯-约当消去法，也优于串行crout方法。因此，本发明具有很高的实用价值与广泛的应用前景。
权利要求
1.基于线程构造块的超大规模矩阵多核并行克劳特分解方法，其特征在于，该分解方法包括以下步骤1)线程构造块(TBB)并行计算平台的安装与环境设置；2)分析并抽取传统克劳特(Crout)方法中可以并行计算的部分，利用TBB并行模板类改写为符合TBB需要的规范类；3)设置问题初始值；4)进入Crout方法的行循环；5)调用parallelreduce并行模块类，将每一行的最大主元保存在临时向量中，同时计算并保存比例因子；6)进入Crout方法的列循环；7)确定新的主元及比例因子，修改TINY，每一行除以主元，判断列循环是否结束，如果判断为否则转到步骤6)，判断行循环是否结束，如果判断为否则转到步骤4)；8)运行结束，完成超大规模矩阵分解，得到因部分主元法而改变了的行列排列次序。
2.根据权利要求1所述的基于线程构造块的超大规模矩阵多核并行克劳特分解方法， 其特征在于，所述的步骤幻中的问题初始值包括矩阵的规模、矩阵的行数、列数，输入矩
全文摘要
本发明涉及基于线程构造块的超大规模矩阵多核并行克劳特分解方法，该分解方法包括以下步骤1)将传统Crout方法中可以并行计算的部分改写为符合TBB需要的规范类；2)设置问题初始值；3)进入Crout方法的行循环；4)调用parallel_reduce并行模块类，将每一行的最大主元保存在临时向量中，同时计算并保存比例因子；5)进入Crout方法的列循环；6)确定新的主元及比例因子，修改TINY，每一行除以主元；8)完成超大规模矩阵分解，得到因部分主元法而改变了的行列排列次序。与现有技术相比，本发明具有能够大幅度提高矩阵LU分解的运行效率、且能跨平台、可扩展、应用广泛等优点。
文档编号G06F9/44GK102486727SQ20101057185
公开日2012年6月6日 申请日期2010年12月3日 优先权日2010年12月3日
发明者孙剑, 张丽岩, 李克平, 马健 申请人:同济大学