CN111881408A - 基于预处理共轭梯度法的并行数据处理方法及处理系统 - Google Patents

Abstract

本发明属于计算机数据处理技术领域，公开了一种基于预处理共轭梯度法的并行数据处理方法及处理系统，特殊矩阵处理模块用于预定义特殊矩阵M；计算模块，用于计算r₀、z₀、p₀；残差处理模块，用于迭代计算残差。检测模块，用于检查处理后的数据进行正确性及精度验证本发明源代码是通过修改PAC2019决赛代码来进行优化。源程序为PCG算法的并行程序，本发明在这个基础上实现了CAPCG并行算法，特别是对解大型线性方程组中最难以并行的前推回代部分，提出了一种新的算法解决，也就是基于9点模板的预处理共轭梯度法进行算法结构改动，充分利用多核计算资源，时间相比于PCG算法的速度提升了2倍有余。

Description

基于预处理共轭梯度法的并行数据处理方法及处理系统

技术领域

本发明属于计算机数据处理技术领域，尤其涉及一种基于预处理共轭梯度法的并行数据处理方法及处理系统。

背景技术

目前，线性方程组是线性代数的核心内容之一,其解法研究是现代数学中经典且重要的研究课题。线性方程组的直接解法有消元法,克拉默法则,直接三角形法,平方根法,追赶法等。线性方程组的直接解法是最早出现的求解方法,但是由于直接解法的存储量和计算量一般很大,所以迭代法在近年来越来越受到重视。迭代法有Jacobi迭代法,Gauss-Seidel迭代法,JGS迭代法,SOR迭代法,参数迭代法等。而对于求解大型线性方程组以共轭梯度法和变形共轭梯度法为主,其具有存储量少,计算量少等优点。共轭梯度法是最基本的Krylov子空间方法。共轭梯度法的收敛速度与系数矩阵的条件数紧密相关,条件数愈小,收敛性愈好,该算法可以在很少的几步就会获得高精度的近似解。但当系数矩阵的条件数很大时,收敛速度就很慢。于是出现了预处理共轭梯度法(简称PCG法),它是通过适当的预处理方法引入预处理矩阵M,使矩阵的特征值分布更为集中,降低矩阵条件数,以达到提高收敏速度的目的。

20世纪50年代，计算机科学领域开始向并行技术进行探索。2000年以来，分布主存并行技术发展迅速，同时共享主存并行机不断进步，特别是多处理机服务器的推广应用，増强了共享主存并行处理技术的生命力。目前，采用的主流并行编程技术有OpenMP、MPI、CPU/GPU及MapReduce等。OpenMP(OpenMulti-Processing)属于共享内存编程模型的技术。通过在源代码中加入制导性注释，称为编译制导指令的专用#pragma行来指明程序的并发属性。由于它基于编译制导，具有简单、移植性好、可扩展性髙以及支持增量并行化开发等优点，OpenMP己经成为共享存储系统中的并行编程标准。随着多核处理器的大量使用，OpenMP日益体现其重要性。

现有技术方案:

1952年M.R.Hestenes和E.Stiefel从极小化的观点来讨论代数方程组Ax＝b的解,给出了著名的共轭梯度法(Conjugate Gradient,简称CG),若A是N阶对称正定实矩阵,记向量x和y的内积为(x,y),定义x的二次函数

则对于Ax＝b的精确解x^*有

由(1.1)和(1.2)可得

由A的正定性可知Q(x),Q(x*)非负且仅当x＝x*时为零,因而x*是Q(x)的唯一极小点.为了在N维空间中逼近x*,CG法采用了搜索方向P_k,在该方向上选取x_k的增量,使||x_k+1-x*||达到极小.CG法的计算公式为下述的迭代形式。

其中a_k的选择使||x_k+ap_k-x*||_A对所a达到极小,||y||_A定义为

),β_k的选择使P_k+1与P_k为A正交,即(P_k+1,AP_k)＝0.CG法所得到的搜索方向集合{P_k}是A正交系,余向量集合{r_k}是正交系,从而可知,当k>N-1时必有r_k＝0,即在N次迭代以内得到精确解,因而CG法是具有迭代形式的精确解法。它还是完全向量化的算法,只们极少量的非向量运算.当A为l对角的稀阵时,CG法的每次迭代约需(l+5)N次乘法,(l+4)N个存储单元.

现有技术方案存在的不足：

然而,尽管CG法有不少突出的优点,在相当长的时间里它并未得到广泛的应用。主要是由于它的收敛速度慢,而且由于实际计算中舍入误差的积累,{r_k}的正交性逐渐丧失.N阶的问题达代N次后并不能得到精确解,当N>1000时计算量就很大.既然CG法在实用中效率不高,尤其是五十年代至六十年代出现了一些效率较高的算法,如超松弛法(SOR),交替方向隐式法(ADI)和强隐式法(SIP)等等,人们的兴趣转向于其它方法,CG法则逐渐被人冷落,例如,在1971年国际稀阵及其应用会议的文集中竟然没有一篇论文是讨论CG法的,这种情形一直延续到70年代中期,终于发生了被某些数学家称为“突然的、戏剧性的变化”,预处理共轭梯度法出现了。

通过上述分析，现有技术存在的问题及缺陷为：(1)现有技术数据处理中，线性方程组直接解法的存储量和计算量很大。线性方程组求解算法是一种相对比较复杂、对计算资源要求相对较高的算法，无论是对计算时间还是对内存空间都有较高的要求。计算机行业的迅速发展，带来数据量的迅猛增长，且线性方程组求解的过程中常常涉及巨大的数据阵列和大规模的数据运算，这对大型矩阵算法的效率提出了更高要求。普通计算机己经不能承受单次数据处理计算量，通用的串行处理已不能满足日益増长的数据量与对数据的快速实时处理运算的要求。需要利用多核，当系数矩阵的条件数很大时，现有技术并没有通过并行方法，降低矩阵条件数，造成数据收敛速度满。

发明内容

为了解决现有技术存在的问题，本发明提供了一种基于预处理共轭梯度法的并行数据处理方法及处理系统。

本发明是这样实现的，一种基于预处理共轭梯度法的并行数据处理系统，应用于计算机，包括：

特殊矩阵处理模块，用于预定义特殊矩阵M；

计算模块，用于计算r₀、z₀、p₀；

残差处理模块，用于迭代计算残差。

进一步，所述基于预处理共轭梯度法的并行数据处理系统进一步包括：

检测模块，用于检查处理后的数据进行正确性及精度验证。

本发明的另一目的在于提供一种基于预处理共轭梯度法的并行数据处理方法，运行在计算机上，包括：

步骤一，预定义特殊矩阵M；

步骤二，计算r₀＝b-Ax₀，z₀＝M^-1r₀，p₀＝z₀，k＝0；

步骤三，迭代计算残差。

进一步，所述步骤一进一步包括：

M为预处理矩阵，M^-1Ax＝M^-1b；

改写为：

x＝(I-M^-1A)x+M^-1b；

迭代形式为

x_k+1＝(I-M^-1A)x_k+M^-1b；

其中B＝I-M^-1A称为迭代矩阵,B的谱半径小I时迭代法收敛；{u_i}是B的特征值集合,{9,i}是B的相应特征向量集合,如果u_i都是单的，误差向量为e_k＝x*-x_k且有

则ek满足

M^-1A的特征值越接近于1,u_i就越接近于0,于是e_k收敛到零向量的速度越快。

进一步，所述步骤二进一步包括：

对(LL^T)^-1Ax＝(LL^T)^-1b作CG,采用的计算公式为:

X₀任意给定,r₀＝b-Ax₀,p₀＝(LL^T)^-1r₀,

若A是l₁对角,L是l₂对角,则每次迭代需(l₁+2l₂+5)N次乘除法；当A为对称的M阵时，i≠j时a_ij≤0,A非奇且A^-1的元素均为非负，则有：

1)序列{x_k}在N次选代以内收敛到x*；

2)设E_k(T)为T＝(LL^T)^-1A的k次多项式集合,在全部形式为

x_k＝x₀+E_k-1(T)·T(x^*-x₀)；

使||x_k-x^*||_A达到极小；

3)若c为T的条件数,即c＝λ_max(T)/λ_min(T),则有

进一步，所述对(LL^T)^-1Ax＝(LL^T)^-1b作CG中，(LL^T)^-1r_k计算方法包括：

首先解Lw_k＝r_k,采用脚标递增的次序算出w_k的各分量；

再解L^Ty_k＝w_k,采用脚标递减的次序算出y_k的各分量。

进一步，所述步骤三迭代计算残差后，检查正确性及精度。

本发明的另一目的在于提供一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行所述的基于预处理共轭梯度法的并行数据处理方法，包括：

步骤1，预定义特殊矩阵M；

步骤2，计算r₀＝b-Ax₀，z₀＝M^-1r₀，p₀＝z₀，k＝0；

步骤3，迭代计算残差。

步骤4，检查正确性及精度。

本发明的另一目的在于提供一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施所述的基于预处理共轭梯度法的并行数据处理方法所述的方法。

本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述的基于预处理共轭梯度法的并行数据处理方法。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：

本发明源代码是通过修改PAC2019决赛代码来进行优化。源程序为PCG算法的并行程序，本发明在这个基础上实现了CAPCG并行算法，特别是对解大型线性方程组中最难以并行的前推回代部分，提出了一种新的算法解决，也就是基于9点模板的预处理共轭梯度法进行算法结构改动，充分利用多核计算资源，时间相比于PCG算法的速度提升了2倍有余。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的基于预处理共轭梯度法的并行数据处理方法流程图。

图2是本发明实施例提供的则L^T或U的元素模示意图；图中：(a)A的近似LU分解图；(b)无填入的L₀分解图；(c)有三条非零对角线填入的L₃分解图。

图3是本发明实施例提供的特征值分布的比较图。图中：

的特征值；

的特征值；

的特征值。

图4是本发明实施例提供的不同方法的比较图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

现有技术数据处理中，线性方程组直接解法的存储量和计算量很大。线性方程组求解算法是一种相对比较复杂、对计算资源要求相对较高的算法，无论是对计算时间还是对内存空间都有较高的要求。计算机行业的迅速发展，带来数据量的迅猛增长，且线性方程组求解的过程中常常涉及巨大的数据阵列和大规模的数据运算，这对大型矩阵算法的效率提出了更高要求。普通计算机己经不能承受单次数据处理计算量，通用的串行处理已不能满足日益増长的数据量与对数据的快速实时处理运算的要求。需要利用多核，当系数矩阵的条件数很大时，现有技术并没有通过并行方法，降低矩阵条件数，造成数据收敛速度满。

针对现有技术存在的问题，本发明提供了一种基于预处理共轭梯度法的并行数据处理方法及处理系统，下面结合附图对本发明作详细的描述。

如图1所示，本发明实施例提供的基于预处理共轭梯度法的并行数据处理方法包括：

S101，预定义一个特殊矩阵M。

S102，计算r₀＝b-Ax₀，z₀＝M^-1r₀，p₀＝z₀，k＝0。本步骤为后续循环求解作初始化准备。

S103，迭代计算残差。

S104，检查正确性及精度。

具体的，包括：

在CG法的计算中矩阵A是未经过加工的,一般说来其条件数都比较大。如果对A进行某种形式的加工改造,可以使它具有更适宜于迭代的特点,使迭代收敛速度大大提高。实际上,当L非奇时,Ax＝b等价于(LL^T)^-1Ax＝(LL^T)^-1b.若LL^T是A的一个近似LU分解,则(LL^T)^{- 1}A应在某种模的意义下近似于单位阵1.若用它来作CG法,可以加速迭代的收敛。如求解下述椭圆型偏微分方程:

采用

的差分网格和五点差分近似时,在未知数为自然排序的情形下有

A＝tridiag(-I,D_i,-I)；

显然,A具有弱对角优势,它的精确LU分解的对角线元素的绝对值分布具有单调性。若A＝LU,L为下三角阵,U为上三角阵,则L^T或U的元素模按图2所示的方向递降。这种单调性尚无理论证明,但已是众所周知的事实,从这种单调性出发,忽略那些绝对值较小的对角线,图2(a)可以得到A的近似LU分解,由此出发,引申出图2(b)无填入的L₀和图2(c)有三条非零对角线填入的L₃。

L₀与A的下三角部分有同样的稀疏结构,而L₃则比A的下三角部分多了三条非零对角线(以虚线表示)。用(L₀L₀ ^T)^-1作预处理,使(L₀L₀ ^T)^-1A的特征值比A的特征值更接近于1,而(L₃L₃ ^T)^-1A又要比(L₀L₀ ^T)^-1A更好。

至于L阵的各非零对角线的元素,则按下述方法来计算。以L₃为例,L₃L₃ ^T是15对角的。除其中有4条对角线是不可避免的误差以外,有5条对角线与A的相应非零对角线相同,另外还有6条对角线可令其为零,这样共有11个条件,可以决定L₃和L₃ ^T的11条非零对角线。由于A对称,简化为6个条件决定L₃的6条非零对角线。

对于N＝36的一个算例,给出其特征值分布,如图3所示。

步骤S101中，特征值分布的这种变化，对于迭代法的收敛速度影响极大。以下述的一阶线性定常迭代法为例。设M为预处理矩阵，把

M^-1Ax＝M^-1b (2.2)

改写成如下的形式：

x＝(I-M^-1A)x+M^-1b (2.3)

其迭代形式为

x_k+1＝(I-M^-1A)x_k+M^-1b (2.4)

其中B＝I-M^-1A称为迭代矩阵,当B的谱半径小I时迭代法(2.4)收敛。设{u_i}是B的特征值集合,{9,i}是B的相应特征向量集合,假定u_i都是单的。令误差向量为e_k＝x*-x_k且有

则ek满足

由(2.6)可知,M^-1A的特征值越接近于1,u_i就越接近于0,于是e_k收敛到零向量的速度就越快。

步骤S102中，对(LL^T)^-1Ax＝(LL^T)^-1b作CG,它所采用的计算公式为:

X₀任意给定,r₀＝b-Ax₀,p₀＝(LL^T)^-1r₀,

若A是l₁对角的,L是l₂对角的,则(2.7)每次迭代约需(l₁+2l₂+5)N次乘除法。当A为对称的M阵(i≠j时a_ij≤0,A非奇且A^-1的元素均为非负)时：

1)序列{x_k}在N次选代以内收敛到x*。

2)设E_k(T)为T＝(LL^T)^-1A的k次多项式集合,则在全部形式为

x_k＝x₀+E_k-1(T)·T(x^*-x₀)

的算法中,(2.7)使||x_k-x^*||_A达到极小。

3)若c为T的条件数,即c＝λ_max(T)/λ_min(T),则有

当x₀为(0,2)内的随机向量时,取n＝30,N＝900,计算结果如图4所示。

由图4可以看出,采用

作预处理,使计算效率远远高于超松弛法和强隐式法.采用

写作预处理,又使计算效率进一步提高。若不作预处理,则计算效率低于强隐式法。

在(2.7)的计算中有(LL^T)^-1r_k出现。它的计算可以分成两步。先解Lw_k＝r_k,由于L是下三角阵,采用脚标递增的次序即可算出w_k的各分量。再解L^Ty_k＝w_k,采用脚标递减的次序即可算出y_k的各分量.这样就避免了直接计算(LL^T)^-1。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到本发明可借助软件加必需的硬件平台的方式来实现，当然也可以全部通过硬件来实施。基于这样的理解，本发明的技术方案对背景技术做出贡献的全部或者部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。

以上所述仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种基于预处理共轭梯度法的并行数据处理系统，搭载在计算机上，其特征在于，所述基于预处理共轭梯度法的并行数据处理系统包括：

特殊矩阵处理模块，用于预定义特殊矩阵M；

计算模块，用于计算r₀、z₀、p₀；

残差处理模块，用于迭代计算残差。

2.如权利要求1所述的基于预处理共轭梯度法的并行数据处理系统，其特征在于，所述基于预处理共轭梯度法的并行数据处理系统进一步包括：

检测模块，用于检查处理后的数据进行正确性及精度验证。

3.一种基于预处理共轭梯度法的并行数据处理方法，运行在计算机上，其特征在于，所述基于预处理共轭梯度法的并行数据处理方法包括：

步骤一，预定义特殊矩阵M；

步骤二，计算r₀＝b-Ax₀，z₀＝M^-1r₀，p₀＝z₀，k＝0；

步骤三，迭代计算残差。

4.如权利要求3所述的于预处理共轭梯度法的并行数据处理方法，其特征在于，所述步骤一进一步包括：

M为预处理矩阵，M^-1Ax＝M^-1b；

改写为：

x＝(I-M^-1A)x+M^-1b；

迭代形式为

x_k+1＝(I-M^-1A)x_k+M^-1b；

其中B＝I-M^-1A称为迭代矩阵,B的谱半径小I时迭代法收敛；{u_i}是B的特征值集合,{9,i}是B的相应特征向量集合,如果u_i都是单的，误差向量为e_k＝x*-x_k且有

则ek满足

M^-1A的特征值越接近于1,u_i就越接近于0,于是e_k收敛到零向量的速度越快。

5.如权利要求3所述的于预处理共轭梯度法的并行数据处理方法，其特征在于，所述步骤二进一步包括：

对(LL^T)^-1Ax＝(LL^T)^-1b作CG,采用的计算公式为:

X₀任意给定,r₀＝b-Ax₀,p₀＝(LL^T)^-1r₀,

若A是l₁对角,L是l₂对角,则每次迭代需(l₁+2l₂+5)N次乘除法；当A为对称的M阵时，i≠j时a_ij≤0,A非奇且A^-1的元素均为非负，则有：

1)序列{x_k}在N次选代以内收敛到x*；

2)设E_k(T)为T＝(LL^T)^-1A的k次多项式集合,在全部形式为

x_k＝x₀+E_k-1(T)·T(x^*-x₀)；

使||x_k-x^*||_A达到极小；

3)若c为T的条件数,即c＝λ_max(T)/λ_min(T),则有

6.如权利要求3所述的于预处理共轭梯度法的并行数据处理方法，其特征在于，所述对(LL^T)^-1Ax＝(LL^T)^-1b作CG中，(LL^T)^-1r_k计算方法包括：

首先解Lw_k＝r_k,采用脚标递增的次序算出w_k的各分量；

再解L^Ty_k＝w_k,采用脚标递减的次序算出y_k的各分量。

7.如权利要求3所述的基于预处理共轭梯度法的并行数据处理方法，其特征在于，所述步骤三迭代计算残差后，检查正确性及精度。

8.一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行权利要求3～7任意一项所述的基于预处理共轭梯度法的并行数据处理方法，包括：

步骤1，预定义特殊矩阵M；

步骤2，计算r₀＝b-Ax₀，z₀＝M^-1r₀，p₀＝z₀，k＝0；

步骤3，迭代计算残差。

步骤4，检查正确性及精度。

9.一种计算机设备，其特征在于，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行权利要求3～7任意一项所述的基于预处理共轭梯度法的并行数据处理方法。

10.一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施如权利要求3～7任意一项所述的基于预处理共轭梯度法的并行数据处理方法。

Images