对块三对角M矩阵的一个并行不完全分解预条件算法

数值计算_矩阵的三角分解算法矩阵的三角分解是一种将矩阵分解为下三角矩阵和上三角矩阵的方法。

三角分解在数学和计算机科学中都有广泛的应用，特别是在线性代数、数值计算和优化问题中。

在本篇文章中，我们将介绍几种常见的矩阵三角分解算法。

一、LU分解LU分解是矩阵三角分解中最常见的一种方法。

它将一个矩阵分解为一个下三角矩阵L和一个上三角矩阵U，使得原始矩阵A可以表示为A=LU。

其中，L矩阵的主对角线元素全为1，而U矩阵的主对角线元素是A矩阵的主对角线元素。

实际上，LU分解可以看作是高斯消元法的矩阵形式。

在进行LU分解时，我们可以通过对原始矩阵A进行一系列的行变换来得到上三角矩阵U。

同时，我们可以记录每一次行变换的乘积以及主元元素的倒数，从而得到下三角矩阵L。

因此，LU分解可以通过高斯消元法来直接实现。

二、Cholesky分解Cholesky分解是一种仅适用于对称正定矩阵的三角分解方法。

它将一个对称正定矩阵分解为一个下三角矩阵L和其转置矩阵的乘积，即A=LL^T。

Cholesky分解非常有效率，尤其适用于解线性方程组和进行矩阵的逆运算。

由于分解结果是一个下三角矩阵，因此Cholesky分解可以减少计算量并提高计算速度。

三、QR分解QR分解是一种将矩阵分解为一个正交矩阵Q和一个上三角矩阵R的方法，即A=QR。

其中，Q矩阵是正交矩阵，其列向量是正交的，而R矩阵是上三角矩阵。

QR分解可以看作是对矩阵A进行一系列的正交变换，使其变为上三角形式。

其中，每一次正交变换可以通过Givens旋转来实现，即通过矩阵的乘积来实现矩阵的旋转。

QR分解在多元线性回归分析、奇异值分解和特征值分解等领域有广泛的应用。

四、LUP分解LUP分解是LU分解的一个变种，并增加了行交换的步骤。

LUP分解将一个矩阵分解为一个下三角矩阵L、一个上三角矩阵U和一个置换矩阵P，使得PA=LU。

其中，L和U的构造方式与LU分解相同，而置换矩阵P是一个与单位矩阵相似的矩阵，用于记录行交换的信息。

题目：三对角矩阵追赶法在Matlab中的应用一、概述三对角矩阵是一类特殊的矩阵，在科学计算领域中有着广泛的应用。

而追赶法（也称为托马斯算法）是解三对角线性方程组的一种高效方法。

在Matlab中，可以通过内置函数或者编写自定义代码来实现三对角矩阵追赶法，本文将介绍该方法的原理及在Matlab中的应用。

二、三对角矩阵及其特点三对角矩阵是指除了主对角线外，只有上下相邻两条对角线上有非零元素的矩阵。

具体来说，三对角矩阵的形式如下：```a1 b1 0 0 0c1 a2 b2 0 00 c2 a3 b3 00 0 c3 a4 b40 0 0 c4 a5```其中a1, a2, ..., an为主对角线上的元素，b1, b2, ..., bn-1为上对角线上的元素，c2, c3, ...,为下对角线上的元素。

三对角矩阵具有以下特点：1. 矩阵的维数较大时，存储和计算均占据较少的空间和时间；2. 存在高效的解法，例如追赶法；3. 在有限差分法、对流扩散方程等实际问题中有广泛应用。

三、追赶法原理及算法流程追赶法是一种用于解三对角线性方程组的常用方法，其原理基于LU 分解和回代求解。

算法流程如下：1. 首先对三对角矩阵进行LU分解，将其分解为一个单位下三角矩阵L和一个上三角矩阵U；2. 然后进行前代和回代，分别求解Ly=d和Ux=y。

这两个步骤可以通过追赶法来高效求解。

四、Matlab中的三对角矩阵追赶法实现在Matlab中，可以通过编写自定义代码来实现三对角矩阵追赶法，也可以利用内置函数进行求解。

下面分别介绍这两种方法：1. 自定义代码实现可以编写自定义函数来实现三对角矩阵追赶法。

下面以一个简单的示例代码来说明实现过程：```matlabfunction x = tridiagonal_solver(a, b, c, d)n = length(d);for i = 2:nm = a(i) / b(i-1);b(i) = b(i) - m * c(i-1);d(i) = d(i) - m * d(i-1);endx(n) = d(n) / b(n);for i = n-1:-1:1x(i) = (d(i) - c(i) * x(i+1)) / b(i);endend```通过上述代码，可以实现对三对角矩阵的追赶法求解。

分块对角矩阵法分块对角矩阵（BD）法是一种矩阵算法，用于进行高效的矩阵运算。

分块对角矩阵法利用对角矩阵和其它矩阵的相乘性质，将矩阵拆分成多个矩阵块，并对每个块矩阵进行单独的计算。

这种方法不仅可以显著减少计算量，还可以提高计算效率，从而大大加快矩阵运算的速度。

分块对角矩阵法是一种适用于各种不同类型的矩阵的通用方法，可以用于求解线性方程组、计算特征值和特征向量、求解最小二乘问题等。

该算法的优点是非常明显的：它提高了计算速度，同时不会降低精度。

下面我们来了解一下这个算法的原理和具体实现方法。

1. 原理分块对角矩阵法的基本思想是将矩阵分成对角块和非对角块两种类型，并对这些块进行分块处理。

其中对角块是由矩阵的主对角线上的元素组成的子矩阵；非对角块则是由其它元素组成的子矩阵。

对角块可以直接进行矩阵的乘法运算，而非对角块可以采用一些优化技巧，如基于矩阵的分解和重组的方法，来加速计算。

分块对角矩阵法的优点在于它可以显著减少计算量。

这是因为对角块之间的乘积是一个对角矩阵，因此我们只需要计算对角块之间的乘积即可；而非对角块之间的乘积可以用一些更快的方法来计算。

这个方法的核心就是将大矩阵的计算问题分解成若干个小问题，从而使得计算速度变得更快。

2. 实现方法为了实现分块对角矩阵法，我们需要先将大矩阵分成若干个小矩阵块。

一般而言，对角块的大小是 $n_b\times n_b$，非对角块的大小是 $n_b \times m$ 或 $m \times n_b$，其中 $n_b$ 是块矩阵的大小， $m$ 是非对角块的长度。

接下来，我们需要将矩阵分解成块矩阵的形式，并用一个分块矩阵来表示它们。

分块矩阵 $A$ 的形式为：$$A = \begin{bmatrix}A_{11} & A_{12} & \cdots & A_{1k} \\A_{21} & A_{22} & \cdots & A_{2k} \\\vdots & \vdots & \ddots & \vdots \\A_{k1} & A_{k2} & \cdots & A_{kk}\end{bmatrix},$$接下来，我们需要将块矩阵乘到一起，得到大矩阵的乘积。

三对角矩阵计算1997年6胃一社E擘高等学校计算数学奢芊盱第z期,聊\/P,三对角矩阵计算唐达(上海电机每科学校数学系,上海200240)TRIDIAGONALMA TRIXCALCULUSTangDa(ShanghaiElec~iealMachineryCollege)AbstractInthispaper,weeoas/derthefollow~.ngeompudngmethodoftridiag onalmatrix:(1)decouplingoflinearequation;(2)determinationofdeterminantandeigenvalue—eigenvectorrelation(3)matrixinversion. Thealgorithminthispaperpossessesatleastofleoftheadvantagessuchasfew- eroperating,lessmemorycapacity,higherprecisionorprogrammingsimplifieomparwiththecurrentmethods.KeywordsTrldiagonalmatrix,~‟-vector,Z-~Ctordecoupling,inversematri x?AMS(1g91)Subjectelasslflcatlons15A57.中圈法分类号O151.21引言在效值计算中,有许多问题最后归结为三对角矩阵的计算,因此研究它们的计算方法是有意义的.此外,有些三对角阵的计算方法可以做为带状阵计算的借鉴.本文讨论三对角线性方程组的解耦算法.矩阵的LR-分解,求行列式,Jacobi矩阵的特征值与特征向量的关系以及三对角阵求逆等方面的问题.与现有的算法比较,本文的算法具有计算量或存贮量较少,或计算精度较高,或编程较简单等某些特点.设A为阶非奇实三对角阵:收藕日期:1995一O5—25.吖七,f,)-/,,r¨,,,,●J●一,f唐选;三对角矩阵计算第2期￡=(,,…,).当c≠0(1,2|..?,n一1)时,向量t按下述递推关系算得ft.一口0=1或不为零之实数),t2=一bitl/cI,(2)【t=一(乌一l一l+61￡‟一I)16-1(1=3,4...-,n).这样的向量称之为对应的)f向量.令At一2.(3)式中==(,≈,…,),称为对应的)z向量.由(2)式知,除≠0外,z向量的其余分量均为零.我们称z向量中不为零的分量为尾分量.上面定义的￡向量是右乘A.同样,如果≠0(1=2,3?,n),也可以定义左乘A的屯?A=.(4)式中屯(￡El,,…,)t2E(zh,‟...,)均是行向量.向量中只有尾分量￡≠0.求屯向量的递推关系与(2)式相仿,此不多述.为了与屯向量对应,有时也将(3)式中的t向量记为.如果令()向量中第f(1≤f<)十分量以后的分量全为零,这种向量记为().此时相应的2()向量中只有二十相邻的分量(即尾分量)墨,而+l≠0,(‰.,.≠0)?由(2)式知,所对应的f向量,z向量除可相差一比饲常数外,是唯一确定的.隶t向量的计算量(抬乘除法的计算次数,下同)为3一.对某些对角占优的三对角阵t文[1]指出,t向量的分量按指数规律增大.更一般地讲.如果的元索为区间(--5,5)中均匀分布的随机数,数值试验表明.t向量的分量平均也按指数规律增大?当试验次数充分多时,t中后一十分量比相邻前一十分量平均增大刭1.57倍左右(绝对值)?这是由大量数据四则运算的统计规律决定的.因此在计算t向量时要防止机器溢出.t向量是一个有用的工具,它可以使算式和演绎大为简化.下面利用它来进行一些矩阵计算.\,,.,●一1lA量向维n为1997年6月高等学校计算数学2三对角线性方程组的解耦,系数阵的LR-1分解文[2]指出,对大型三对角方程组采用解耦算法可避免产生太的误差.故这里也给出一种解耦算法.设血=(5)式中为维未知向量,f为方程组的右端常数项向量,A非奇且次对角线上无零元?令矩阵R(厶)为n阶单位阵将第i列(行)换为f()后所得的阵.即尼=,厶=11ttz:…第行(6)令=厶AR-A与R相乘的结果,除第f列外,其余各元素均与A相同I 在第f列上,第i一1个分量变成了零,第f及f+1个分量变成了尾分量.同样,然后厶左秉的结果,使矩阵A在f行上第f一1个分量变为零.故=L4R=01c1第f行…第列(7)可见除第i行第i列外,其余元素均与A相同.这样便把A变换成具有z4”~t角块的三对角阵.为了使(7)式中新生成的6r,及与中的其它元素大小接近,可按比倒调整f向量使(6)式的及接近1-夸一;』,1-...编,,,.,“,_4;O,100唐达.三对角矩阵计算第2期P—LJ,则方程(5)可写成,=,.因此便把方程(5)解耦成二个独立的方程来求解.若A为对称阵,显然有￡=().故R=,所以变换后的矩阵A亦为对称阵. 以上是将(5)变换成二个独立的方程.同样,若要将方程(5)变换戚,,-+1个独立的方程,则可令AL1…L一.L_AR_一,…R.?(8)式中l<i<<…<‟一<‟≤.而().文[2]中给出的三对角方程组的解耦算法,需要附带解一个降了阶的,带宽为5的稀疏带状方程组.在(8)式中令=(一2,3…,n一1.n),则就变成对角阵D.这样便得到了A 的￡I1DRI1分解.在(10)的右端如果不乘厶并令AR=L,可知L是二对角的下三角阵.这样便得到A的LR分解;A一(AR)R-;LR～.(11)3求行列式,Jacobi矩阵的特征位与特征向量的关聂求detA时,不妨设c≠0(l,2‟.-?,n--1),否则可按~plac,定理展开成子行,日式来忤算.此时应先求出对应的t向量.令Q为单位阵将第一列换为t后的矩阵.印夸Q;l1●:B—AQ.(12)]●●●●1997年6月高等学校计算数学10l则B除第一列外,其余元素均与A相同.在B的第一列上,只有第个元素不为零,是尾分量.对(12)式两端求行列式.detB=detA?detQ.将detB按第一列展开,并注意到detQ=1,就得A之行列式的计算公式detA=(一1)+l?.Tic,.(13)1列递推关系求得:.;口0≠O),{,.一l=一b./a.?,.,(16)『I皇一(6.+l?t件l+a+l?tf1),dI+1(1≤i≤一2)与前述一样,令At=zJ,则2,向量中只有尾分量(第一个分量)≠O.在向量中如令第f(2≤:n)个分量以前的分时均为零.这种向量记为,,其对应的尾分量为.-.. 向量l”…及t”.所对应的尾分量分别为,及,令一..._.,作向量1O2唐选;三对角矩阵计算第2期则向量”…只有第i个分量不为零,它为面是距主对角线衰减的.如果A不是强主元矩阵,饲如A的元素为区间(一.5,5)中均匀分布的随机数,则如前所述,此时I肛l 的平均值约为1.57,即逆阵的元素距主对角线的衰减率平均为1.57.因此一般来讲,求A-1时,只要计算主对角线附近的元素,离对角线较远的元素之值将趋于零.当很大时,这就节省了很多运算量及存贮量.由(21)式可见,A_1的结构是很简单的.(因而也就是t及,)的性态,决定了A-1的性态?A_1的第行(列)及第f+1行(列)的元素,除an,嘶-f+l,nm口ff+】四个紧靠主对角线的元素外,其余元素均两两成固定比倒.这说明如果任意给出A的3n-2个元素,则A的个元素不能任意给出,它们间必须遵循上述的比倒规律.由(21)式计算A的个元时,求i;.及的四则运算量为o(n),它们相除(乘)的运算量为..故求一的全部元素,其乘除法运算量为+o(n),加减法运算量仅为0().文[4],[5]及其列出的参考文献中也曾讨论过三对角阵的逆阵的算法.文[4]算法的运算量为2n.+.(一).文[5]的算法,对于非对称阵,先按该文定理2求出逆阵的下三角部分元素,如果安排得好,其运算量为+o(n).然后按引理5并化俺后隶出逆阵的上三角部分元素,其运O一…一“,●●●●,●【一∞一}lgg?年6月高等学校计算教学103算量亦为+.(n).故总的运算量为2n+.(n).无论是文[4]或文[5],均未给出A-I的结构5误差分析,数值实验前面述及的种种算法,均建立在l向量的基础上.因此,这里主要讨论计算￡向量的浮点误差分析.由(z)式可得f一f查一?tl-z+一?一的位散即可.如果位致明显减少,则说明产生了数字的对消.如果在计算时t向时有一个分量产生严重对消,不一定使尾分量不精确I如果向量有相邻二个分量不精确,则一定导致尾分量的不精确.当发生这种情况时,可以类似于解方程组的迭代改善方法…来使有效位增加.本文提出的方程组解耦,求行列式及求逆矩阵等的算法,只是在l向量算出后再增加一些四则运算.可以证明它们都是向后稳定的算法.同理,文[5]求逆阵的算法也是向后稳定的.l04唐达:三对角矩阵计算第说明相对误差是小的..参考文献l萨乌里耶夫(袁兆鼎译),抛物型方程的冈格积分法,北京科学出版社.1963,l72一l73.2宋晓秋,袁兆鼎,刘穗贵,求解线性三对角方程组的解耦分解方法,系统工程与电子技术,12(199O),No.5,29—34.3蒋尔雄,对称矩阵计算,上海:上海科学技术出版社,1984,52—55.4陈增荣,任意三对角阵求逆,数值计算与计算机应用,8(1987),No.3,158--164.5Usmain,R.A.,Inverslon0fJaeobi‟sTKdiagomlmatrix,Comput.Math.App 1.,27(1994).No.8.59—66.6曹志浩等,矩阵计算和方程求根,第二版.北京;高等教育出版社,1987,30—46.。

毕业论文(设计) 题目三对角矩阵及其应用学生姓名吴莉学号20081314025院系数学与统计学院专业信息与计算科学指导教师杨兴东教授二Ｏ一二年五月八日目录O 引言 (1)1 三对角矩阵的行列式 (2)1.1 主要结论 (3)1.2 数值例子 (4)2 三对角矩阵的逆 (5)2.1 有关引理 (5)2.2 主要结论 (6)2.3 数值例子 (8)3 三对角矩阵的特征值及特征向量 (9)3.1 主要结论 (9)3.2 数值例子 (12)4 三对角矩阵的应用 (12)参考文献 (13)英文摘要 (15)致谢 (16)三对角矩阵及其应用吴莉南京信息工程大学数学与统计学院，南京210044摘 要 文章主要讨论了有关三对角矩阵的几个不同方面，首先讨论了一种特殊的三对角矩阵的行列式的计算方法，并且简要地给出了证明，接着通过几个引理和定理的证明给出了一种三对角矩阵的求逆算法，然后得到了计算一种三对角矩阵的特征值及其特征向量的公式，最后简要介绍了三对角矩阵法在粗甲醇精馏塔计算中的应用。

关键字 三对角矩阵；逆；特征值；特征向量;应用O 引言三对角矩阵在线性代数、计算数学、组合数学和应用数学上有着非常广泛的应用，在物理和工程方面也有着重要的应用，比如在数值分析中的差分方程、三次样条插值、及热传导方程的数值解等许多问题都涉及到了三对角矩阵，所以三对角矩阵的研究一直受到人们的关注。

本文主要介绍了一种特殊的三对角矩阵的行列式计算，三对角矩阵的求逆方法，三对角矩阵的特征值及其特征向量以及三对角矩阵在精馏中的应用，下文将一一给出说明。

1三对角矩阵的行列式1.1 主要结论为了便于计算，我们考虑如下的一种特殊三对角矩阵⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=αγβαγβαγβα n T ， 0≠βγ （1）定理1.1 设n T 为如（1）式所示的特殊三对角矩阵，则有（1）当042≠-βγαdet n T =βγαβγααβγαα42)4()4(211212-----++++n n n（2） 当042=-βγα时，det nn n n T α21+= 证明：令n n T D det =，按第一行展开便得到21---=n n n D D D βγα 设p ，q 为方程02=+-βγαx x 的两个解，即24,2422βγααβγαα--=-+=q p则（1）式可以写成)(211----=-n n n n pD D q pD D 。

三对角矩阵的快速算法三对角矩阵是一种特殊的矩阵，其非零元素只在主对角线和其相邻的两条对角线上。

由于其特殊性，我们可以采用一些特殊的算法来加速其求解。

下面就介绍三对角矩阵的快速算法。

1. Thomas算法Thomas算法是最常用的三对角矩阵求解算法。

该算法的基本思想是将原始矩阵转化为一个下三角矩阵和一个上三角矩阵的乘积形式，在此基础上使用迭代求解的方法。

具体实现过程如下：设三对角矩阵为A，向量b为待求解的线性方程组右端项，向量x为解向量，则原始的线性方程组可以表示为：Ax = b将A分解为下三角矩阵L和上三角矩阵U的乘积形式：则有：令y = Ux，可以得到：这是一个下三角矩阵和向量的乘积，可以使用前代法求解。

然后，再使用后代法求解Ux = y即可得到解向量x。

2. Sherman-Morrison算法设当前矩阵为A，其逆矩阵为A-1，将A作出如下形式的修正：A' = A + uvT其中u和v是向量，T表示向量的转置。

则有：A'-1 = A-1 - A-1uvTA-1 / (1 + vTA-1u)将上式表示为如下形式：显然，如果u和v分别为e1和en，即：u = [1, 0, …, 0]T这个式子就是Sherman-Morrison算法的核心推导。

然后，可以递归地使用这个式子进行修正即可。

3. Doolittle算法Doolittle算法是一种用于求解三对角矩阵的LU分解的算法，与一般的LU分解算法类似，其也将矩阵分解为下三角矩阵L和上三角矩阵U的乘积形式。

总结三对角矩阵是一种特殊的矩阵，其求解具有一些特殊性质，例如可以使用Thomas算法、Sherman-Morrison算法和Doolittle算法等算法进行快速求解。

由于其特殊性，这些算法都比一般的矩阵求解算法具有更高的效率，因此在实际应用中得到了广泛的应用。

矩阵奇异值分解算法的并行实现在现代数据分析和机器学习领域中，矩阵奇异值分解（Singular Value Decomposition, SVD）算法是一种重要且广泛应用的数学方法，用于特征提取、数据降维、矩阵压缩等任务。

然而，传统的SVD算法计算复杂度高，尤其是对于大规模矩阵，计算时间长，不适用于实时或高效处理。

因此，研究者们开始关注并行计算的技术，以提高SVD 算法的计算效率。

并行计算是将计算任务划分为多个子任务，同时进行计算以提高效率的一种方法。

在矩阵奇异值分解算法中，可以通过并行计算来加速计算过程，提高算法的效率。

下面将介绍一种基于并行计算的矩阵奇异值分解算法的实现方法。

首先，我们需要了解矩阵奇异值分解算法的基本原理。

矩阵奇异值分解是将一个矩阵分解为三个矩阵的乘积，即A = UΣV^T，其中A是待分解的矩阵，U和V是正交矩阵，Σ是对角矩阵，对角线上的元素称为奇异值。

传统的SVD算法需要对矩阵A进行特征值分解，计算复杂度较高。

而并行计算的矩阵奇异值分解算法可以将复杂的计算任务分解为多个子任务，同时进行计算，以提高计算效率。

一种常用的并行计算矩阵奇异值分解算法是基于分布式内存并行计算的。

在这种算法中，矩阵A被划分为多个小矩阵，每个进程或节点负责计算其中的一部分矩阵。

具体步骤如下：1. 将矩阵A分解为多个子矩阵：将矩阵A按行或列等分为多个子矩阵，每个子矩阵分配给一个进程或节点。

2. 并行计算SVD的近似值：每个进程或节点在分配到的子矩阵上计算SVD的近似值，得到局部的奇异值和奇异向量。

3. 收集局部的奇异值和奇异向量：各个进程或节点将局部计算得到的奇异值和奇异向量发送给主进程或节点，主进程或节点将它们进行合并。

4. 进一步迭代：基于合并后的局部奇异值和奇异向量，进行进一步的迭代计算，以获得更精确的SVD结果。

可以看到，通过并行计算的方法，我们可以实现对大规模矩阵的快速奇异值分解。

这种并行计算的矩阵奇异值分解算法已经在众多领域取得了广泛应用，例如图像处理、推荐系统、自然语言处理等。

大型实对称特征值问题的块Jacobi-Davidson方法的不精确求解谭静;汪晓红【摘要】The block Jacobi - Davidson method is effective for computing large scale real symmetric eigenvalue problems, the issues addressed being the multiple or clustered eigenpairs. The block Jacobi - Davidson method includes outer and inner iterative. The outer iterativc is used to compute the pairs of eigenvalues while the inner iterative is used for the correction equations. The more time - consuming computation lies in solving the correction equations. To handle the imprecise solution of the correction equation, we propose several block incomplete factorization methods to obtain the pre - conditioning matrix. Numerical experiments were also carried out to compare the effect of these methods.%块Jacobi—Davidson方法是计算大型实对称矩阵特征值问题的有效方法，可解决矩阵存在重特征值和密集特征值情况时的计算问题．块Jacboi—Davidson算法分为内外两层迭代，外层迭代计算矩阵特征对，内层迭代求解校正方程组，计算量主要花费是校正方程组的求解．针对校正方程的不精确求解，提出了几种构造预条件子的块不完全分解方法，并通过数值试验，对多种预条件子的效果进行比较．【期刊名称】《西安文理学院学报（自然科学版）》【年(卷),期】2012(015)004【总页数】7页(P38-44)【关键词】块Jacobi—Davidson方法;校正方程;不精确求解;预条件子块【作者】谭静;汪晓红【作者单位】南京航空航天大学理学院,南京210016 淮安信息职业技术学院基础部,江苏淮安223003;南京航空航天大学理学院,南京210016【正文语种】中文【中图分类】O241.6矩阵特征值问题Ax=λx，是数值代数的一个重要研究领域.无论是从理论还是应用角度，特征值问题的研究都具有重要的价值.对于分离较好的极端特征值问题，Krylov子空间迭代法，如Lanczos［1－2］方法和Arnoldi［3］方法，都有很好的收敛性.但当特征值的分离度不好，或要求内部特征值时，每步迭代中需要求解方程组［4］(A－σI)x=y.对于大型特征值问题，求大型方程组(A－σI)x=y的精确解有困难，或用迭代方法求近似解.但Lanczos方法和Arnoldi方法对于方程组不精确求解较敏感，影响了计算的收敛性.1996年，Sleijpen和Van der Vorst将Jacobi方法的校正思想和Davidson方法的内外迭代格式融合起来，提出了Jacboi－David son方法［5］.Jacboi－Davidson方法是求解大型特征值问题的有效迭代方法，特别是对求解中间特征值有优势，其核心是在近似特征向量的正交补上对其进行修正，通过求解校正方程得到修正向量，然后利用修正向量扩充搜索子空间.良好的收敛性质使得我们只需适度地近似求解校正方程，用少量的迭代保持良好的收敛速度.针对所要求的特征值问题是重特征值或者是密集特征值时，克服特征值的分离度差造成问题的性态可能会变坏的困难，有学者提出了块Jacboi－Davidson方法，它可以同时计算多个特征对.算法1 块Jacobi－Davidson方法(1)输入矩阵A，投影子空间最大维数m，块大小l，选取n×l的列正交规范矩阵V1;(2)对k=1，2，…，m.(i)计算(ii)求Hk的l个最大特征对(iii)计算Ritz向量(iv)计算残量，若满足收敛条件则停止;否则进行下一步;(v)解校正方程组.求新的方向向量矩阵Tk=［t(k)1，…，t(k)l］，其中块Jacboi－Davidson算法分为内外两层迭代，外层迭代计算矩阵特征对，内层迭代求解校正方程组，计算量主要花费是校正方程组的求解.本文主要研究块Jacboi －Davidson方法中校正方程的不精确求解方法.下面首先研究不精确求解校正方程的增广形式的预处理解法.针对求解校正方程近似解过程中方程的系数矩阵条件数逐渐变坏问题，然后研究求解校正方程的预处理方法，构造校正方程的预条件子，给出了构建预条件子块不完全LU分解方法.最后通过数值实验分析、研究不精确求解校正方程时算法的有效性，并尝试对多种预条件子的效果进行比较.1 校正方程的不精确求解在块Jacobi－Davidson方法中，每个迭代步需要求解如下形式的校正方程:其中，ΦTΦ=I.块Jacobi－Davidson方法的工作量主要在于对校正方程(1)的求解.对于校正方程，除了可以精确求解以外，更多的时候使用不精确求解.随着逼近程度的提高，方程的系数矩阵的条件数会变坏，我们可以采用预处理克服问题.设Ki是A－λiI的左预条件子.将预条件子K限制到Φ的垂直补空间，即用作为校正方程的预处理算子.下面给出一种不精确求解校正方程的方法——增广形式的预处理解法.定理1［7］若记其中Fi=(I－ΦΦT)(A－λiI)(I－ΦΦT)，i=1，2，…，l且ΦTΦ=I，则(1)式等价于即ti是(2)式的解当且仅当ti是(3)式的解.其中fi=－ΦT(A－λiI)ti.定理2［7］对于Z∈Rn×p和Y∈Rn×p，有其中V≡ZTY，当且仅当ZTY非奇异时可逆.当Ki是(A－λiI)易于求逆的近似矩阵时.因为则又由式(4)得所以于是，将(6)式代入(5)式，从而可得:且因此，代入式(3)，可以得到左预处理的增广校正方程如下:那么，在块Jacobi－Davidson方法的实现中，求解校正方程就可以转化为求解上述(9)式.而对于(9)式的求解，可以运用Krylov子空间方法，如CG方法，GMRES 方法等.2 预条件子的生成不精确求解校正方程，预优化必不可少.因为坏条件数正是在所要求的特征方向上产生的，故预条件子可以对特征向量的计算起促进作用.因此，构造适合的预条件子即构造A－λiI的近似矩阵Ki，是不精确求解校正方程的重要环节.构造预条件子常用的方法有很多，如块Jacobi方法、不完全LU分解、不完全Cholesly分解等.2.1 块不完全LU分解对于大型的问题，直接使用不完全LU分解时一次参与运算的数据很大.为了提高效率，我们考虑将矩阵进行分块，通过计算分块子矩阵来提高运算效率.下面介绍LU分解分块算法［8］.首先，将矩阵A分为四块，则可将LU分解写为:其中，A11为d1×d1阶矩阵，A22为(n－d1)×(n－d1)阶矩阵.由(10)式有:算法2 分块不完全LU分解算法(1)对式(11)中的分块矩阵A11进行ILU分解，可得L11，U11;(2)求单位下三角阵L11=(lij)d1×d1的逆矩阵B1=(bij)d1×d1.B1仍为单位下三角阵即主对角元素均为1，且(3)代入式(12)得，U12=L－111A12=B1A12;(4)求上三角阵U11=(uij)d1×d1的逆矩阵C1=(cij)d1×d1.C1仍为上三角阵且(5)代入式(13)得(6)改写式(14)，设，则式(14)成为的LU分解，即算法中对的LU分解，可继续对其进行分块计算，如此一直分下去.只要每次选取适当的分块大小，分块后的各子矩阵的相应运算就有可能在计算机的高缓中进行，从而大幅提高运算速度.算法中，求逆矩阵我们并没有使用高斯约旦消元法，而是充分利用三角矩阵的特点，省去了不必要的计算.2.2 三对角块不完全LU分解对于一些实际问题，其偏微分方程离散后得到的线性方程组系数矩阵A是如下形式的块三对角矩阵.其中，Ti为ni阶矩阵，Fi和Ei为对角矩阵.将A分解为，其中由此得到当Ti为带状矩阵时，可以使用一种块不完全分解法［9］，即用窄带宽矩阵近似带状矩阵的逆.我们用(18)式作为式(17)的近似:其中，(B)(p)表示一个上下半带宽均为p的近似于B的带状矩阵.算法3 三对角块不完全LU分解(1)D1=T1(2)对i=2，3，…，m(i)计算(ii)取Ci－1为中上下半带宽均为p的带状部分;(iii)计算(iv)计算Di=Ti－EiCi－1Fi;(v)舍去Di中以p为上下半带宽以外的部分.不完全分解因子L和U分别由(16)式中，将和替换成Di和Ui得到.3 数值试验本算法由MATLAB2009a实现.A为n阶实对称矩阵，初始矩阵V1为随机产生的列正交规范阵.不完全LU分解法记为ILU，不完全Cholesly分解记为IC，块不完全LU分解记为BILU.对于校正方程组的求解使用增广形式的预处理解法.A=其中B=为20阶矩阵，I为20阶单位阵，取δ=0.6.首先，利用Matlab软件包中自带的程序计算A的4个最大特征值如下:当n=500时，λ1=6.775 6，λ2=6.775 4，λ3=6.732 1，λ4=6.731 5当n=1 000时，λ1=6.786 4，λ2=6.786 3，λ3=6.775 0，λ4=6.774 9当n=2 000时，λ1=6.789 2，λ2=6.789 2，λ3=6.786 3，λ4=6.786 3如果计算A－λI的预条件子，那么每次外迭代都要计算l个预条件子，工作量较大.于是，我们尝试计算矩阵A的预条件子，而不是整个A－λI.这样，预条件子只在一开始计算，不考虑λ的值.分别用对A－λI和A构造预条件子(取对角)的算法计算500阶和1 000阶两种情况，结果见表1.由表1的数值结果可以看出，直接对矩阵A构造预条件子，是行的通的，而且效果很好.下面的计算中将直接对矩阵A构造预条件子.当分别使用对角、IC、ILU和BILU预条件子，计算结果见表2、表3和表4.表1 分别对A－λI和A构造预条件子的算法计算结果(n=500和n=1 000)构造预条件子的对象 n 特征值残量范数迭代次数 CPU时间/s 500 6.775 6 6.775 4 6.723 1 6.731 5 2.899 3e－007 5.739 2e－007 5.287 0e－006 1.440 5e－004 13 11.232 686 A－λI 1 000 6.786 4 6.786 3 6.775 0 6.774 9 1.099 7e－005 6.666 2e－006 9.791 0e－005 1.997 2e－004 16 91.500 589 50013 11.169 595 A 6.775 6 6.775 4 6.723 1 6.731 5 1.395 1e－006 5.779 9e－006 1.340 5e－004 4.609 6e－005 1 000 6.786 4 6.786 3 6.775 0 6.774 9 6.321 1e－005 5.195 3e－004 7.832 2e－004 2.694 4e－004 15 84.018 093表2 各种预条件子下的迭代次数注:*表示不收敛n迭代次数对角IC ILU BILU 500 13 20 21 19 1 000 15 26 24 22 2 000 20 *34 36表3 各种预条件子下的CPU时间注:*表示不收敛n时间对角CPU IC ILU BILU 500 11.169 595 71.768 032 15.197 446 15.595 842 1 000 84.018 093 684.914 107 93.418 466 109.926 004 2 000 881.568 323 *964.020 571 1 397.875 919表4 各种预条件子下的残量范数注:*表示不收敛n特征值残量范数对角IC ILU BILU 500 6.775 6 6.775 4 6.723 1 6.731 5 1.395 1e－006 5.779 9e－006 1.340 5e－004 4.609 6e－005 3.253 6e－005 1.400 0e－005 6.892 2e－0043.549 6e－004 2.702 2e－0074.446 4e－006 4.149 4e－004 1.191 1e－0044.123 9e－006 1.641 7e－0055.612 9e－005 7.340 7e－004 1 0006.786 4 6.786 3 6.775 0 6.7749 6.321 1e－005 5.195 3e－0047.832 2e－004 2.694 4e－004 1.226 6e－004 2.364 7e－004 9.465 2e－0058.214 9e－005 1.523 1e－004 7.846 6e－005 1.604 7e－004 5.978 2e－0049.984 1e－005 1.016 6e－004 1.713 8e－004 8.704 2e－004 2 000* 1.331 8e－004 3.732 0e－004 8.980 7e－004 2.728 9e－005 6.789 2 6.789 2 6.786 3 6.786 3 8.229 1e－005 8.931 5e－005 3.113 2e－004 6.994 2e－005 8.990 2e－005 1.979 1e－004 2.431 6e－004 6.146 1e－004从表2、表3、表4中的数据可以看出，块Jacobi－Davidson方法的效果，与预条件子的选取有较大的关系.此外，结果显示，越高效的预条件子并不必然提高块Jacobi－Davidson方法的收敛性.下面以500阶矩阵为例，作出各种预条件子收敛情况的图像见图1.从图1中可以明显看出，使用对角预条件子，块Jacobi－Davidson方法收敛的最快.图1 各种预条件子下的收敛情况(n=500)4 结论对于块Jacobi－Davidson方法的校正方程不精确求解时，可以直接计算矩阵A 的预条件子，而不是整个A－λI.此外，块Jacobi－Davidson方法的运算效果，与预条件子的选取有较大的关系，并且越高效的预条件子并不必然提高块Jacobi－Davidson方法的收敛性.因此在实际计算块Jacobi－Davidson方法时，可以选取对角预条件子、ILU预条件子等较简单的预处理子.［参考文献］［1］ LANCZOS C.An iteration method for the solution of the eigenvalue problem of linear differential and integral operators［J］.J.Res.Nat.Bur.Standard，1950，45:255－282.［2］ PAIGE C putational variants of the Lanczos method of the eigenproblem［J］.J.Inst.Math.Appl，1972，10:373－381.［3］ SAAD Y.Variations on Arnoldi's method for computing eigenelements of large unsymmetric matrices［J］.Linear Algebra Appl，1980，34:269－295.［4］ HEINRICH V.A new justification of the Jacobi－Davidson method forlarge eigenproblems［J］.Linear Algebra and its Applications，2007，424:448－455.［5］ SLEIJPEN G L G，VAN Der Vorst H A.A Jacobi－Davidson method for linear eigenvalue problems［J］.SIMA J Matrix Anal Appl，1996，17:401－425.［6］王岩青.实对称矩阵特征值问题的迭代块Jacobi－Davidson方法［J］.解放军理工大学学报:自然科学版，2002，3 (1):93－96.［7］康艳艳.求解大型对称特征值问题的改进块Jacobi－Davidson方法［J］.西安文理学院学报:自然科学版，2010，13 (2):44－49.［8］纪坤，陈建平，石振国，等.矩阵三角分解分块算法的研究与实现［J］.计算机应用与软件，2010，27(9):72－74.［9］李晓梅，吴建平.稀疏线性方程组不完全分解预条件方法［J］.计算机工程与科学，2006，28(8):59－62.。