基于分治法求解对称三对角矩阵特征问题的混合并行实现
求实对称三对角矩阵的特征值和特征向量
求实对称三对角矩阵的特征值和特征向量要求求解一个实对称三对角矩阵的特征值和特征向量。
在介绍如何求解之前，首先我们来了解一下实对称三对角矩阵的定义。
实对称三对角矩阵是指矩阵的非零元素主对角线上的元素为a，副对角线上的元素为b，而其他元素均为0。
可以表示为如下形式：[a1b100...0][b1a2b20...0][0b2a3b3...0][00b3a4...0][..................][ 0 0 0 ... bn-1 an ]下面我们将介绍如何求解实对称三对角矩阵的特征值和特征向量。
求解实对称三对角矩阵的特征值和特征向量有多种方法，其中一种常用的方法是通过迭代法，特别是Householder迭代法。
下面我们将介绍这种方法的主要步骤。
1. 首先，将实对称三对角矩阵转化为对称上Hessenberg矩阵。
对称上Hessenberg矩阵是一个具有类似三对角矩阵结构的对称矩阵。
2. 在转化得到的对称上Hessenberg矩阵上应用QR迭代，不断迭代直到矩阵的对角线元素基本上收敛于特征值。
3. 在每次QR迭代中，我们通过施密特正交化方法（Gram-Schmidt orthogonalization）来构建Q矩阵，然后计算出新的矩阵R，并将其与Q相乘，得到下一次迭代的矩阵。
4.在QR迭代的最后一步，我们得到了一个上三角矩阵，其对角线上的元素即为所求的特征值。
5. 然后，我们可以通过反复应用幂迭代法（power iteration method）来求解对应于这些特征值的特征向量。
幂迭代法是一种求解线性代数特征向量的数值方法。
通过上述方法，我们可以求解实对称三对角矩阵的特征值和特征向量。
这种方法具有较高的数值稳定性和计算效率，因此在实际求解中被广泛采用。
需要注意的是，在特征值和特征向量的计算过程中，可能会出现一些特殊情况。
比如矩阵中的主对角线元素不是严格递增或递减的时候，对于这种情况，我们需要进行一些额外的处理。
一类三对角矩阵特征对的分段快速算法
一类三对角矩阵特征对的分段快速算法唐达【摘要】三对角矩阵特征对的计算复杂性一般为O(n2)(n为矩阵的阶).利用一类三对角矩阵特征对的局限性质,采用分段快速算法,其计算复杂性仅为O(n).该算法适用于特征对具局限性的一类大型非对称三对角矩阵,且具有较高的精度;适合于并行计算.最后给出了数值算例.【期刊名称】《上海电机学院学报》【年(卷),期】2014(017)002【总页数】5页(P120-124)【关键词】三对角矩阵;特征对;特征值;特征向量;分段;t向量【作者】唐达【作者单位】上海电机学院数理教学部,上海201306【正文语种】中文【中图分类】O241.6在许多科学与工程计算中，常需要计算矩阵的特征值。
而对于对称矩阵，一般是将其约化为三对角矩阵后再求其特征值；另外，三对角矩阵的特征问题也常作为原始问题出现。
因此，三对角矩阵的特征问题长期来为许多学者所关注，如美籍华人数学家顾明关于三对角矩阵分－治算法的研究成果［1］在SIAM第六届应用线性代数会议上获最佳论文。
当今，计算三对角矩阵特征问题的方法很多，有QL或 QR 算法［2］372－373［3－4］、二（多）分法［2］391－397［5－9］、分－治算法［2］391－397［1，7，10－11］、同伦法［7，12］以及其他的一些迭代算法［7，13］。
其中，二（多）分法、分－治算法、同伦法等都能用于并行计算。
一般来讲，计算一个n阶三对角矩阵的n个特征对，其计算量总不小于O（n2）；而本文利用一类三对角矩阵特征对的局限性质，来计算大型非对称矩阵，其计算量仅为O（n）。
本文的算法也非常适合并行计算，其并行效率较高。
1 三对角矩阵t向量的衰减性质设A 为n 阶实三对角矩阵，t＝（t1，t2，…为n维实向量。
若式（1）中右端向量仅第n个分量w 不为零，则称t为A 所对应的t向量[14-15]。
可以证明，对角占优三对角矩阵（i＝2，3，…，n）。
也就是说，t向量之分量的绝对值是随i的减小而衰减的。
对称三对角矩阵特征值的二分法
数学软件实验任务书课程名称数学软件实验班级数0901实验课题对称三对角矩阵特征值的二分法实验目的熟悉对称三对角矩阵特征值的二分法运用Matlab/C/C++/Java/Maple/Mathematica等其中实验要求一种语言完成对称三对角矩阵特征值的二分法实验内容成绩教师实验1对称三对角矩阵特征值的二分法1 实验原理对于对称三对角矩阵：111222111n n n n n c b b c b C b c b b c ----⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ 设0i b ≠，1,2,,1i n =- ，记特征矩阵I C λ-的左上角的k 阶子式为()k p λ，设0()1p λ=，利用行列式的展开式，可得()k p λ的递推公式：0112112()1() 1,2,,()()()k k k k k p p a k n p a p b p λλλλλλ---=⎧⎪=-⎪=⎨⎪⎪=--⎩ ()det()n p I C λλ=-为C 的特征多项式，n 个零点为矩阵C 的n 个特征值 2 实验程序程序：int ebstq(int n,double b[],double c[],double q[],double eps,int l) {int i,j,k,m,it,u,v;double d,f,h,g,p,r,e,s;c[n-1]=0.0;d=0.0;f=0.0;for (j=0; j<=n-1; j++){it=0;h=eps*(fabs(b[j])+fabs(c[j]));if (h>d){d=h;}m=j;while ((m<=n-1)&&(fabs(c[m])>d)){m=m+1;}if (m!=j){do{if (it==l){printf("fail\n");return(-1);}it=it+1;g=b[j];p=(b[j+1]-g)/(2.0*c[j]); r=sqrt(p*p+1.0);if (p>=0.0){b[j]=c[j]/(p+r);}else{b[j]=c[j]/(p-r);}h=g-b[j];for (i=j+1; i<=n-1; i++) {b=b-h;}f=f+h;p=b[m];e=1.0;s=0.0;for (i=m-1; i>=j; i--){g=e*c;h=e*p;if (fabs(p)>=fabs(c)) {e=c/p;r=sqrt(e*e+1.0);c[i+1]=s*p*r;s=e/r;e=1.0/r;}else{e=p/c;r=sqrt(e*e+1.0);c[i+1]=s*c*r;s=1.0/r;e=e/r;}p=e*b-s*g;b[i+1]=h+s*(e*g+s*b); for (k=0; k<=n-1; k++) {u=k*n+i+1;v=u-1;h=q;q=s*q[v]+e*h;q[v]=e*q[v]-s*h; }}c[j]=s*p;b[j]=e*p;}while (fabs(c[j])>d);}b[j]=b[j]+f;}for (i=0; i<=n-1; i++){k=i; p=b;if (i+1<=n-1){j=i+1;while ((j<=n-1)&&(b[j]<=p)){k=j;p=b[j];j=j+1;}}if (k!=i){b[k]=b;b=p;for (j=0; j<=n-1; j++){u=j*n+i;v=j*n+k;p=q; q=q[v];q[v]=p;}}}return(1);}。
计算广义实对称三对角矩阵特征值问题的分治算法
2002 年第 24 卷第 5 期 Vol 24, No 5, 2002
文章编号: 1007 130X( 2002) 05 0015 03
计算广义实对称三对角矩阵特征值问题 的分治算法
The Divide and Conquer Algorithm for Generalized Symmetric Tridiagonal Eigenvalue Problems
L ) ( x ) = x + n/ (- f ∗( x ) / f ( x ) )
( n - 1) [ ( n - 1) (- f ∗( x ) / f ( x ) ) 2 - n(f +( x ) / f ( x ) ) ] )
其中 f ( x ) 、f ∗( x ) 和 f +( x ) 均采用三项递推式, 且 计算量基本相同。当 x ( i , i + 1) 时, L + ( x ) 比 较接近 i + 1, L_( x ) 比较接近 i 。并且 Laguerre 迭 代是三次收敛的, 割线迭代法的渐进收敛率为 1. 618[ 3] 。显然, 要收敛到相同的精度, Laguerre 迭代 次数比割线迭代法少。但每迭代一次, Laguerre 迭 代算法的计算量显然要比割线法多。那么, 两者 比较, 到底有什么样的关系呢?
T^ - S^ 0
0T
W( )
定义 1[ 4] 对称矩阵 A 的惯性( v, !, ∀) 分别
表示矩阵 A 的小于 0、等于 0 和大于 0 的特征值 个数。
定理 1( Sylvester s 惯性定理[ 1] ) 如果 A Rn # n 是对称的, X Rn# n 是非奇异的, 那么矩阵
16
A 和 XT AX 有相同的惯性。 由引理 1 和定理 1 可得: 定理 2 设 不是( T^ , S^ ) 的特征值, 并令 V
数学实验“对称三对角矩阵特征值的二分法”实验报告(内含matlab程序)
西京学院数学软件实验任务书课程名称数学软件实验班级数0901 学号0912020107 姓名李亚强实验课题对称三对角矩阵特征值的二分法实验目的熟悉对称三对角矩阵特征值的二分法运用Matlab/C/C++/Java/Maple/Mathematica等其中实验要求一种语言完成实验内容对称三对角矩阵特征值的二分法成绩教师实验十四实验报告一、实验名称：对称三对角矩阵特征值的二分法。
二、实验目的：熟悉对称三对角矩阵特征值的二分法。
三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成程序设计。
四、实验内容：%对称三对角矩阵特征值的二分法program ex0001integer nreal x1,x2,x,f1,f2,fx,epsreal,allocatable::d(,e(write(*,*) "Please enter n:"read(*,*) nallocate(d(n),e(n-1))eps=1.0E-6fx=1.0do j=1,nd(j)=-2end dodo j=1,n-1e(j)=1end do!write(*,*) "Please enter array d and e:"!read(*,*) d,ea1=d(1)-e(1)a2=d(1)-2*e(1)a3=d(1)+e(1)a4=d(1)+2*e(2)y1=min(a1,a2)y2=max(a4,a3)y=(y2-y1)/nx1=y1x2=y1+ywrite(*,*)"矩阵的特征值为:" do m=1,ntemp=x210 if(abs(fx)>eps) thenx=(x1+x2)/2f1=MValue(x1,n,d,e)f2=MValue(x2,n,d,e)fx=MValue(x,n,d,e)if (fx*f1>0) thenx1=xelsex2=xend ifgo to 10end ifprint*,"--------------------"write(*,*) xx1=tempx2=temp+yz=(x1+x2)/2fx=MValue(z,n,d,e)end doprint*,"--------------------" containsreal function mValue(x,n,d,e) Integer n,ireal xreal d(n),e(n-1),s(n)S(1)=x-d(1)S(2)=(x-d(1))*(x-d(2))-e(1)**2Do i=3,nS(i)=(x-d(i))*s(i-1)-e(i-1)**2*s(i-2) End domValue=s(n)returnEnd function mValueend。
三对角矩阵fortran语言
三对角矩阵fortran语言什么是三对角矩阵，如何在Fortran语言中实现它，以及它在计算机科学和工程中的应用。
一、什么是三对角矩阵？三对角矩阵是一种具有特殊结构的矩阵。
它的主对角线上的元素是非零的，而其上方和下方的对角线上元素均为零。
换句话说，只有主对角线和两个相邻的对角线上的元素不为零，其余元素均为零。
三对角矩阵可以表示为以下形式：[A B 0 0 …0][C A B 0 …0][0 C A B …0][………………][0 0 …C A B][0 0 …0 C A]其中，A、B和C代表矩阵中的非零元素。
三对角矩阵通常是方阵。
二、如何在Fortran语言中实现三对角矩阵？在Fortran语言中，我们可以使用数组来表示矩阵，并利用循环结构来计算和操作矩阵的元素。
为了实现三对角矩阵，我们可以利用Fortran中的二维数组来存储矩阵，同时使用循环结构来对主对角线和相邻对角线上的元素进行赋值。
下面是一个简单的示例代码：fortranprogram tridiagonal_matriximplicit noneinteger, parameter :: n = 5 ! 矩阵的维度integer :: i, jreal :: A(n,n) ! 存储三对角矩阵的数组! 初始化矩阵do i = 1, ndo j = 1, nA(i, j) = 0.0 ! 将所有元素初始化为零end doend do! 赋值主对角线和相邻对角线上的元素do i = 1, nA(i, i) = 2.0 ! 主对角线上的元素if (i < n) thenA(i, i+1) = 1.0 ! 上方相邻对角线上的元素A(i+1, i) = 1.0 ! 下方相邻对角线上的元素end ifend do! 打印矩阵write(*, '(5F5.2)') (A(i, j), i = 1, n, j = 1, n)end program tridiagonal_matrix通过上述示例代码，我们可以创建一个5x5的三对角矩阵，并将其打印出来。
6.2三对角矩阵法
2
150
305.4.
3
150
319.3
4
150
333.2
5
150
347.0
V2=150 mol/h
平衡常数：
按(6-12)~(6-15)计算常数 A、B、C、D， 得到方程组的矩阵(6-16)的形式:
⎡−150
⎢ ⎢
100
⎢0
⎢ ⎢
0
⎢⎣ 0
244.5 − 344.5
100 0 0
0 325.5 525.5 200
j
Bi, j = −[Vj+1 + ∑ (Fm −Um −Wm ) −V1 + U j + (Vj +Wj )Ki, j ] m=1
Ci, j = Vj+1Ki, j+1 1 ≤ j ≤ N -1
Di, j = −Fj zi, j 1≤ j ≤ N
(6-12)
1≤ j ≤ N (6-13) (6-14) (6-15)
⎥⎢VN−1⎥ ⎢ γ N−2 ⎥
⎢⎣
αN-1 βN-1⎥⎦⎢⎣ VN ⎥⎦ ⎢⎣ γ N−1 ⎥⎦
(6-33)
逐级求解的通式：
Vj
=
γ
j−1
−α jVj−1 β j-1
⑥迭代终止标准：
(6-36)
ε T = ∑[(T j ) k − (T j ) k−1 ]2 ≤ 0.01N
(6-37) (6-38)
对于具有三对角线矩阵的线性方程组， 常用追赶法(或称托玛斯法)求解。该
法仍属高斯消元法。变换(6-16)第一行(式)：
j=1：
Bi,1xi,1 + Ci,1xi,2 = Di,1 ，
矩阵特征值求解
多项式 f () det( A I ) 的根可能对多项式的系数非常敏感.因此,这个方法只
能在理论上是有意义的,实际计算中对一般矩阵是不可行的.首先,若矩阵 A 的阶
数较大,则行列式 det( A I ) 的计算量将非常大;其次,根据 Galois 理论,对于次
数大于四的多项式求根不存在一种通用的方法,基于上述原因,人们只能寻求其
如下分割
Tn T (0) T (1) E
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，通系电1，力过根保管据护线生高0不产中仅工资2艺料22高试2可中卷以资配解料置决试技吊卷术顶要是层求指配，机置对组不电在规气进范设行高备继中进电资行保料空护试载高卷与中问带资题负料2荷试2，下卷而高总且中体可资配保料置障试时2卷，32调需3各控要类试在管验最路；大习对限题设度到备内位进来。行确在调保管整机路使组敷其高设在中过正资程常料1工试中况卷，下安要与全加过，强度并看工且25作尽52下可22都能护可地1关以缩于正小管常故路工障高作高中；中资对资料于料试继试卷电卷连保破接护坏管进范口行围处整，理核或高对者中定对资值某料，些试审异卷核常弯与高扁校中度对资固图料定纸试盒，卷位编工置写况.复进保杂行护设自层备动防与处腐装理跨置，接高尤地中其线资要弯料避曲试免半卷错径调误标试高方中等案资，，料要编试求5写、卷技重电保术要气护交设设装底备备置。4高调、动管中试电作线资高气，敷料中课并设3试资件且、技卷料中拒管术试试调绝路中验卷试动敷包方技作设含案术，技线以来术槽及避、系免管统不架启必等动要多方高项案中方；资式对料，整试为套卷解启突决动然高过停中程机语中。文高因电中此气资，课料电件试力中卷高管电中壁气资薄设料、备试接进卷口行保不调护严试装等工置问作调题并试，且技合进术理行，利过要用关求管运电线行力敷高保设中护技资装术料置。试做线卷到缆技准敷术确设指灵原导活则。。：对对在于于分调差线试动盒过保处程护，中装当高置不中高同资中电料资压试料回卷试路技卷交术调叉问试时题技，，术应作是采为指用调发金试电属人机隔员一板，变进需压行要器隔在组开事在处前发理掌生；握内同图部一纸故线资障槽料时内、，设需强备要电制进回造行路厂外须家部同出电时具源切高高断中中习资资题料料电试试源卷卷，试切线验除缆报从敷告而设与采完相用毕关高，技中要术资进资料行料试检，卷查并主和且要检了保测解护处现装理场置。设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
矩阵特征值求解
矩阵特征值求解的分值算法12组1.1矩阵计算的基本问题(1)求解线性方程组的问题•即给定一个n阶非奇异矩阵A和n维向量b，求一个n维向量X,使得Ax =b (1. 1. 1 )(2)线性最小二乘问题，即给定一个mx n阶矩阵A和m维向量b ,求一个n维向量X,使得|AX -b| =min{ | Ay -比严R n} (1.1.2 )(3)矩阵的特征问题，即给定一个n阶实(复)矩阵A,求它的部分或全部特征值以及对应的特征向量，也就是求解方程(1. 1. 3 )一对解(4 X),其中R(C), x- R n(C n),即A为矩阵A的特征值，X为矩阵Ax = ZxA的属于特征值A的特征向量。
在工程上，矩阵的特征值具有广泛的应用，如大型桥梁或建筑物的振动问题：机械和机件的振动问题；飞机机翼的颤振问题；无线电电子学及光学系统的电磁振动问题;调节系统的自振问题以及声学和超声学系统的振动问题•又如天文、地震、信息系统、经济学中的一些问题都与矩阵的特征值问题密切相关。
在科学上，计算流体力学、统计计算、量子力学、化学工程和网络排队的马尔可夫链模拟等实际问题，最后也都要归结为矩阵的特征值问题.由于特征值问题在许多科学和工程领域中具有广泛的应用，因此对矩阵的特征值问题的求解理论研究算法的开发软件的制作等是当今计算数学和科学与工程计算研究领域的重大课题，国际上这方面的研究工作十分活跃。
1.2矩阵的特征值问题研究现状及算法概述对一个nxn阶实(复)矩阵A,它的特征值问题，即求方程(1.1.3)式的非平凡解，是数值线性代数的一个中心问题•这一问题的内在非线性给计算特征值带来许多计算问题•为了求(1.1.3)式中的A ,—个简单的想法就是显式地求解特征方程det (A 一几I)二0 (121 ) 除非对于个别的特殊矩阵，由于特征方程的系数不能够用稳定的数值方法由行列式的计算来求得，既使能精确计算出特征方程的系数，在有限精度下，其特征多项式f〃)二det(A-ZJ)的根可能对多项式的系数非常敏感能•因此，这个方法只在理论上是有意义的，实际计算中对一般矩阵是不可行的数 _ . _ . 人较大，则行列式det (A -几I)的计算量将非常大；其次，根据•首先，右矩即AfbJ阳数大于四的多项式求根不存在一种通用的方法，基于上述原Galois理论对于次因，人们只能寻求其它途径•因此，如何有效地！精确地求解’矩阵特征值问题，就成为数值线性代数领域的一个中心问题.目前，求解矩阵特征值问题的方法有两大类：一类称为变换方法，另一类称为向量迭代方法•变换方法是直接对原矩阵进行处理，通过一系列相似变换，使之变换成 一个易于求解特征值的形式，如Jacobi 算法,Givens 算法,QR 算法等。
对称矩阵三对角化的混合并行算法设计
对称矩阵三对角化的混合并行算法设计
赵永华;迟学斌;陈江
【期刊名称】《计算机工程》
【年(卷),期】2005(31)22
【摘要】基于Householder转换,给出了稠密对称矩阵三对角化的
MPI+OpenMP混合并行算法.内容集中在SMP集群系统环境下算法的负载平衡、通信开销和性能评价.OpenMP共享内存并行采用了粗粒度方法,解决了MPI算法
中的负载平衡问题,降低了通信开销.在深腾6800上的试验结果表
明,MPI+OpenMP版本比纯MPI版本具有更好的性能和可扩展性.
【总页数】4页(P39-41,53)
【作者】赵永华;迟学斌;陈江
【作者单位】中国科学院计算机网络信息中心超级计算中心,北京100080;中国科
学院软件所,北京100080;中国科学院研究生院,北京100080;德州学院计算机科学系,德州253000;中国科学院计算机网络信息中心超级计算中心,北京100080;中国
科学院计算机网络信息中心超级计算中心,北京100080
【正文语种】中文
【中图分类】TP301.6
【相关文献】
1.解实三对角对称矩阵特征值问题的并行LLT—QR算法 [J], 郭照立;王能超
2.对称稀疏矩阵三对角化并行算法 [J], 谷艺;谷元
3.基于分治法求解对称三对角矩阵特征问题的混合并行实现 [J], 朱京乔;赵永华
4.基于分治法求解对称三对角矩阵特征问题的MPI/Cilk混合并行算法 [J], 朱京乔; 赵永华
5.并行对称矩阵三对角化算法在GPU集群上的有效实现 [J], 刘世芳;赵永华;于天禹;黄荣锋
因版权原因，仅展示原文概要，查看原文内容请购买。
《线性方程组的直接解法及其应用》研究综述
学院：建筑工程学院专业：结构工程组号：16 成绩：报告题目：《线性方程组的直接解法及其应用》研究学院：建工学院专业：结构工程组号：16号成员：xxx学院： 建筑工程学院 专业：结构工程 组号：16 成绩：《线性方程组的直接解法及其应用》研究第一章对象描述一、 《线性方程组的直接解法及其应用》描述在科技、工程、医学、经济等各个领域中，经常遇到求解n 阶线性方程组⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++,,,22112222212*********m n mn m m n n n n b x a x a x a b x a x a x a b x a x a x a （1.1） 的问题.方程组（1.1）的系数),,2,1.(n j i a ij =和右端项),,2,1(n i b i =均为实数，且1b 、2b ，……，n b 不全为零。
方程组（1.1）可简记为b Ax =其中 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=mn m m n n a a a a a a a a a A 212222111211⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n x x x x 21 ， ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=m b b b b 21. 线性方程组的数值解法有两大类，一类是直接法，另一类是迭代法。
本次主要研究的是直接解法。
所谓直接法就是经过有限步算术运算，可求得线性方程组精确解的方法（若计算过程中没有舍入误差）。
但实际计算中由于舍入误差的存在和影响，这种方法也只能求得现行方程组的近似解。
这类算法中最基本的是高斯消元法及其某些变形，它是解决低阶稠密矩阵方程组及某些系数矩阵方程组的有效方法学院： 建筑工程学院 专业：结构工程 组号：16 成绩：二、 《线性方程组的直接解法及其应用》的相关概念1.特征值和特征向量设A 是一个n n ⨯阶实矩阵，若对于数λ，存在非零向量x ，使得x Ax λ=成立。
则称λ是A 的特征值(Characteristic Value)，x 为A 的对应于λ的特征向量(Characteristic Vector)。
三对角系统并行算法的研究概况
地 Ｌ 本地 Ｌ Ｄ 和本地循 环约化 法求解 ， 在 １９ Ｕ、 Ｕ 并 ９ ５年 提 出 基 于矩 阵 乘分 解 的并行 Ｑ Ｒ算 法 。 Ｈ．ｉｈｅ ｅ和 ＡＶ ｎｄ ｒ Ｍ ｃｉｌ ｓ ．ａ ｅ Ｖ ｒｔ ｏｓ改变 Ｗａ ｇ 法的消元次序 。 出了通信 量减少 的算法 。 ｎ算 提 李 晓梅等将 ＨＭｉｈ ｅ和 ＡＶ ｎｄ ｒ ｏｓ算法 中的通信模式 ． ｃｉ ｄ８ ．ａ ｅ ￣ ｒｔ 从 单 向 串行 改 为双 向 并行 ， 出 Ｄ Ｐ算 法 ， 目前 最 好 的 三 对 提 Ｐ 是 角 方 程 组 分 布 式 算 法 之 一 。 ２ｏ ０ ０年 骆 志 刚 等 中 依 据 Ｄ Ｐ算 Ｐ 法， 利用 计 算 与 通 信 重 叠 技 术 ， 少处 理 机 空 闲 时 间 取 得 了更 减 好 的并行效果。此类 算法要求解 Ｐ １ 一 阶缩减系统。 ２不重 叠分解。例如 Ｌ ｗｉ Ｓｍ ｈ算法 、 ｈｏｎ算法 、 ． ａｒ ｅ＆ ａ ｅ Ｊ ｓ ｏ ｏ
（ 林 电子科技 大学计 算科学与数 学 系， 桂 广西 桂 林 ５ １０ ） ４ ０ ４
【 摘 要】 在科学和工程 计算中， 多问题 往往 归结为三对角线性 方程组的求解 ， 并行算 法的研 究具有重要 意义。文章全 面 许 其
总 结 了 当前 求 解三 对 角 线性 方 程 组 的 两类 并 行 算 法 ： 直接 解 法和 迭 代 解 法 ． 介 绍 了其 特 点 。 并
＝ 问题 的提 出
设 三 对 角 线性 方程 组 为
Ａ ＝ ＸＹ （） １
式中： Ｒ ＡＥ 非奇异， ０ ＝＞。ｘ （，， ０Ｙ ｄ－ ， ， １ ； ｌ １ ＝ － Ｉ ｘ ｘ …ｘ
（ ，２…ｙ）。 ｙ Ｙ， ｎ １ 此 系统 在许多算法 中被提出 ．因此研究其高性能并行算
分治策略设计并行算法
分治策略设计并行算法
唐策善;梁维发
【期刊名称】《微电子学与计算机》
【年(卷),期】1990(7)4
【摘要】本文通过给出适应D&C策略的几种典型体系结构,用一些具体例子说明D&C策略在设计并行算法方面的重要性.
【总页数】4页(P17-20)
【关键词】分治法;并行算法;计算机
【作者】唐策善;梁维发
【作者单位】中国科技大学计算机系
【正文语种】中文
【中图分类】TP301.6
【相关文献】
1.基于分治策略的数据库查询优化设计与实现 [J], 颜丽
2.分治策略在归并排序中的算法设计 [J], 李六杏
3.并行算法的好教材──评陈国良教授的《并行算法的设计与分析》一书 [J], 陈凌;
4.基于分治策略的数据库查询优化设计与实现 [J], 颜丽
5.收割机远程监测系统的设计—基于云平台数据挖掘并行算法 [J], 吴东林; 张玉华因版权原因，仅展示原文概要，查看原文内容请购买。
对称三对角矩阵 python
对称三对角矩阵是一种特殊的矩阵结构，它在数学和计算机科学领域中有着重要的应用。
对称三对角矩阵具有许多独特的性质，可以通过Python编程语言来实现对称三对角矩阵的计算和操作。
本文将介绍对称三对角矩阵的基本概念、性质、以及如何使用Python来处理对称三对角矩阵。
一、对称三对角矩阵的基本概念1. 对称三对角矩阵的定义对称三对角矩阵是指正定的对称矩阵和负定的对称矩阵的逆矩阵。
它的主对角线上的元素全部非零，其余非零元素都出现在主对角线的三条带状区域内。
2. 对称三对角矩阵的形式对称三对角矩阵的形式可以表示为：[ a1 b1 0 0 0b1 a2 b2 0 00 b2 a3 b3 00 0 b3 a4 b40 0 0 b4 a5 ]3. 对称三对角矩阵的性质对称三对角矩阵具有许多重要的性质，包括对角线上的元素均为实数、任意两个对角线上的非零元素之间有且仅有一个相邻的非零元素等。
二、对称三对角矩阵的Python实现Python是一种优秀的编程语言，它具有强大的数学计算能力和丰富的科学计算库，可以方便地实现对称三对角矩阵的相关操作。
下面将介绍如何使用Python来实现对称三对角矩阵的计算和操作。
1. 生成对称三对角矩阵可以使用numpy库来生成对称三对角矩阵，可以通过以下代码来实现：```pythonimport numpy as npn = 5a = np.random.rand(n)b = np.random.rand(n-1)A = np.diag(a) + np.diag(b, -1) + np.diag(b, 1)print(A)```2. 对称三对角矩阵的特征值和特征向量可以使用numpy库来计算对称三对角矩阵的特征值和特征向量，可以通过以下代码来实现：```pythonimport numpy as npn = 5a = np.random.rand(n)b = np.random.rand(n-1)A = np.diag(a) + np.diag(b, -1) + np.diag(b, 1)eigval, eigvec = np.linalg.eig(A)print("Eigenvalues:", eigval)print("Eigenvectors:", eigvec)```3. 对称三对角矩阵的LU分解可以使用scipy库来进行对称三对角矩阵的LU分解，可以通过以下代码来实现：```pythonimport numpy as npfrom scipy.linalg import lun = 5a = np.random.rand(n)b = np.random.rand(n-1)A = np.diag(a) + np.diag(b, -1) + np.diag(b, 1)P, L, U = lu(A)print("P:", P)print("L:", L)print("U:", U)```4. 对称三对角矩阵的求解可以使用numpy库来求解对称三对角矩阵的线性方程组，可以通过以下代码来实现：```pythonimport numpy as npn = 5a = np.random.rand(n)b = np.random.rand(n-1)A = np.diag(a) + np.diag(b, -1) + np.diag(b, 1)b = np.random.rand(n)x = np.linalg.solve(A, b)print("Solution:", x)```三、对称三对角矩阵的应用对称三对角矩阵在科学和工程计算中有着重要的应用，例如在有限差分法、有限元法、最小二乘法等领域都有着广泛的应用。
求解实对称三对角矩阵特征值的二分法
求解实对称三对角矩阵特征值的二分法特征值分解是矩阵理论中一个重要且有用的概念，它可以用来对矩阵进行分析、研究和求解，特别是在实物应用中，它通常用于研究一个系统的特性，甚至可以用于求解它。
因此，求解实对称三对角矩阵的特征值是一个非常重要的研究课题。
实对称三对角矩阵是一种特殊的矩阵，它的结构比一般的矩阵要简单得多。
如果矩阵是实对称三对角矩阵，那么它可以被简化为三对角矩阵，其中第一行为主对角线，主对角线上的元素称为主对角元素，其余的元素可以分为左右的非对角线，由于其有序的特点，求解三对角矩阵的特征值比一般矩阵的特征值更容易，可以采用更有效的方法。
求解实对称三对角矩阵特征值的二分法是一种求解实对称三对角矩阵特征值的算法，它可以有效地求解实对称三对角矩阵的特征值。
在二分法中，将实对称三对角矩阵分成两部分，以半空间作为边界，通过比较半空间内的三角矩阵的特征根，计算出一个新的解，再通过比较新的解的特征根，来更新当前的解。
一直这样循环，直到取得特征值为止。
二分法求解实对称三对角矩阵的优点在于它具备收敛性，即它可以以收敛到某一特定特征值为条件，从而有效地求解出矩阵的特征值。
此外，它还比较简单，容易理解，易于实现，不需要使用复杂的数学方法，可以使用计算机来解决实际问题。
然而，二分法存在一些缺点。
首先，该方法的收敛速度比较慢，特别是在求解有行列式的实对称三对角矩阵时，可能会出现过拟合的情况。
其次，它变化一般，将矩阵分成两部分时，可能会遗漏有用的信息。
最后，它也可能不能保证收敛到最优解。
由于实对称三对角矩阵的特殊性，特征值的求解有时可以通过特殊的方法得到简单的公式解。
例如，如果矩阵元素有效且是定值，可以采用直接法求解。
在实际应用中，也可以采用其它更有效的方法，如硬件支持的经典方法和软件支持的算法等。
综上所述，求解实对称三对角矩阵特征值的二分法是一种有用的算法，它可以帮助我们快速求解出实对称三对角矩阵的特征值，在实际应用中，我们可以根据实际情况选择合适的求解方法，以满足我们的需求。
基于分治法求解对称三对角矩阵特征问题的MPI
㊀第52卷第1期郑州大学学报(理学版)Vol.52No.1㊀2020年3月J.Zhengzhou Univ.(Nat.Sci.Ed.)Mar.2020收稿日期:2019-01-14基金项目:国家重点研发计划项目(2017YFB0202202,2016YFB0201302);国家自然科学基金重点项目(91430214);中国科学院 十三五 信息化建设专项项目(XXH13506-405)㊂作者简介:朱京乔(1993 ),男,安徽安庆人,硕士研究生,主要从事并行和分布式计算研究,E-mail:zjqucas@.;通信作者:赵永华(1966 ),男,山东德州人,研究员,主要从事高性能计算和软件㊁并行与分布式处理与算法研究,E-mail:yhzhao@㊂基于分治法求解对称三对角矩阵特征问题的MPI /Cilk 混合并行算法朱京乔1,2,㊀赵永华1(1.中国科学院计算机网络信息中心高性能部㊀北京100190;2.中国科学院大学计算与控制学院㊀北京100049)摘要:对称稠密矩阵特征问题的求解通常转化为三对角矩阵特征问题的求解,基于对称三对角矩阵特征求解的分而治之方法,提出了一种基于消息传递接口(message passing interface,MPI)+Cilk 多任务并行模型的混合求解算法,采用进程间数据并行和进程内多线程任务并行的方法,实现了对分而治之算法中分治阶段和合并阶段的多任务划分和动态调度㊂进程内利用Cilk 任务执行的有向无环图模型,解决了线程级并行的数据依赖和饥饿等待等问题,提高了程序的并行性;进程间通过粗粒度计算任务的划分,优化了MPI 部分的数据通信流程和负载均衡问题㊂数值实验表明,混合并行算法在计算性能和可扩展性方面都要优于纯MPI 并行算法㊂关键词:三对角矩阵;对称特征问题;并行计算;分治算法;Cilk;MPI 中图分类号:TP391㊀㊀㊀㊀㊀文献标志码:A㊀㊀㊀㊀㊀文章编号:1671-6841(2020)01-0033-06DOI :10.13705/j.issn.1671-6841.20190140㊀引言科学计算和工程问题中经常涉及稠密对称矩阵的特征求解问题,通常的处理方法是先通过正交变换把原矩阵化为三对角矩阵,如Householder 或Givens 方法,再求解三对角阵的特征值和特征向量㊂三对角阵的特征求解算法有很多种,常用的有QR 法㊁二分法㊁分而治之方法和多近似健壮表示(multiple relatively robust representations,MRRR)方法等等,本文工作基于分治法展开㊂对称三对角阵特征求解的分治算法,最早由文献[1]提出,文献[2]作出了改进,提高了该算法的求解精度和特征向量的正交性㊂分治算法采用分而治之思想,将原矩阵划分为若干小的子矩阵,先求出子矩阵的特征分解,接着通过修正计算,逐级合并子矩阵的结果,回代得到原矩阵的特征分解㊂分治法具有很强的灵活性,对于大规模三对角矩阵特征问题的并行求解是一种比较理想的算法㊂现有很多模型都实现了针对分治算法的并行化,如基于分布式存储的消息传递接口(message passing in-terface,MPI)进程级并行[3],基于共享内存的openMP 线程级并行,还有基于对称多处理结构(symmetrical multi procession,SMP)的MPI /openMP 混合并行模型[4]等㊂以openMP 为代表的线程级并行模型主要基于数据并行的考量,通过把数据划分到多个核上实现并行,该模型在计算负载均衡的条件下能取得较好的效果㊂但对于分治和递归类的场景,当数据依赖呈复杂网状结构时,则不能很好地实现负载均衡㊂Cilk 是一个基于任务的并行编程模型,为串行程序的并行化提供了向量化㊁视图等机制(http:ʊ /cilk /manual-5.4.6.pdf),通过关键字来声明操作,改造后的代码具有语义化㊁可读性强的特点㊂Cilk 为共享内存系统提供了一个高效的任务窃取调度器,适合用来实现高效的多核任务并行㊂1㊀三对角矩阵特征求解分治算法的基本原理1.1㊀三对角矩阵的划分对任意n 阶实对称三对角矩阵T ,基于秩1修正作如下划分㊂郑州大学学报(理学版)第52卷T =α1β1β1α2β2β2α3⋱⋱⋱βn βnαn éëêêêêêêêùûúúúúúúú⇨T =T 100T 2éëêêùûúú+ββββéëêêùûúú,其中:T 1㊁T 2是三对角子矩阵块,β元素所在矩阵构成原矩阵的修正矩阵㊂为使修正矩阵元素一致,T 1的末位元素和T 2的首元素都减去相同的β㊂经过划分原矩阵变为T =T ᶄ+βzz T ,其中Z =(0, ,0,1,1,0, ,0)T ㊂1的索引代表对原矩阵进行划分的位置㊂对于划分后的矩阵T ᶄ,若规模未达到指定阈值,则继续划分㊂1.2㊀分治法求解过程若划分后的子矩阵规模达到设定阈值,调用QR 或者其他算法,得到子矩阵T 1㊁T 2的特征值分解为T 1=Q 1D 1Q T 1;T 2=Q 2D 2Q T 2,(1)此时T 可以表示为T =Q 1Q 2éëêêùûúú(D +βxx T )[Q T 1Q T 2],(2)其中:D =D 1D 2éëêêùûúú;x T =(q T 1,q T 2),q 1㊁q 2分别是Q 1的最后一行和Q 2的第一行元素构成的向量㊂要求T 的特征分解T =Q ΛQ T ,先求D +βxx T 的特征分解,使得D +βxx T =P ΛP T ,再令Q =Q 1Q 2éëêêùûúú,即可求得T 的特征分解㊂对于D +βxx T 的特征分解,需要根据划分后矩阵的特征值和原矩阵特征值的间隔性[5],通过求解对应的secular 方程1+βði x 2i d i -λ=0,d i =D ii ,(3)得到相应的特征值λ㊂文献[5]讨论了方程(3)根的分布特性和解法,在特征区间内使用牛顿法或拉格朗日法迭代逼近区间内的特征值,线性时间内即可收敛㊂解出所有的特征值后,再通过式(D -λI )-1x [1]求得每个特征值对应的特征向量,回代到式(1)中得到T 的特征分解㊂图1㊀分治算法结构Figure 1㊀The structure of divide and conquer method整个求解过程可看作逻辑上的树型结构,如图1所示㊂在叶子结点进行特征求解,在非叶子结点上进行特征合并㊂1.3㊀求解注意事项降阶:合并过程中如果D 的主对角线上出现了相同的元素或者X 向量中出现0元素时,文献[6]中提出要先进行de-flation 操作,通过Givens 正交旋转变换进行迭代逼近前的预处理,部分特征值和特征向量就能原地得到,从而避免了特征区间内迭代不收敛的情况㊂正交性:如果特征值前后间隔过小,使用原公式求解得到的特征向量会逐渐丢失正交性㊂文献[2]给出了特征向量的改进求法,引入了新的计算方式来修正计算过程中结果的正交性㊂2㊀基于Cilk 的节点内并行2.1㊀节点内并行算法对于单节点内的三对角矩阵块,使用分治法求解主要包含两个过程:调用特征求解算法求出最小划分矩阵的特征分解,逐步合并特征值和特征向量㊂由于每个子矩阵的特征计算以及同级子矩阵之间的合并过程是相互43㊀第1期朱京乔,等:基于分治法求解对称三对角矩阵特征问题的MPI /Cilk 混合并行算法独立的,下面结合Cilk 任务并行机制,给出分治算法在节点内多核环境下的并行化递归实现㊂㊀㊀算法1㊀分治算法基于Cilk 的节点内并行㊂dc _solver(T ):if (size T ==threshold):㊀exec _solve(T );㊀return result;else:㊀cilk _spawn ㊀dc _solver(T 1);㊀cilk _spawn ㊀dc _solver(T 2);㊀cilk _sync;㊀exec _merge();㊀return result;end㊀㊀先判断输入矩阵的规模,如果矩阵规模小于等于指定阈值,调用QR 或其他算法进行特征求解,返回结果;否则执行Cilk 任务划分,划分矩阵T 得子矩阵T 1㊁T 2,并行对T 1㊁T 2递归调用算法1,求出子矩阵的特征分解后,进行排序和迭代逼近等操作,合并子矩阵的特征值和特征向量,得到T 的特征分解,返回结果㊂Cilk 在程序递归执行过程中不断划分新的任务,产生Cilk 线程去处理,并把其分配到空闲的核上,和父线程并行执行㊂Cilk 通过任务窃取(Work-Stealing)的方式来执行调度:当一个核完成自己的任务而其他核还有很多任务未完成,此时会将忙的核上的任务重新分配给空闲的核㊂算法1计算过程中涉及矩阵乘法㊁排序㊁迭代逼近等操作,可以利用Cilk 提供的数据并行机制加速,Cilk 通过数据划分形成一个个可供调度的线程级任务并行执行㊂空闲线程通过工作密取从其他线程获取一部分工作量,Cilk 利用贪心策略来进行任务调度分配,能够避免单核上出现任务过载和饥饿等待的问题,同时能将密取的次数控制在最低水平,减少密取带来的性能开销㊂图2㊀节点内任务并行流程Figure 2㊀Task flow of divide and conquer method inside nodes2.2㊀算法分析一个Cilk 线程可看作在同步㊁生成新线程或返回结果之前所能执行的最长序列化指令集(http:ʊ /cilk-documentation-full)㊂根据此定义,求解过程中,递归调用T 1之前的过程可以看作一个Cilk 线程,记为任务A,递归调用T 2的过程可以看作任务B,同步及合并过程可以看作任务C㊂现假设程序只进行了4次嵌套调用,则全局任务划分后生成的任务流程如图2所示㊂程序从有向无环图的最上面开始执行,对应递归的最外层㊂每一个出度大于1的结点都会派生出新的Cilk 线程,执行新的任务㊂设程序执行过程产生的所有的Cilk 线程数为W (图2中节点数),关键路径长度为S (图2中虚线部分路径长度)㊂则程序的并行性Parallelism 可表示为Parallelism =W /S ㊂(4)㊀㊀程序的执行时间Time 满足Time >max(T (W /P ),T (S )),(5)其中:P 是处理器核数;T (W /P )表示所有任务平均到每个核上运行的时间;T (S )表示关键路径上的任务运行花费的总时间㊂3㊀基于MPI /Cilk 的多节点混合并行3.1㊀混合并行算法实现分治算法的MPI 多节点并行,需要把原始矩阵按图1结构划分成小的子矩阵分发到叶子结点上进行初步计算,然后通过MPI 通信,每棵子树的父结点收集子结点的计算结果后进行汇总,把整理后的结果重新发送到子结点进行合并操作,重复此过程,直到返回到图1中的根结点㊂为了减少进程间通信,充分利用5363郑州大学学报(理学版)第52卷Cilk任务并行机制,需要对原矩阵进行较粗粒度的划分,以榨取单节点处理器的计算性能㊂并且在消息传递时,只使用一个线程作为主控线程,以减少通信开销㊂对于n阶对称三对角矩阵T在N个进程上的特征求解,假设N=2q,n=2p,算法初始时把原矩阵划分为N个子矩阵,每个子矩阵阶数为2p-q,设为2k,下面给出分治算法的在多进程间的混合并行实现㊂算法2㊀分治算法基于MPI/Cilk的节点间并行㊂第一步,对于划分到N个进程上的每个初始子矩阵,调用算法1求出各部分特征值和特征向量㊂第二步,进行p次循环(p对应图1中树的深度),每次循环先把进程分组并确定主控进程来进行结果收发㊂对每个MPI进程,需要根据自己的进程类别,按序选择执行下列步骤:1)组内进程向主控进程发送本进程求得的所有局部特征值,如果本进程是上一轮循环的主控进程,则还要向本轮主控进程发送上一轮计算得到的特征向量矩阵㊂2)主控进程收集属于同一组进程中所有求得的局部特征值和特征向量,对特征值进行排序,并保留排序后的位置索引数组㊂按照1.2节的过程拼接Z1㊁Z2成一个新向量X,根据位置索引将X排序,将排好序的特征值均匀地划发到组内各进程,同时也把向量X传到组内各进程㊂3)组内进程接收主控进程发送过来的排好序的部分特征值和向量X,求解secular方程和特征向量,并发送结果给主控进程㊂4)主控进程更新当前循环的特征向量矩阵,接着执行下一轮循环㊂3.2㊀算法分析在对特征值进行排序后,同时记录排序后各位置元素的原来位置的索引,避免对特征向量重复排序,方便下一轮特征向量的计算㊂算法2中步骤3)和4)中涉及大量计算密集型操作,占据程序执行的主要时间,同样可以通过Cilk数据并行机制划分多个线程级子任务并行执行㊂在第p轮循环中,从0号进程开始,每2p+1个进程分为一组,每隔2p+1个进程被选为主控线程,每个小组内共进行3(2p+1-1)次通信,每轮循环的通信次数是3N(2p+1-1)/2p+1㊂分治算法总的时间复杂度介于O(n2)~O(n3)之间,空间复杂度为O(n2),误差和矩阵规模的关系为O(nε)㊂若单节点内划分的数据粒度适当粗时,就能利用Cilk机制充分发挥多核处理器的计算和调度能力,得到较好的计算性能,同时减少进程间通信次数,降低通信开销㊂4㊀实验结果与分析我们在中科院 元 级超算上针对该混合模型进行了数值实验,实验运行在CPU II计算队列上,节点上搭载的是Intel E5-2680V3芯片,每块芯片包含12个主频为2.5GHz的计算核心㊂MPI程序使用Intel mpi-icpc编译器编译,编译时链接Cilk Plus运行,计算过程中使用了部分MKL库中的函数㊂实验对象是30000阶实对称三对角矩阵,矩阵划分求解的规模阈值设为200阶,最终求出全部特征值和特征向量㊂实验通过环境变量设置使用4和8个线程进行实验,为方便对比,已测出实验在单节点单线程环境下平均串行时间为882.53s㊂表1是纯MPI方法和MPI/Cilk方法在4~64个核参与计算下所用时间比较㊂从表1可以看出,对于同一矩阵,使用相同数目的核参与计算时,MPI/Cilk混合方法所用时间要少于纯MPI方法,性能要优于纯MPI 方法,而8线程的MPI/Cilk方法性能要优于4线程的,这表明8线程的混合模型有更好的负载均衡㊂在较少节点参与计算时,获得加速效果较为明显,这说明对数据进行粗粒度划分时能取得较好性能㊂表2是MPI模型和MPI/Cilk模型程序的加速比对比㊂表2数据表明,MPI/Cilk方法具有比纯MPI方法更高的加速比,而8线程MPI/Cilk方法的加速比开始时和4线程接近,随着可供调度的核数增多,加速效果逐渐超过了4线程方法㊂MPI/Cilk方法的加速比在一开始上升较快,随着进程数的增加,加速效果随之下降,但仍比纯MPI方法变缓更慢,拥有更好的扩展性㊂对于给定规模的矩阵,当参与计算的核数超过一定数量时,计算效果的提升越来越有限,出现这种现象主要是因为计算任务粒度变细和进程间通信开销增加导致的㊂㊀第1期朱京乔,等:基于分治法求解对称三对角矩阵特征问题的MPI/Cilk混合并行算法表3和表4分别给出了Cilk和openMP这两种模型的计算时间和加速比的比较,其中openMP模型采用数据并行方式,通过指导语句对涉及计算的部分进行了并行化㊂从结果可以看出,在参与计算的核数较少时,openMP模型的效果要优于Cilk模型,原因是计算粒度较大时,数据并行占主导地位㊂而随着参与计算的核数增多,Cilk模型的效果逐渐超过openMP模型,此时可供调度的核数增多,Cilk动态调度的优势逐渐体现㊂表1㊀MPI模型和MPI/Cilk模型程序的时间Table1㊀Time of pure MPI method and MPI/Cilk method核数时间/sMPI/Cilk(8线程)MPI/Cilk(4线程)纯MPI4 232.62245.31 8135.91144.26155.62 1674.3681.4392.47 3251.2358.2967.54 6444.7650.5858.87表2㊀MPI模型和MPI/Cilk模型程序的加速比Table2㊀Speedup of pure MPI method and MPI/Cilk 核数加速比MPI/Cilk(8线程)MPI/Cilk(4线程)纯MPI 4 3.79 3.59 8 6.49 6.11 5.67 1611.8610.849.54 3217.2215.1413.07 6419.7117.4414.99表3㊀MPI/openMP模型和MPI/Cilk模型程序的时间Table3㊀Time of MPI/openMP method and MPI/Cilk method核数时间/sMPI/Cilk(8线程)MPI/openMP(8线程)MPI/Cilk(4线程)MPI/openMP(4线程)4 232.62223.15 8135.91124.86144.26138.37 1674.3669.3081.4377.21 3251.2355.1358.2962.43 6444.7648.5750.5855.28表4㊀MPI/openMP模型和MPI/Cilk模型程序的加速比Table4㊀Speedup of MPI/openMP method and MPI/Cilk method核数加速比MPI/Cilk(8线程)MPI/openMP(8线程)MPI/Cilk(4线程)MPI/openMP(4线程)4 3.79 3.958 6.497.07 6.11 6.381611.8612.7310.8411.433217.2216.0115.1414.146419.7118.1717.4415.965㊀总结本文针对实对称三对角矩阵特征求解的分治算法,提出了一种基于MPI/Cilk的混合并行实现算法㊂在节点间,通过粗粒度计算任务的划分,在相同的求解效率下减少了所需要的计算进程和进程间的通信次数㊂在节点内,利用Cilk的任务并行机制提高CPU的利用率,改善了负载均衡,提高了并行度㊂实验结果体现了该算法的性能㊂该算法还可进一步研究,如运用数据并行机制加速特征合并过程,研究Cilk线程启动数和数据划分粒度对性能的影响以及同openMP的任务并行机制对比等㊂参考文献:[1]㊀CUPPEN J J M.A divide and conquer method for the symmetric tridiagonal eigenproblem[J].Numerische mathematik,1980,7383郑州大学学报(理学版)第52卷36(2):177-195.[2]㊀GU M,EISENSTAT S C.A stable and efficient algorithm for the rank-one modification of the symmetric eigenproblem[J].SI-AM journal on matrix analysis and applications,1994,15(4):1266-1276.[3]㊀TISSEUR F,DONGARRA J.A parallel divide and conquer algorithm for the symmetric eigenvalue problem on distributed mem-ory architectures[J].SIAM journal on scientific computing,1999,20(6):2223-2236.[4]㊀ZHAO Y H,CHEN J,CHI X B.Solving the symmetric tridiagonal eigenproblem using MPI/OpenMP hybrid parallelization[M].Heidelberg:Springer,2005:164-173.[5]㊀AMMAR G S,REICHEL L,SORENSEN D C.An implementation of a divide and conquer algorithm for the unitary eigen prob-lem[J].ACM transactions on mathematical software,1992,18(3):292-307.[6]㊀MELMAN A.Numerical solution of a secular equation[J].Numerische mathematik,1995,69(4):483-493.[7]㊀王进,晏世凯,高延雨,等.基于MPI的ML-kNN算法并行[J].郑州大学学报(理学版),2018,50(3):34-38.WANG J,YAN S K,GAO Y Y,et al.Parallelization of ML-kNN based on MPI[J].Journal of Zhengzhou university(natural science edition),2018,50(3):34-38.A Hybrid Parallel Algorithm Based on Divided and Conquer MethodUsing MPI/Cilk for the Symmetric Tridiagonal EigenproblemZHU Jingqiao1,2,ZHAO Yonghua1(1.Department of High Performance Computing,Computer Network Information Center,ChineseAcademy of Sciences,Beijing100190,China;2.School of Computer Science and Technology,University of Chinese Academy of Sciences,Beijing100049,China) Abstract:Divide and conquer method was a good algorithm with high performance used for solving tridi-agonal matrix eigenproblems,while with the enlargement of the scale,computing efficiency and storage limitation became the bottleneck.A new eigenproblem solver algorithm based on a hybrid parallel para-digm with MPI/Cilk was proposed to optimize the divide and conquer algorithm both at data and task lev-els.In order to take good advantage of multi-core CPU,task-based parallelization mechanism was intro-duced inside computing nodes,which solved the problem in data dependence and thread starvation by di-rected acyclic graph model.The overhead of data communication among MPI nodes was optimized by coarse-grained partition of tasks,which also improved load balance.The numerical test was carried out and the result was compared with the pure MPI and MPI/openMP parallel algorithm,which showed the performance and efficiency of the algorithm.Key words:tridiagonal matrix;symmetric eigenproblem;parallel computing;DC method;Cilk;MPI(责任编辑:方惠敏)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
基于分治法求解对称三对角矩阵特征问题的混合并行实现①朱京乔1,2, 赵永华11(中国科学院 计算机网络信息中心, 北京 100190)2(中国科学院大学, 北京 100049)通讯作者: 赵永华, E-mail: yhzhao@摘 要: 基于对称三对角矩阵特征求解的分而治之方法, 提出了一种改进的使用MPI/Cilk 模型求解的混合并行实现, 结合节点间数据并行和节点内多任务并行, 实现了对分治算法中分治阶段和合并阶段的多任务划分和动态调度. 节点内利用Cilk 任务并行模型解决了线程级并行的数据依赖和饥饿等待等问题, 提高了并行性; 节点间通过改进合并过程中的通信流程, 使组内进程间只进行互补的数据交换, 降低了通信开销. 数值实验体现了该混合并行算法在计算效率和扩展性方面的优势.关键词: 并行计算; 对称特征问题; 分治算法; Cilk引用格式: 朱京乔,赵永华.基于分治法求解对称三对角矩阵特征问题的混合并行实现.计算机系统应用,2019,28(9):246–250. /1003-3254/7078.htmlHybrid Parallel Algorithm Using MPI/Cilk for Symmetric Tridiagonal EigenproblemsZHU Jing-Qiao 1,2, ZHAO Yong-Hua 11(Computer Network Information Center, Chinese Academy of Sciences, Beijing 100190, China)2(University of Chinese Academy of Sciences, Beijing 100049, China)Abstract : Divide and conquer algorithm is widely used for tridiagonal matrix eigenproblems while computing efficiency and storage limitation are always bottlenecks for large scale problems. In this study, the proposed eigenproblem algorithm based on hybrid parallel paradigm with MPI/Cilk optimizes the divide and conquer algorithm both at data and task levels.The introduced task-based parallelization mechanism inside computing nodes solves the problem in data dependence and thread starvation by directed acyclic graph model. By coarse-grained partition of tasks the overhead of data communication among MPI nodes is also optimized, which helps to improve load balance. The numerical test is carried out and the result is compared with the pure MPI and MPI/openMP parallel algorithm, which shows the performance and efficiency of the algorithm.Key words : parallel computing; symmetric eigenproblem; DC method; Cilk对称矩阵的特征求解问题在科学计算领域普遍存在, 包括计算物理、计算化学、结构力学、纳米材料等前沿学科. 通常的处理方法是先通过正交变换把原矩阵化为三对角矩阵, 如使用Householder 或Givens 方法, 再求解三对角阵的特征值和特征向量, 然后映射回原矩阵的特征值和特征向量. 三对角阵的特征求解计算机系统应用 ISSN 1003-3254, CODEN CSAOBNE-mail: csa@ Computer Systems & Applications,2019,28(9):246−250 [doi: 10.15888/ki.csa.007078] ©中国科学院软件研究所版权所有.Tel: +86-10-62661041① 基金项目: 国家重点研发计划(2017YFB0202202, 2016YFB0201302); 中国科学院“十三五”信息化建设专项(XXH13506-405)Foundation item: National Key Research and Development Program of China (2017YFB0202202, 2016YFB0201302); CAS Special Fund for Informatization Construction in 13th Five-Year Plan (XXH13506-405)收稿时间: 2019-03-05; 修改时间: 2019-04-02; 采用时间: 2019-04-10; csa 在线出版时间: 2019-09-05算法有很多种, 常见的有QR法、二分法、MRRR方法和分而治之方法等等. 由于QR法和二分法求解特征向量时往往需要显式的正交过程, 时间复杂度为O(n3), 常用于中小规模矩阵的特征求解. 而分治法和MRRR法能避免这一过程, 分治法的平均时间复杂度仅为O(n2.3), 其分而治之思想对于大规模矩阵特征问题的并行求解求解有着天然的优势.分而治之求解对称三对角矩阵特征问题的算法最早由Cuppen提出[1], 后M. Gu和S. C. Eisenstat等人作出改进, 提高了该算法求解的特征向量的正交性[2], 使得该算法能够在实际求解中应用. 该算法采用分而治之思想, 将原矩阵划分为若干子矩阵, 先求出子矩阵的特征分解, 再通过修正计算, 逐级合并子矩阵的结果,回代得到原矩阵的特征分解. 分治法具有很强的灵活性, 适合大规模三对角矩阵特征问题求解的并行化求解.分治算法的并行化工作很早就开始展开[3–5], 如基于分布式存储的MPI进程级并行, 基于共享内存的openMP线程级并行, 还有基于SMP的MPI/openMP 混合并行等模型. 并行化的实现主要基于数据并行的考量, 通过划分数据到多个核或节点上计算, 在负载均衡的条件下能取得较好的并行效果. 但对于分治和递归类的场景, 当数据依赖呈现复杂的网状结构时, 难以实现理想的负载均衡效果.基于任务的并行编程模型按照功能将一个应用划分成多个能并行执行的模块, 可以是细粒度的, 也可以是粗粒度的. 常见的任务并行模型有TBB,Cilk和SMPSS等, 使用这些模型时无需关心任务的具体分配和调度细节.Cilk作为一个基于任务的并行编程模型, 为共享内存系统提供了一个高效的任务窃取调度器, 适合用来为分治和递归类算法实现高效的多核任务并行. 对于并行编程, Cilk为串行程序的并行化提供了向量化、锁、视图等机制[6], 通过关键字来声明操作, 改造后的代码具有语义化、可读性强的特点.1 分治算法求解的基本原理1.1 矩阵划分对任意n阶实对称三对角矩阵T, 基于秩1修正从第k行作如下划分:其中, T'1, T'2是三对角子矩阵块, e k代表第k个元素为1的单位向量, 元素β是T的第k个副对角元素. 为使修正矩阵元素一致, 对T1的末位元素和T2的首元素都减去了相同的元素.分治过程一般采用折半划分,对于划分后的矩阵T', 若规模未达到指定阈值, 则继续划分.1.2 分治法求解过程若划分后的子矩阵规模达到设定阈值, 调用QR 或者其他算法得到子矩阵T’1, T’2 的特征值分解:令修正向量则矩阵T的特征分解如下:令D+βzz T=W, 因为diag(Q1, Q2)为正交矩阵, 所以T和W相似有相同的特征值, 求T的特征分解T=QΛQ T, 可以通过先求出W的特征分解: W=PΛP T, 再令Q=diag(Q1, Q2)P, 即可求得T的特征分解.对于W的特征分解, 需要根据划分后矩阵的特征值和原矩阵特征值之间的间隔性[7], 通过求解对应的secular方程(3):得到相应的的特征值λ. 文献[6]讨论了方程(3)根的分布特性和解法, 在特征区间内使用牛顿法或拉格朗日法迭代逼近区间内的特征值, 线性时间内即可收敛. 解出所有的特征值后, 再通过式(D−λI)−1x[1]求得每个特征值对应的特征向量, 回代到式(1)中得到T的特征分解.整个求解过程在逻辑上可看作树型结构, 在叶子结点上进行特征求解, 在非叶子结点上进行特征合并.2019 年 第 28 卷 第 9 期计算机系统应用1.3 注意事项文献[5]中指出合并过程中如果D的主对角线上出现了相同的元素或者z向量中出现0元素时, 在迭代前进行deflation操作可以避免特征区间内迭代不收敛的情况. 通过Givens正交旋转变换, 部分特征值和特征向量就能原地得到. 如果求得的特征值前后过于接近, 按原方法求得的特征向量会丢失正交性, 文献[8,9]讨论了特征值间隔过小时的特征向量的求法, 避免了显式正交过程. 文献[2]根据特征问题的反问题给出了特征向量的改进求法, 保证了计算结果的正交性.2 分治算法的并行化2.1 节点内并行实现单节点内求解三对角矩阵块的特征分解主要包含两个过程, 分别是求出最小划分矩阵的特征分解以及逐步合并特征值和特征向量. 由于每个子矩阵的特征计算以及同级子矩阵之间的合并过程是相互独立的,下面结合Cilk模型给出递归实现的分治算法在单节点多核环境下的并行化算法.算法1. 分治算法基于Cilk的节点内并行1) 判断输入矩阵的规模, 如果矩阵规模小于等于设定 阈值, 执行2);否则执行3);2) 调用特征求解算法进行特征分解, 返回结果;3) 划分矩阵T得子矩阵T1, T2, 通过Cilk执行任务划分, 并行地对T1, T2递归调用步骤1), 求出子矩阵的特征值和特征向量后, 进行排序和迭代逼近等操作合并子矩阵的特征值和特征向量, 得到T的特征分解, 返回结果.Cilk在程序递归执行过程中不断划分新的任务,产生Cilk线程去处理, 并把其分配到空闲的核上, 和父线程并行执行. 一个Cilk任务是在同步、生成新线程或返回结果之前执行的最长序列化指令集或代码段[10].求解过程中, 递归调用T1之前的过程可看作任务A, 递归调用T2的过程可看作任务B, 同步及合并过程可看作任务C. 现假设程序只进行了四次嵌套调用, 则全局任务划分后生成的任务流如图1所示.2.2 多节点混合并行实现实现分治算法的多节点并行, 需要先把原始矩阵按逻辑树结构划分成小的子矩阵分发到叶子结点上进行初步计算, 然后通过MPI通信汇总每棵子树的计算结果, 把整理后的结果重新发送到子结点进行合并操作, 重复此过程, 直到返回到树的根结点. 为了减少进程间通信, 充分利用Cilk任务并行机制, 需要对原矩阵进行较粗粒度的划分, 以榨取单节点处理器的计算性能.图1 分治算法的任务并行流程对于n阶对称三对角矩阵T在N个进程上的特征求解, 假设N=2k, n=2j, 算法初始时把原矩阵划分为N个子矩阵, 每个子矩阵阶数为2j–k, 下面给出分治算法在多节点上的混合并行实现.算法2. 分治算法基于MPI/Cilk的节点间并行首先, 对于划分到N个进程上的每个初始子矩阵, 调用算法1求出各部分特征值和特征向量;接着, 进行p轮循环(p对应逻辑树的高度). 每轮循环先把进程分组并确定主控进程. 每个MPI进程根据自己的进程类别按序选择执行下列步骤:1) 组内进程通过互补通信交换求得的局部特征向量矩阵, 接着向主控进程发送本进程求得的局部特征值和拼接的修正向量;2) 主控进程收集所有组内进程求得的局部特征值进行排序, 并保留排序后的位置索引数组, 并根据此数组将修正向量排序, 将排好序的特征值和修正向量划分到组内各进程;3) 组内进程接收主控进程发送过来的排好序的部分特征值和修正向量, 求解secular方程和特征向量;4) 各组内进程按列执行分块矩阵乘法, 更新当前轮循环的部分特征向量矩阵, 接着执行下一轮循环.2.3 算法分析图1最上方左面对应递归的最外层, 每个出度大于1的结点都会派生出新的Cilk线程, 执行新的任务.程序执行过程产生的所有的Cilk线程数为图1中节点数, 关键路径长度为图1中虚线部分路径长度. 程序的并行性为Cilk线程数与关键路径长度之比, 执行时间大于等于所有任务平均到每个核上运行的时间和关键路径上任务运行花费的总时间.算法2中在对特征值进行排序后, 同时记录排序后各位置元素的原来位置的索引, 避免对特征向量重复排序, 方便下一轮特征向量的计算. 在第p轮循环中,从0号进程开始每2p+1个进程分为一组, 每隔2p+1个进程被选为选为主控线程. 更新局部特征向量矩阵时,计算机系统应用2019 年 第 28 卷 第 9 期为了减少数据通信开销, 避免主控进程进行统一收发,使组内进程间只进行必要的局部特征向量的互补交换,然后按列执行分块矩阵乘法, 分别计算更新的特征向量矩阵的部分, 通信量为O(n2/N).算法1和算法2中的合并过程涉及大量计算密集型操作, 占据程序执行的主要时间, 可以使用Cilk提供的数据并行机制加速. Cilk通过数据划分形成一个个可供调度的线程级任务并行执行, 通过基于贪心策略的任务窃取方式执行调度: 当一个核完成自己的任务而其它核还有很多任务未完成, 此时会将忙的核上的任务重新分配给空闲的核.由于在迭代计算不同特征值的近似值时, 所执行的迭代次数可能不同, 均匀分配任务可能导致线程间的负载不平衡. Cilk通过工作密取机制从其他线程获取一部分工作量, 能够避免单核上出现任务过载和饥饿等待的问题, 同时能将密取的次数控制在最低水平,减少密取带来的性能开销. 若单节点内划分的数据粒度适当, 就能利用Cilk机制充分发挥多核处理器的计算和调度能力, 得到较好的负载均衡.3 实验结果与分析我们在中科院“元”超算上的cpuII计算队列进行了数值实验, 节点上搭载的是Intel E5-2680 V3芯片,每块芯片包含12颗主频为2.5 GHz的计算核心. MPI 程序使用Intel mpiicpc编译并链接到Cilk Plus运行时.实验对象是30 000阶主对角元素是4, 副对角元素是1实对称三对角矩阵. 矩阵划分求解的规模阈值设为200阶, 实验最终求出全部特征值和特征向量(通过环境变量设置每个节点使用2~12个Cilk线程进行实验).表1是MPI模型与MPI/Cilk混合模型在4~128个核、2~12个Cilk线程参与计算下所用时间汇总. 从表中可以看出, 相同计算条件下MPI/Cilk混合并行计算时间要少于纯MPI方法, 而随着参与计算的Cilk线程数的增加, 计算用时也越来越少, 这表明当计算任务密集时, 计算效率随着可供调度的线程的增加而增加.当参与计算的Cilk线程数超过10时, 对于当前规模的矩阵计算获得的加速提升效果逐渐变得有限, 这是因为Cilk是基于贪心策略执行调度的, 当任务粒度较小时不会频繁的执行任务窃取, 从而减小系统开销. 而在较少节点参与计算时, 获得的加速效果较为明显, 这说明对数据进行粗粒度划分时能取得较好性能.图2是MPI模型和MPI/Cilk模型在不同参数下的加速比变化趋势曲线(已测出在单节点单线程环境下的平均串行时间为880秒左右). 可以看出, 随着参与计算的核数和Cilk线程数增加, 三条曲线的上升趋势都逐渐变得平缓, 出现这种现象的原因主要为计算任务粒度变细和进程间通信开销增加导致的. MPI/Cilk 方法的曲线在一开始上升较快, 随着进程数的增加, 计算任务粒度逐渐减小, 加速效果也随之下降, 但仍比纯MPI方法变缓得更慢. 这表明MPI/Cilk方法的并行效果要优于纯MPI方法, 拥有更好的可扩展性. 随着可供调度的核数增多, 使用更多的Cilk线程数参与计算能获得更好的计算性能, Cilk动态调度的优势也渐渐体现了出来.表2给出了MPI/Cilk模型(使用8个Cilk线程)同ELPA、ScaLAPACK软件包实现的分治算法求解30000阶矩阵特征的计算时间比较. E L P A和ScaLAPACK是求解线性代数问题领域中应用广泛的并行软件包. 从表中数据可以得出, 当参与计算的核数较少时, MPI/Cilk方法与ELPA、Scalapack中分治算法计算用时接近, 随着参与计算的核数增多, MPI/Cilk 方法、ELPA方法与ScaLAPACK方法的计算性能逐渐拉开, 但本文的MPI/Cilk方法与ELPA方法的扩展性仍然存在着一定差距. 这是由于ELPA按照块级循环分布矩阵, 实现了对计算流程的更细粒度的控制.表1 MPI模型和MPI/Cilk模型运行时间核数时间(秒)MPI/Cilk(12 threads)MPI/Cilk(10 threads)MPI/Cilk(8 threads)MPI/Cilk(6 threads)MPI/Cilk(4 threads)MPI/Cilk(2 threads)PureMPI4----229.14234.72241.47 8--123.72130.28139.26143.86155.62 1664.2166.3868.6772.8977.2782.2692.47 3243.1545.7949.9054.0458.2363.3969.54 6437.4739.0341.7643.3846.6852.4459.71 12832.1734.8535.4237.6641.5145.3552.12 2019 年 第 28 卷 第 9 期计算机系统应用25201510548163264128计算核数Pure MPI mode1MPI/Cilk with 2 threads MPI/Cilk with 4 threads MPI/Cilk with 6 threads MPI/Cilk with 8 threads MPI/Cilk with 10 threadsMPI/Cilk with 12 threads 图2 MPI 模型和MPI/Cilk 模型的加速比曲线表2 MPI/Cilk 方法与ELPA 、ScaLAPACK 求解器时间对比核数时间(秒)ELPA DC solver MPI/Cilk(8 threads)Scalapack DC solver4214.17-269.968119.66123.72171.341663.4768.6794.513242.9349.9062.296430.9441.7651.8312822.1235.4244.364 结论与展望本文针对对称三对角矩阵特征求解的分治算法,提出了一种改进的基于MPI/Cilk 的混合并行实现. 在节点间, 进行粗粒度计算任务的划分, 更新局部特征矩阵时通过弱化主控进程的中心作用, 使组内进程之间只进行必要的互补数据交换, 避免了统一收发流程, 减少了进程间的通信开销; 在节点内, 利用Cilk 的任务并行机制提高了CPU 的利用率和特征求解过程中的负载均衡性. 实验结果体现了该算法的性能. 本文实现的分治算法与最新的ELPA 软件包仍然存在一定的性能差距, 还有可以提升的空间. 此外, Cilk 线程启动数和数据划分粒度对计算性能的影响以及同openMP 等其它任务并行模型的效果对比等等是值得进行下一步研究的方向.参考文献Cuppen JJM. A divide and conquer method for thesymmetric tridiagonal eigenproblem. NumerischeMathematik, 1980, 36(2): 177–195. [doi: 10.1007/BF01396757]1Gu M, Eisenstat SC. A stable and efficient algorithm for therank-one modification of the symmetric eigenproblem. SIAM Journal on Matrix Analysis and Applications, 1994, 15(4):1266–1276. [doi: 10.1137/S089547989223924X ]2Tisseur F, Dongarra J. A parallel divide and conqueralgorithm for the symmetric eigenvalue problem on distributed memory architectures. SIAM Journal on Scientific Computing, 1998, 20(6): 2223–2236.3Zhao YH, Chen J, Chi XB. Solving the symmetrictridiagonal eigenproblem using MPI/OpenMP hybrid parallelization. Cao J, Nejdl W, Xu M. Lecture Notes in Computer Science. Berlin: Springer, 2005. 3756. 164–173.4Ammar GS, Reichel L, Sorensen DC. An implementation ofa divide and conquer algorithm for the unitary eigen problem. ACM Transactions on Mathematical Software,1992, 18(3): 292–307. [doi: 10.1145/131766.131770]5Intel Cilk Plus. Intel Cilk reference manual anddocumentation, https:///cilk-plus-tutoriall .6Melman A. Numerical solution of a secular equation.Numerische Mathematik, 1995, 69(4): 483–493. [doi:10.1007/s002110050104]7Dhillon IS, Parlett BN. Orthogonal eigenvectors and relativegaps. SIAM Journal on Matrix Analysis and Applications,2003, 25(3): 858–899. [doi: 10.1137/S0895479800370111]8Dhillon IS, Parlett BN. Multiple representations to computeorthogonal eigenvectors of symmetric tridiagonal matrices.Linear Algebra and its Applications, 2004, 387: 1–28. [doi:10.1016/a.2003.12.028]9Supercomputing Technologies Group, MIT Laboratory forComputer Science. Cilk-5.4.6 reference manual. Cambridge,Massachusetts, USA: Massachusetts Institute of Technology./cilk/manual-5.4.6.pdf .10计算机系统应用2019 年 第 28 卷 第 9 期。