并行计算-矩阵特征值计算 -
矩阵特征值计算
矩阵特征值计算
1 矩阵特征值计算
矩阵特征值计算是数学分析中一种重要的概念，它与矩阵的属性
密切相关。
一般情况下，矩阵特征值是指与矩阵有关的实值函数。
它
可以提供有关矩阵行列式和特征向量的有关信息。
矩阵特征值计算的定义为：设A是m×n的实矩阵，若存在一个实
数λ，使得A的每个元素都等于λ的倍数，则称λ为A的特征值（eigenvalue）。
所有不同的特征值构成A的特征值多项式，即特征
多项式的多项式，由
λ ^ n+λ ^ (n-1)+...+λ ^ 1+λ ^ 0
组成。
这里λ重复出现n次，证明A有n个线性无关的特征向量。
特征值是矩阵分析中最重要的因素。
一个矩阵的特征值有助于分
析数学问题，特别是矩阵和相对应的特征向量，以及A的分解及A的秩。
特征值及其特征向量也有助于我们解决具体问题，如求解最大值、最小值、最优解等。
计算矩阵特征值一般有两种方法：一种是近似方法，即迭代法；
另一种是精确方法，即行列式法和特征多项式法。
迭代法可以快速求
得与原矩阵几乎相等的特征值，但精确的特征值只能通过行列式法和
特征多项式法求出。
由于矩阵特征值的重要性，众多学者和研究者致力于提出更好的特征值计算方法。
而在数学和计算机科学领域，研究者更多地关注矩阵特征值，研究其计算方法和性能，以求取更好的计算效率提高计算精度。
矩阵特征值求解
矩阵特征值求解的分值算法12组1. 1矩阵计算的基本问题(1) 求解线性方程组的问题.即给定一个n 阶非奇异矩阵A 和n 维向量b ,求 一个n 维向量X,使得Ax =b (1.1.1 )(2) 线性最小二乘问题，即给定一个mx n 阶矩阵A 和m 维向量b ,求一个n 维向量使得 |A X -b | =min{ |Ay -比严 R n }(3) 矩阵的特征问题，即给定一个n 阶实(复)矩阵A ,求它的部分或全部特 征值以及对应的特征向量，也就是求解方程Ax = Z xA 的属于特征值A 的特征向量。
在工程上,矩阵的特征值具有广泛的应用，如大型桥梁或建筑物的振动问题： 机械和机件的振动问题；飞机机翼的颤振问题 ；无线电电子学及光学系统的电磁 振动问题;调节系统的自振问题以及声学和超声学系统的振动问题 .又如天文、地 震、信息系统、经济学中的一些问题都与矩阵的特征值问题密切相关。
在科学上，计算流体力学、统计计算、量子力学、化学工程和网络排队的马 尔可夫链模拟等实际问题，最后也都要归结为矩阵的特征值问题.由于特征值问 题在许多科学和工程领域中具有广泛的应用，因此对矩阵的特征值问题的求解理 论研究算法的开发软件的制作等是当今计算数学和科学与工程计算研究领域的 重大课题,国际上这方面的研究工作十分活跃。
1.2矩阵的特征值问题研究现状及算法概述对一个nxn 阶实(复)矩阵A,它的特征值问题，即求方程(I.1.3)式的非平凡 解,是数值线性代数的一个中心问题.这一问题的内在非线性给计算特征值带来 许多计算问题.为了求(1.1.3)式中的A , —个简单的想法就是显式地求解特征方 程det(A —几I) = 0除非对于个别的特殊矩阵，由于特征方程的系数不能够用稳定的数值方法由 行列式的计算来求得，既使能精确计算出特征方程的系数，在有限精度下，其特征 多项式f ") =det(A-ZJ)的根可能对多项式的系数非常敏感 能在理论上是有意义的，实际计算中对一般矩阵是不可行的 数较大，则行列式det(A -几I)的计算量将非常大；其次，根据 数大于四的多项式求根不存在一种通用的方法 ，基于上述原因，人们只能寻求其 它途径.因此,如何有效地！精确地求解矩阵特征值问题，就成为数值线性代数领 域的一个中心问题.目前,求解矩阵特征值问题的方法有两大类：一类称为变换方法，另一类称为 向X,(1.1.2 )(1.1.3 ) 一对解(4 X),其中R(C),x- R n (C n ),即A 为矩阵A 的特征值，X 为矩阵(121 ).因此，这个方法只 .首先,若矩阵A 的阶 Galois 理论,对于次量迭代方法.变换方法是直接对原矩阵进行处理，通过一系列相似变换，使之变换成一个易于求解特征值的形式，如Jacobi算法,Givens算法,QR算法等。
基于GPU的多项式矩阵特征值并行计算
基于GPU的多项式矩阵特征值并行计算
杨智诚;黄友钦;傅继阳
【期刊名称】《中国科技论文》
【年(卷),期】2016(011)008
【摘要】以多项式矩阵特征值求解问题为背景,推导任意项式特征值计算方法,阐述了基于图形处理单元（GPU）的并行计算机理并建立Matlab快速并行计算平台。
以二次特征值计算为例,分别测试基于GPU的单精度类型计算与双精度类型计算的效率与误差,同时验证公式的正确性。
数值计算表明,单精度计算与双精度计算能够大幅度降低计算时间,提高计算效率,误差对比分析也表明GPU的计算结果能够满足一定的工程计算要求。
【总页数】4页(P957-960)
【作者】杨智诚;黄友钦;傅继阳
【作者单位】广州大学工程结构抗风与结构安全国际科技合作基地,广州510006【正文语种】中文
【中图分类】TU311.1
【相关文献】
1.大型矩阵特征值问题并行计算研究概况 [J], 李晓梅;罗晓广
2.一类特殊矩阵特征值的并行计算 [J], 尹小红
3.在PVM上并行计算实对称三对角矩阵特征值 [J], 罗晓广;李晓梅
4.基于 GPU 的多项式矩阵特征值并行计算 [J], 杨智诚;黄友钦;傅继阳
5.CPU-OpenMP和GPU-CUDA并行计算技术对矩阵乘法运算的加速效果分析[J], 张岩
因版权原因，仅展示原文概要，查看原文内容请购买。
并行算法在矩阵计算中的应用研究
２并行计算机理论基础
（２）阵列处理机阵列处理机通常系由一个控制单元ＣＵ（Ｃｏｎｔｒｏｌ Ｕｎｉｔ）、ｎ个处理单 元ＰＥ（Ｐｒｏｃｅｓｓｉｎｇ Ｅｌｅｍｅｎｔ）、Ｍ个存储模块Ｍ和一个互连网络１Ｎ（Ｉｎｔｅｒｃｏｒｍｅｃｔｉｏｎ Ｎｅｔｗｏｒｋ） 所组成。ＣＵ将单一指令（单指令流）播送至各个ＰＥ，而所有活动的ＰＥ将从相应的Ｍ中 取出各自所需的数据元素（多数据流），以同步的方式执行该条指令。在阵列机中，ＩＮ也 常叫做对准（Ａｌｉｇｎｍｅｍ）或置换（Ｐｅｒｍｕｔａｔｉｏｎ）网络，它用来提供各ＰＥ之间或ＰＥ与Ｍ之 间的通信连接。
ｓｔｕｄｉｅｄ ｔｈｅ ｍａｔｒｉｘ ｍｕｌｔｉｐＨｃｍｉｏｎ ｐａｒａｌｌｅｌ ｐｒｏｃｅｓｓｉｎｇ ｔｅｃｈｎｉｑｕｅｓ ｉｎ ｄｉｓｔｒｉｂｕｔｅｄ
ｅｎｖｉｒｏｎｍｅｎｔ．Ｗｅ ａｎａｌｙｚｅｄ ｔｈｅ ｐａｒａｌｌｅｉ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｍａｔｒｉｘ ｍｕｌｔｉｔ）ｌｉｃａｔｉｏｎ ｉｎ ＭＩＭＤ，
１绪论
用。主要分为三个大的部分。 第一部分：对并行处理机和并行处理技术进行了概述。介绍了并行算法的基本概念和
相关的评价指标。 第二部分：研究了矩阵乘法在分步式处理机上的并行处理技术。在原有Ｓｔｒａｓｓｅｎ算法
和Ｗｉｎｏｇｒａｄ算法的基础上，提出了对ＳＩＭＤ．Ｌｃ结构的并行机，可适用的按行分块的并行 处理算法，给出了算法的实现步骤，并且分析了该算法的加速比。对ＭＩＭＤ上的矩阵乘法 的并行算法进行了分析，进一步分析了算法的可扩展性和并行处理的效率。
２．２．１ ＳＩＳＤ型计算机
目前所设计和建造的大多数计算机都是这类机器，它不属于并行机，而属于传统的 串行计算机。它有两种非常重要的特征：串行的和确定的。
２．２．２ ＳＩＭＤ型并行机
矩阵特征值分解算法的并行实现
矩阵特征值分解算法的并行实现特征值分解（Eigenvalue Decomposition）是线性代数中的一个重要概念，广泛应用于数学、物理、工程等领域。
以矩阵的特征值和特征向量为基础，可以提取出矩阵的结构信息，对于解决线性系统、优化问题、信号处理等具有重要意义。
然而，传统的特征值分解算法在处理大规模矩阵时面临计算量巨大、执行时间长的问题。
为了加速特征值分解的计算过程，研究者们逐渐引入并行计算的思想，对矩阵特征值分解算法进行并行实现。
一、并行计算的概述并行计算是指将一个复杂的计算任务划分成多个小任务，通过多个处理单元同时进行计算，从而提高计算效率的一种方法。
在特征值分解算法中，通过合理划分矩阵的计算任务，并行计算可以有效降低计算时间。
二、传统特征值分解算法传统特征值分解算法主要有幂法、QR方法等。
幂法通过迭代矩阵的特征向量来逼近最大特征值，并进一步得到特征向量。
QR方法则通过将矩阵转化为上Hessenberg矩阵，然后通过反复QR分解来逼近矩阵的特征值。
这些算法都存在计算复杂、耗时长的问题，在处理大规模矩阵时效率较低。
三、并行特征值分解算法为了提高特征值分解算法的计算效率，研究者们开始尝试将并行计算技术引入到特征值分解算法中。
一种常用的并行特征值分解方法是块特征值分解法（Block Eigenvalue Decomposition）。
该方法将矩阵划分为多个较小的块，每个块由一个处理单元处理。
通过并行计算每个块的特征值和特征向量，最后组合得到整个矩阵的特征值和特征向量。
块特征值分解法可以显著提高计算效率，适用于处理大规模矩阵。
四、算法实现与优化在并行特征值分解算法的实现过程中，需要考虑如何合理划分矩阵的计算任务，确保各个处理单元之间的负载均衡。
通常可以将矩阵按行或者按列等方式划分成多个子矩阵，并行计算各个子矩阵的特征值和特征向量。
此外，还可以采用多级划分的方法，将矩阵分解成更小的块，进一步提高并行计算效率。
第章矩阵特征值的计算
第章矩阵特征值的计算矩阵特征值是矩阵理论中的一个重要概念，它在很多领域中都有广泛的应用，如物理、化学、工程等。
本文将从特征值的定义、计算方法和应用举例等方面进行阐述。
一、特征值的定义对于一个n阶方阵A，如果存在一个非零向量x，使得Ax=kx，其中k 是一个常数，那么k称为A的特征值，x称为A的对应于特征值k的特征向量。
从定义可以看出，矩阵A的特征值和特征向量是成对出现的，特征向量可以是一个实数或是一个向量，特征值可以是实数或是复数。
二、特征值的计算方法1.直接计算法此方法适合于较小的矩阵。
给定一个n阶矩阵A，首先构造特征方程det(A-λI)=0，其中I是n阶单位矩阵，λ是未知数，然后求解特征方程得到特征值，将特征值代入(A-λI)x=0求解对应的特征向量。
2.幂法幂法是一种迭代方法，适用于大型矩阵。
假设特征值的绝对值最大，那么从非零向量b开始迭代过程，令x0=b，求解x1=A*x0，然后再将x1作为初始值，求解x2=A*x1，以此类推，直到收敛为止。
最后，取最终得到的向量xn，其模即为特征值的近似值。
3.QR方法QR方法是一种迭代方法，可以用于寻找特征值和特征向量。
首先将矩阵A分解为QR，其中Q是正交矩阵，R是上三角矩阵，然后对R进行迭代，重复进行QR分解，直到收敛。
最后，得到的上三角矩阵的对角元素即为特征值的近似值，在QR分解的过程中，特征向量也可以得到。
三、特征值的应用举例1.物理学中的量子力学量子力学中的哈密顿算符可以表示为一个矩阵，物理量的测量值就是对应的特征值。
例如，电子的自旋可以有上自旋和下自旋两种状态，上自旋对应的特征值为1，下自旋对应的特征值为-12.工程中的振动问题在工程中，矩阵特征值可以用来求解振动问题。
例如，振动系统的自由度决定了特征向量的个数，而特征值则表示了振动的频率。
通过计算矩阵的特征值和特征向量，可以预测系统的振动频率和振型。
3.网络分析中的中心性度量在网络分析中，矩阵特征值可以用来计算节点的中心性度量。
矩阵特征值计算矩阵的特征值和特征向量
矩阵特征值计算矩阵的特征值和特征向量矩阵是线性代数中的重要概念之一，它在众多学科领域中都有广泛的应用。
而矩阵的特征值和特征向量则是矩阵分析与应用中的核心内容之一。
本文将详细介绍矩阵特征值的计算方法，以及如何求解矩阵的特征向量。
1. 特征值和特征向量的定义首先，我们来了解一下什么是矩阵的特征值和特征向量。
给定一个n阶方阵A，如果存在一个数λ以及一个非零n维列向量X，使得满足下述条件：AX = λX那么，λ就是矩阵A的一个特征值，而X则是对应于特征值λ的特征向量。
特征值和特征向量的求解在很多应用中都具有重要的意义。
2. 特征值的计算方法接下来，我们介绍几种常见的特征值计算方法。
2.1 特征多项式法特征多项式法是求解特征值的一种常用方法。
它利用方阵A减去λ乘以单位矩阵I之后的行列式为零的性质，构造出特征多项式，并求解多项式的根即可得到特征值。
举个例子，对于二阶方阵A = [a, b; c, d]，其特征多项式为：| A - λI | = | a-λ, b; c, d-λ | = (a-λ)(d-λ) - bc = 0解这个方程可以得到A的特征值。
2.2 幂迭代法幂迭代法也是一种常见的特征值计算方法。
它利用特征向量的性质，通过迭代计算来逼近矩阵的特征值。
其基本思想是，给定一个初始向量X0，不断迭代计算：Xk+1 = AXk然后对得到的向量序列进行归一化处理，直到收敛为止。
最后得到的向量X就是对应的特征向量，而特征值可以通过如下公式计算：λ = X^TAX / X^TX2.3 QR方法QR方法是一种数值稳定性较好的特征值计算方法。
它利用矩阵的QR分解的性质来逐步逼近矩阵的特征值。
首先，对矩阵A进行QR分解，得到一个正交矩阵Q和一个上三角矩阵R。
然后，将分解后的矩阵R与矩阵Q逆序相乘，得到一个新的矩阵A'。
重复进行QR分解和相乘的操作，直到收敛为止。
最后，得到的矩阵A'的对角线上的元素即为矩阵A的特征值。
超级计算技术中的并行算法与矩阵运算
超级计算技术中的并行算法与矩阵运算在现代科学和工程领域中，超级计算技术发挥着至关重要的作用。
为了解决复杂问题，超级计算机采用了并行算法和矩阵运算等技术，以实现高效的计算和分析。
本文将探讨超级计算技术中的并行算法和矩阵运算，并分析其应用与发展趋势。
首先，我们来了解一下超级计算技术中的并行算法。
并行算法是指将复杂的计算任务分解成多个子任务，然后并发地运行于多个计算单元上，以提高计算效率和性能。
并行算法的设计需要考虑任务分解、通信和同步等关键问题。
常用的并行算法包括分治法、并行排序和并行搜索等。
对于矩阵运算来说，超级计算技术也发挥着重要的作用。
矩阵运算是指对矩阵进行基本的数学运算，诸如加法、减法、乘法和求逆等。
这些矩阵运算在科学计算、图像处理和人工智能等领域中广泛应用。
超级计算技术能够利用并行算法加速矩阵运算的速度，提高计算效率。
在超级计算技术中，矩阵乘法是一种常见且重要的矩阵运算。
矩阵乘法的基本思想是将两个矩阵相乘，得到一个新的矩阵。
然而，传统的矩阵乘法算法在大规模矩阵计算时效率较低，因此需要并行算法的支持。
并行矩阵乘法算法采用了多种策略，如按块分配、循环分配和网络划分等，以充分利用计算资源。
除了矩阵乘法，超级计算技术还可以应用于其他矩阵运算，如矩阵分解和特征值计算等。
矩阵分解是将一个矩阵分解成多个部分矩阵的过程，常见的矩阵分解包括LU分解、QR分解和SVD分解。
这些分解可以帮助我们更好地理解和处理复杂的数学问题。
而特征值计算则是计算一个矩阵的特征值和特征向量，对于解决线性方程组和优化问题都具有重要意义。
随着科学技术的不断发展，超级计算技术中的并行算法和矩阵运算也在不断演进。
首先，随着计算单元和存储器的不断增加，我们可以使用更大规模的矩阵进行计算，从而解决更加复杂的科学问题。
其次，新的并行算法和矩阵运算技术的提出，使得超级计算机能够更快速地处理大规模数据，实现更高效的计算。
最后，超级计算技术的发展也推动了云计算和人工智能等领域的进步，使得计算资源得到充分利用。
矩阵特征值与特征向量的计算与应用
添加标题
矩阵特征值与特征向量的计算：Octave提供了计算矩阵特征值和特征向量的函数，如 `eig()`和`svd()`等，可以方便地进行矩阵分解和特征值计算。
添加标题
数值计算环境：Octave具有友好的用户界面和交互式命令行，支持脚本编写和函数封 装，方便用户进行数值计算和分析。
添加标题
应用领域：Octave广泛应用于科学计算、数据分析、控制系统等领域，尤其在处理大 规模矩阵运算和数值分析方面具有优势。
雅可比法
定义：雅可比法是一种求解矩阵特征值和特征向量的方法，通过迭代的方式逐步逼近矩阵的特 征向量。
适用范围：适用于实对称矩阵和非实对称矩阵。
计算步骤：通过迭代公式逐步计算特征向量，直到收敛。
优缺点：雅可比法具有简单易行、收敛速度快等优点，但需要选择合适的初始向量和迭代参数， 否则可能无法收敛或收敛到非特征向量。
矩阵特征值与特征向 量的计算与应用
XX,a click to unlimited possibilities
汇报人：XX
目录 /目录
01
点击此处添加 目录标题
04
矩阵特征值与 特征向量的计 算方法
02
矩阵特征值与 特征向量的基 本概念
05
矩阵特征值与 特征向量的应 用实例
03
矩阵特征值与 特征向量的应 用场景
06
矩阵特征值与 特征向量的计 算工具与软件
01 添加章节标题
02
矩阵特征值与特征向量 的基本概念
特征值与特征向量的定义
特征值：矩阵A中与单 位向量相乘后得到一个 与原向量共线的向量
特征向量：矩阵A中与 特征值对应的向量
特征值与特征向量的性质
特征值和特征向量是矩阵的重要属性，它们描述了矩阵对向量空间的作用。
矩阵特征值计算
第五章矩阵特征值计算与线性方程组的求解问题一样，矩阵特征值与特征向量的计算也是数值线性代数的重要内容. 在理论上，矩阵的特征值是特征多项式方程的根，因此特征值的计算可转化为单个多项式方程的求解. 然而对于高阶矩阵，这种转化并不能使问题得到简化，而且在实际应用中还会引入严重的数值误差. 因此，正如第二章指出的，我们一般将多项式方程求解转化为矩阵特征值计算问题，而不是反过来.本章介绍有关矩阵特征值计算问题的基本理论和算法. 与非线性方程求根问题类似，计算矩阵特征值的算法也是迭代方法①.5.1基本概念与特征值分布本节先介绍矩阵特征值、特征向量的基本概念和性质，然后讨论对特征值分布范围的简单估计方法.5.1.1基本概念与性质定义5.1：矩阵A=(a kj)∈ℂn×n，(1) 称φ(λ)=det(λI−A)=λn+c1λn−1+⋯+c n−1λ+c n为A的特征多项式(characteristic polynomial)；n次代数方程φ(λ)=0为A的特征方程(characteristic equation)，它的n个根：λ1,⋯,λn，被称为A的特征值（eigenvalue）. 此外，常用λ(A)表示A的全体特征值的集合，也称为特征值谱(spectrum of eigenvalue).(2) 对于矩阵A的一个给定特征值λ，相应的齐次线性方程组(λI−A)x=0 , (5.1)有非零解（因为系数矩阵奇异），其解向量x称为矩阵A对应于λ的特征向量（eigenvector）.根据方程(5.1)，我们得出矩阵特征值与特征向量的关系，即Ax=λx .(5.2)第三章的定义3.5就利用公式(5.2)对矩阵特征值和特征向量进行了定义，它与定义5.1是等价的. 另外，同一个特征值对应的特征向量一定不唯一，它们构成线性子空间，称为特征子空间(eigenspace).我们一般讨论实矩阵的特征值问题. 应注意，实矩阵的特征值和特征向量不一定是实数和实向量，但实特征值一定对应于实特征向量（方程(5.1)的解），而一般的复特征值对应的特征向量一定不是实向量. 此外，若特征值不是实数, 则其复共轭也一定是特征值（由于特征方程为实系数方程）. 定理3.3表明，实对称矩阵A∈ℝn×n的特征值均为实数，存在n个线性无关、且正交的实特征向量，即存在由特征值组成的对角阵Λ和特征向量组成的正交阵Q，使得:A=QΛQ T.(5.3)例5.1（弹簧－质点系统）：考虑图5-1的弹簧－质点系统，其中包括三个质量分别为m1、m2、m3的物体，由三个弹性系数分别为k1,k2,k3的弹簧相连，三个物体的位置均为时间的函数，①如果用有限次运算能求得一般矩阵的特征值，则多项式方程求根问题也可用有限次运算解决，这与阿贝尔证明的“高于4次的多项式并不都有用初等运算表示的求根公式”的理论矛盾.这里考查三个物体偏离平衡位置的位移，分别记为y 1(t), y 2(t), y 3(t). 因为物体在平衡状态所受的重力已经和弹簧伸长的弹力平衡，所以物体的加速度只和偏离平衡位置引起的弹簧伸长相关. 根据牛顿第二定律以及胡克定律（即弹簧的弹力与拉伸长度成正比）可列出如下微分方程组②： My ′′(t)+Ky(t)=0 ,其中y (t )=[y 1(t)y 2(t)y 3(t)]T ,M =[m 1000m 2000m 3],K =[k 1+k 2−k 20−k 2k 2+k 3−k 30−k 3k 3] . 在一般情况下，这个系统会以自然频率ω做谐波振动，而y 的通解包含如下的分量： y j (t )=x j e iωt ,(j =1,2,3)其中i =√−1，根据它可求解出振动的频率ω及振幅x j . 由这个式子可得出：y j ′′(t )=−ω2x j e iωt ,(j =1,2,3)代入微分方程，可得代数方程:−ω2Mx +Kx =0,或Ax =λx ,其中A =M −1K ，λ=ω2. 通过求解矩阵A 的特征值便可求出这个弹簧－质点系统的自然频率（有多个）. 再结合初始条件可确定这三个位移函数，它们可能按某个自然频率振动（简正振动），也可能是若干个简正振动的线性叠加.例5.2（根据定义计算特征值、特征向量）：求矩阵A =[5−1−131−14−21]的特征值和特征向量.[解]: 矩阵A 的特征方程为：det (λI −A )=|λ−511−3λ−11−42λ−1|=(λ−3)(λ−2)2=0故A 的特征值为λ1=3，λ2=2（二重特征值）.当λ=λ1=3时，由(λI −A)x =0，得到方程[−211−321−422][x 1x 2x 3]=[000]它有无穷多个解，若假设x 1=1, 则求出解为x =[1,1,1]T ，记为x 1，则x 1是λ1对应的一个特征向量.当λ=λ2=2时，由(λI −A)x =0，得到方程[−311−311−421][x 1x 2x 3]=[000]它有无穷多个解，若假设x 1=1, 则求出解为x =[1,1,2]T ，记为x 2，则x 2是λ2对应的一个特② 本书第八章将介绍这种常微分方程组的数值求解方法.图5-1 弹簧－质点系统.征向量.下面概括地介绍有关矩阵特征值、特征向量的一些性质，它们可根据定义5.1，以及公式(5.2)加以证明.定理5.1：设λj (j =1,2,…,n)为n 阶矩阵A 的特征值，则(1) ∑λj n j=1=∑a jj n j=1=tr(A) ;(2) ∏λj n j=1=det(A) .这里tr(A)表示矩阵对角线上元素之和，称为矩阵的迹（trace ）.从上述结论(2)也可以看出，非奇异矩阵特征值均不为0, 而0一定是奇异矩阵的特征值. 定理5.2：矩阵转置不改变特征值，即λ(A )=λ(A T ).定理5.3：若矩阵A 为对角阵或上(下)三角阵，则其对角线元素即为矩阵的特征值.定理5.4：若矩阵A 为分块对角阵，或分块上(下)三角阵，例如A =[A 11A 12⋯A 1m A 22⋯A 2m ⋱⋮A mm] , 其中每个对角块A jj 均为方阵，则矩阵A 的特征值为各对角块矩阵特征值的合并，即λ(A )=⋃λ(A jj )m j=1.定理5.5：矩阵的相似变换(similarity transformation)不改变特征值. 设矩阵A 和B 为相似矩阵，即存在非奇异矩阵X 使得B =X −1AX ，则(1) 矩阵A 和B 的特征值相等，即 λ(A )=λ(B ) ;(2) 若y 为B 的特征向量，则相应地，Xy 为A 的特征向量.通过相似变换并不总能把矩阵转化为对角阵，或者说矩阵A 并不总是可对角化的(diagonalizable). 下面给出特征值的代数重数、几何重数，和亏损矩阵的概念，以及几个定理..定义5.2: 设矩阵A ∈ℝn×n 有m 个（m n ）不同的特征值λ̃1,⋯,λ̃m ，若λ̃j 是特征方程的n j 重根，则称n j 为λ̃j 的代数重数(algebraic multiplicity)，并称λ̃j 的特征子空间(ℂn 的子空间)的维数为λ̃j 的几何重数(geometric multiplicity). 定理5.6：设矩阵A ∈ℝn×n 的m 个不同的特征值为λ̃1,⋯,λ̃m ，特征值λ̃j ,(j =1,⋯,m)的代数重数为n j ，几何重数为k j ，则(1) ∑n j m j=1=n ，且任一个特征值的几何重数不大于代数重数，即∀j ，n j ≥k j .(2) 不同特征值的特征向量线性无关，并且将所有特征子空间的∑k j m j=1个基(特征向量)放在一起，它们构成一组线性无关向量.(3) 若每个特征值的代数重数等于几何重数，则总共可得n 个线性无关的特征向量，它们是全空间ℂn 的基.定义5.3：若矩阵A ∈ℝn×n 的某个代数重数为k 的特征值对应的线性无关特征向量数目少于k （即几何重数小于代数重数），则称A 为亏损阵（defective matrix ），否则称其为非亏损阵（nondefective matrix ）.定理5.7：设矩阵A ∈ℝn×n 可对角化，即存在非奇异矩阵X ∈ℂn×n 使得X −1AX =Λ,其中Λ∈ℂn×n 为对角阵, 的充要条件是A 为非亏损矩阵. 此时，Λ的对角线元素为矩阵A 的特征值，而矩阵X 的列向量为n 个线性无关的特征向量.定理5.7中方程的等价形式为A =XΛX −1, 它被称为特征值分解，也叫谱分解(spectrum decomposition). 特征值分解存在的充要条件是A 为非亏损矩阵. 但现实中还有很多矩阵是亏损矩阵，例如例5.2中的矩阵，它的特征值2的代数重数为2，而几何重数仅为1. 这种矩阵不能相似变换为对角阵，但存在下面的若当分解(Jordan decomposition).定理5.8：设矩阵A ∈ℝn×n , 存在非奇异矩阵X ∈ℂn×n 使得A =XJX −1,矩阵J 为形如[J 1⋱J p ]的分块对角阵（称为若当标准型），其中J k =[ λk 1λk ⋱⋱1λk ] 称为若当块，其对角线元素为矩阵A 的特征值. 设矩阵A 有m 个不同的特征值为λ̃1,⋯,λ̃m ，特征值λ̃j ,(j =1,⋯,m)的代数重数为n j ，几何重数为k j ，则p =∑k j m j=1, λ̃j 对应于k j 个若当块, 其阶数之和等于n j .在若当分解中，如果所有若当块都是1阶的，则J 为对角阵，这种分解就是特征值分解，相应的矩阵为非亏损阵. 若当分解是很有用的理论工具，利用它还可证明下面关于矩阵运算结果的特征值的定理.定理5.9：设λj (j =1,2,…,n)为n 阶矩阵A 的特征值，则(1) 矩阵cA, c 为常数, 的特征值为cλ1,cλ2,⋯,cλn .(2) 矩阵A +pI, p 为常数, 的特征值为λ1+p,λ2+p,⋯,λn +p.(3) 矩阵A k , k 为正整数, 的特征值为λ1k ,λ2k ,⋯,λn k .(4) 设p (t )为一多项式函数，则矩阵p (A )的特征值为p (λ1),p (λ2),⋯ ,p (λn ) .(5) 若A 为非奇异矩阵，则λj ≠0,(j =1,2,…,n), 且矩阵A −1的特征值为λ1−1,λ2−1,⋯,λn −1.5.1.2特征值分布范围的估计估计特征值的分布范围或它们的界，无论在理论上或实际应用上，都有重要意义. 比如，本书前面的内容曾涉及两个问题：(1). 计算矩阵的2-条件数：cond (A )2=√λmax (A T A)λmin (A T A) ;(2). 考察一阶定常迭代法x (k+1)=Bx (k)+f 的收敛性、收敛速度：收敛的判据是谱半径ρ(B)=max 1≤j≤n |λj (B)|<1 ; 收敛速度为R =−log 10ρ(B) .其中都需要对矩阵特征值分布范围的了解.上一章的定理4.4说明谱半径的大小不超过任何一种算子范数，即ρ(A )≤‖A ‖ ,这是关于特征值的上界的一个重要结论.下面先给出定义5.4，再介绍有关特征值的界的另一个重要结论.定义5.4：设A =(a kj )∈ℂn×n ，记r k =∑|a kj |n j=1j≠k ,(k =1,⋯,n)，则集合D k ={z||z −a kk |≤r k ,z ∈ℂ},(k =1,⋯,n)在复平面为以a kk 为圆心、r k 为半径的圆盘，称为A 的Gerschgorin （格什戈林）圆盘.图5-2显示了一个3⨯3复矩阵的格什戈林圆盘.定理5.10 (圆盘定理)：设A =(a kj )∈ℂn×n ，则：(1) A 的每一个特征值必属于A 的格什戈林圆盘之中，即对任一特征值λ必定存在k,1≤k ≤n ，使得：|λ−a kk |≤∑|a kj |nj=1j≠k .(5.4)图5-2 复坐标平面，以及3⨯3矩阵A 的格什戈林圆盘.用集合的关系来说明，这意味着λ(A)⊆⋃D k n k=1.(2) 若A 的格什戈林圆盘中有m 个组成一连通并集S ，且S 与余下的n −m 个圆盘分离，则S内恰好包含A 的m 个特征值（重特征值按重数计）.对图5-2所示的例子，定理5.10的第(2)个结论的含义是：D 1中只包含一个特征值，而另外两个特征值在D 2,D 3的并集中. 下面对定理5.10的结论(1)进行证明，结论(2)的证明超出了本书的范围.[证明]: 设λ为A 的任一特征值，则有Ax =λx ，x 为非零向量. 设x 中第k 个分量最大，即|x k |=max 1≤j≤n|x j |>0 , 考虑方程(5.2)中第k 个方程：∑a kj x j nj=1=λx k , 将其中与x k 有关的项移到等号左边，其余到右边，再两边取模得：|λ−a kk ||x k |=|∑a kj x j n j=1j≠k |≤∑|a kj ||x j |n j=1j≠k ≤|x k |∑|a kj |nj=1j≠k .(5.5)最后一个不等式的推导利用了“x 中第k 个分量最大”的假设. 将不等式(5.5)除以|x k |，即得到(5.4)式，因此证明了定理 5.10的结论(1). 上述证明过程还说明，若某个特征向量的第k 个分量的模最大，则相应的特征值必定属于第k 个圆盘中.根据定理5.2，还可以按照矩阵的每一列元素定义n 个圆盘，对于它们定理5.10仍然成立. 下面的定理是圆盘定理的重要推论，其证明留给感兴趣的读者.定理5.11：设A ∈ℝn×n ，且A 的对角元均大于0，则(1) 若A 严格对角占优，则A 的特征值的实部都大于0.(2) 若A 为对角占优的对称矩阵，则A 一定是对称半正定矩阵，若同时A 非奇异，则A 为对称正定矩阵.例5.3 (圆盘定理的应用)：试估计矩阵A =[41010−111−4]的特征值范围.[解]: 直接应用圆盘定理，该矩阵的三个圆盘如下:D 1: |λ−4|≤1, D 2: |λ|≤2, D 3: |λ+4|≤2.D 1与其他圆盘分离，则它仅含一个特征值，且必定为实数（若为虚数则其共轭也是特征值，这与D 1仅含一个特征值矛盾）. 所以对矩阵特征值的范围的估计是：3≤λ1≤5,λ2,λ3∈D 2∪D 3 .再对矩阵A T 应用圆盘定理，则可以进一步优化上述结果. 矩阵A T 对应的三个圆盘为： D ’1: |λ−4|≤2, D ’2: |λ|≤2, D ’3: |λ+4|≤1.这说明D ’3中存在一个特征值，且为实数，它属于区间[-5, -3]，经过综合分析可知三个特征值均为实数，它们的范围是：λ1∈[3,5],λ2∈[−2,2],λ3∈[−5,−3].事实上，使用Matlab 的eig 命令可求出矩阵A 的特征值为：4.2030, -0.4429, -3.7601.根据定理5.5，还可以对矩阵A 做简单的相似变换，例如取X 为对角阵，然后再应用圆盘定理估计特征值的范围.例5.4 (特征值范围的估计)：选取适当的矩阵X ，应用定理5.5和5.10估计例5.3中矩阵的特征值范围.[解]: 取X−1=[100010000.9] , 则A 1=X −1AX =[41010−109⁄0.90.9−4]的特征值与A 的相同. 对A 1应用圆盘定理，得到三个分离的圆盘，它们分别包含一个实特征值，由此得到特征值的范围估计：λ1∈[3,5],λ2∈[−199,199],λ3∈[−5.8,−2.2]. 此外，还可进一步估计ρ(A)的范围，即3≤ρ(A)≤5.8 .上述例子表明，综合运用圆盘定理和矩阵特征值的性质（如定理5.2, 定理5.5），可对特征值的范围进行一定的估计. 对具体例子，可适当设置相似变换矩阵，尽可能让圆盘相互分离，从而提高估计的有效性.5.2幂法与反幂法幂法是一种计算矩阵最大的特征值及其对应特征向量的方法. 本节介绍幂法、反幂法以及加快幂法迭代收敛的技术.5.2.1幂法定义5.5：在矩阵A 的特征值中，模最大的特征值称为主特征值，也叫“第一特征值”，它对应的特征向量称为主特征向量.应注意的是，主特征值有可能不唯一，因为模相同的复数可以有很多. 例如模为5的特征值可能是5,−5,3+4i,3−4i , 等等. 另外，请注意谱半径和主特征值的区别.如果矩阵A 有唯一的主特征值，则一般通过幂法能方便地计算出主特征值及其对应的特征向量. 对于实矩阵，这个唯一的主特征值显然是实数，但不排除它是重特征值的情况. 幂法(power iteration)的计算过程是，首先任取一非零向量v 0∈ℝn ，再进行迭代计算：v k =Av k−1,(k =1,2,⋯)得到向量序列{v k }，根据它即可求出主特征与特征向量. 下面用定理来说明.定理5.12: 设A ∈ℝn×n ，其主特征值唯一，记为λ1，且λ1的几何重数等于代数重数，则对于非零向量v 0∈ℝn ，v 0不与主特征值对应的特征向量正交，按迭代公式进行计算：v k =Av k−1,(k =1,2,⋯),存在如下极限等式：lim k→∞v k λ1k =x 1 , (5.6) lim k→∞(v k+1)j (v k )j =λ1 , (5.7)其中x 1为主特征向量，(v k )j 表示向量v k 的第j 个分量(k =1,2,⋯).[证明]: 为了推导简便，不妨设主特征值λ1不是重特征值，并且假设矩阵A 为非亏损矩阵. 设A 的n 个特征值按模从大到小排列为: |λ1|>|λ2|≥⋯≥|λn |，它们对应于一组线性无关的单位特征向量x ̂1,⋯,x ̂n . 向量v 0可写成这些特征向量的线性组合：v 0=α1x̂1+⋯+αn x ̂n 根据已知条件，α1≠0，则v k =Av k−1=A k v 0=α1λ1k x ̂1+α2λ2k x̂2+⋯+αn λn k x ̂n =λ1k [α1x ̂1+∑αj (λj λ1)kx ̂j n j=2] =λ1k (α1x̂1+εk ) 其中εk =∑αj (λj λ1)k x ̂j n j=2. 由于|λj λ1|<1,(j =2,…,n), 则 lim k→∞εk =0 ⟹lim k→∞v kλ1k =α1x̂1 . 由于特征向量放大、缩小任意倍数后仍是特征向量，设x 1=α1x̂1，则它是主特征对应的一个特征向量. 上式说明，随k 的增大， v k 越来越趋近于主特征值的对应的特征向量.设j 为1到n 之间的整数，且(v k )j ≠0，则(v k+1)j (v k )j =λ1(α1x ̂1+εk+1)j (α1x̂1+εk )j 由于lim k→∞εk =0，随k 的增大上式等号右边趋于一个常数: λ1. 这就证明了定理的结论.若矩阵A 为亏损矩阵，可利用矩阵的若当分解证明这个定理，这里略去. 在这种情况下，“主特征值的几何重数等于代数重数”这一条件很重要，例如，若A =[310030001] ,它的主特征值为3，但其几何重数为1，不满足条件. 对这个矩阵A 进行实验显示无法用幂法求出主特征值.关于定理5.12，再说明几点：● 当主特征值λ1为重特征值时，应要求其几何重数等于代数重数，此时特征子空间维数大于1，向量序列{v k λ1k ⁄}的收敛值是其特征子空间中的某一个基向量.● 公式(5.7)式的含义是相邻迭代向量分量的比值收敛到主特征值. 因此在实际计算时，可任意取j 的值，只需保证比值的分母不为零.● 证明中假设了α1≠0，在实际应用中往往随机选取v 0，由于存在舍入误差，它一般都能满足. 感兴趣的读者也可思考一下，若初始向量v 0恰好与主特征向量都正交，那么幂法中的迭代向量序列会有什么结果？直接使用幂法，还存在如下两方面问题：(1) 溢出：由于v k ≈λ1k x 1，则|λ1|>1时，实际计算v k 会出现上溢出（当k 很大时）；|λ1|<1时，实际计算v k 会出现下溢出（当k 很大时）.(2) 可能收敛速度很慢. 由于εk =∑αj (λj λ1)kx j n j=2， εk →0的速度取决于求和式中衰减最慢的因子|λ2λ1|，当|λ2λ1|≈1时，收敛很慢. 由此导致v k →λ1k α1x 1, (v k+1)j (v k )j →λ1的收敛速度都将很慢，严重影响计算的效率.下面采用规格化向量的技术防止溢出，导出实用的幂法. 关于加速收敛技术的讨论，见下一小节.定义 5.6：记max ̅̅̅̅̅̅(v )为向量v ∈ℝn 的绝对值最大的分量， max ̅̅̅̅̅̅(v )=v j ，其中j 满足|v j |=max 1≤k≤n |v k |, 若j 的值不唯一，则取最小的那个. 并且，称u =v/max ̅̅̅̅̅̅(v )为向量v 的规格化向量(normalized vector).例5.5（规格化向量）：设v =[3,−5,0]T ，max ̅̅̅̅̅̅(v )=−5，对应的规格化向量为u =[−35,1,0]T .根据定义5.6，容易得出规格化向量的两条性质.定理5.13: 定义5.6中的规格化向量满足如下两条性质：(1) 若u 为规格化向量，则‖u ‖ =1，并且max ̅̅̅̅̅̅(u )=1.(2) 设向量v 1和v 2的规格化向量分别为u 1和u 2，若v 1=αv 2, 实数α≠0，则u 1= u 2.在幂法的每一步增加向量规格化的操作可解决溢出问题. 先看第一步，v 1=Av 0，此时计算v 1的规格化向量u 1=v 1max ̅̅̅̅̅̅(v 1)=Av 0max ̅̅̅̅̅̅(Av 0). 然后使用规格化向量计算v 2:v 2=Au 1=A 2v 0max ̅̅̅̅̅̅(Av 0), (5.8) 再进行向量规划化操作，u 2=v 2max ̅̅̅̅̅̅(v 2)=A 2v 0max ̅̅̅̅̅̅(A 2v 0). (5.9) 公式(5.9)的推导，利用了(5.8)式和定理5.13的结论(2). 依次类推，我们得到： { v k =Au k−1=A k v 0max ̅̅̅̅̅̅(A k−1v 0) u k =v k max ̅̅̅̅̅̅(v k )=A k v 0max ̅̅̅̅̅̅(A k v 0) , k =1,2,⋯. (5.10) 根据定理5.12的证明过程， A k v 0=λ1k [α1x ̂1+∑αj (λj λ1)k x ̂j n j=2] ⟹u k =A k v 0max ̅̅̅̅̅̅(A k v 0)=α1x ̂1+∑αj (λj λ1)k x ̂j n j=2max ̅̅̅̅̅̅(α1x ̂1+∑αj (λj λ1)k x ̂j n j=2)k→∞→ x 1max ̅̅̅̅̅̅(x 1) , 即u k 逐渐逼近规格化的主特征向量. 同理，v k =Au k−1=A k v 0max ̅̅̅̅̅̅(A k−1v 0)=λ1k [α1x ̂1+∑αj (λj λ1)k x ̂j n j=2]max ̅̅̅̅̅̅(λ1k−1[α1x ̂1+∑αj (λj λ1)k−1x̂j n j=2]) =λ1α1x ̂1+∑αj(λj λ1)kx ̂j n j=2max ̅̅̅̅̅̅(α1x ̂1+∑αj (λj λ1)k−1x ̂j n j=2) 因此，根据定理5.13的结论(1)有：lim k→∞v k=λ1x1max̅̅̅̅̅̅(x1)⟹limk→∞max̅̅̅̅̅̅(v k)=λ1.基于上述推导，我们得到如下定理，以及如算法5.1描述的实用幂法.定理5.14: 设A∈ℝn×n，其主特征值唯一（且几何重数等于代数重数），记为λ1，取任意非零初始向量v0=u0，它不与主特征值对应的特征向量正交，按迭代公式(5.10)进行计算，则lim k→∞u k=x1max̅̅̅̅̅̅(x1),(5.11)lim k→∞max̅̅̅̅̅̅(v k)=λ1 ,(5.12)其中x1为主特征向量.算法5.1：计算主特征值λ1和主特征向量x1的实用幂法输入：v,A; 输出：x1,λ1.u:=v;While不满足判停准则dov:=Au;λ1:=max̅̅̅̅̅̅(v); {主特征值近似值}u:=v/λ1; {规格化}Endx1:=u. {规格化的主特征向量}在算法5.1中，可根据相邻两步迭代得到的主特征值近似值之差来判断是否停止迭代. 每个迭代步的主要计算是算一次矩阵与向量乘法，若A为稀疏矩阵则可利用它的稀疏性提高计算效率. 实用的幂法保证了向量序列{v k},{u k}不溢出，并且向量v k的最大分量的极限就是主特征值.最后，针对幂法的适用范围再说明两点：(1). 若实矩阵A对称半正定或对称半负定，则其主特征值必唯一（而且是非亏损阵）. 有时也可以估计特征值的分布范围，从而说明主特征值的唯一性. 只有满足此条件，才能保证幂法的收敛性.(2). 对一般的矩阵，幂法的迭代过程有可能不收敛，此时序列{u k}有可能包括多个收敛于不同向量的子序列，它趋向于成为多个特征向量的线性组合. 但是，一旦幂法的迭代过程收敛，向量序列的收敛值就一定是特征向量，并可求出相应的特征值.例5.6 (实用的幂法)：用实用的幂法求如下矩阵的主特征值：A=[3113] ,[解]: 取初始向量为v0=u0=[01]T . 按算法5.1的迭代过程，计算结果列于表5-1中.表5-1 实用幂法的迭代计算过程从结果可以看出，在每次迭代步中做的规格化操作避免了分量的指数增大或缩小. 经过9步迭代，特征值max ̅̅̅̅̅̅(v k )已非常接近主特征值的准确值4，特征向量也非常接近[1 1]T .5.2.2加速收敛的方法 加速幂法迭代收敛过程的方法主要有两种：原点位移技术和瑞利商（Rayleigh quotient ）加速. 下面做些简略的介绍.一. 原点位移技术原点位移技术，也叫原点平移技术，它利用定理5.9的结论(2)，即矩阵A −pI 的特征值为A 的特征值减去p 的结果. 对矩阵B =A −pI 应用幂法有可能得到矩阵A 的某个特征值λj 和相应的特征向量. 要使原点位移达到理想的效果，首先要求λj −p 是B 的主特征值，其次还要使幂法尽快收敛，即比例|λ2(B)λj −p |要尽量小，这里的λ2(B)表示矩阵B 的（按模）第二大的特征值. 在某种情况下设置合适的p 值，矩阵A,B 可同时取到主特征值. 图5-3显示了这样一个例子，矩阵A 的特征值分布在阴影区域覆盖的实数轴上，λ1为其主特征值. 按图中所示选取的p 值，将使得λ1−p 是矩阵B =A −pI 的主特征值，并且显然有|λ2(B)λ1−p |<|λ2(A)λ1| . 此时用幂法计算B 的主特征值能更快地收敛，进而得到矩阵的A 的主特征值. 图5-3也解释了原点位移法名字的由来，即将原点（或虚数坐标轴）移到p 的位置上，原始矩阵A 的特征值分布变成了矩阵B 的特征值分布.采用原点位移技术后，执行幂法仅带来很少的额外运算，而且仍然能利用矩阵A 的稀疏性. 它的关键问题是，如何选择合适的参数p 以达到较好的效果？这依赖于具体矩阵的情况，以及对其特征值分布的了解. 在后面，我们还会看到原点位移技术的其他用途.二. 瑞利商加速首先给出瑞利商的定义，以及它与特征值的关系，然后介绍瑞利商加速技术.定义5.7：设A ∈ℝn×n ，且为对称矩阵，对任一非零向量x ≠0，称R (x )=〈Ax,x 〉〈x,x 〉为对应于向量x 的瑞利商（Rayleigh quotient ）. 这里符号〈,〉代表向量内积.定理5.15：设A ∈ℝn×n ，且为对称矩阵，其n 个特征值依次为：λ1≥λ2≥⋯≥ λn ，则矩阵A 有关的瑞利商的上下确界分别为λ1和λn . 即∀x ≠0,λn ≤R (x )≤λ1,且当x 为λ1对应的特征向量时R (x )=λ1，当x 为λn 对应的特征向量时R (x )=λn .[证明]: 根据实对称矩阵的特点，即可正交对角化（定理3.3），设特征值λ1,λ2,⋯,λn 对应的单位特征向量为x 1,x 2,⋯,x n ，设x =∑αj x j n j=1，则〈x,x 〉=〈∑αj x j n j=1,∑αj x j n j=1〉=∑αj 2n j=1，而图5-3 原点位移技术示意图.。
并行计算-矩阵特征值计算--
9 矩阵特征值计算在实际的工程计算中，经常会遇到求n 阶方阵 A 的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。
对于一个方阵A，如果数值λ使方程组Ax=λx即(A-λI n )x=0 有非零解向量(Solution Vector)x，则称λ为方阵A的特征值，而非零向量x为特征值λ所对应的特征向量，其中I n 为n阶单位矩阵。
由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。
本章主要介绍求一般方阵绝对值最大的特征值的乘幂(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-side Rotation)法以及求一般矩阵全部特征值的QR 方法及一些相关的并行算法。
1.1 求解矩阵最大特征值的乘幂法1.1.1 乘幂法及其串行算法在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。
乘幂法是一种求矩阵绝对值最大的特征值的方法。
记实方阵A的n个特征值为λi i=(1,2, …,n)，且满足：│λ1 │≥│λ2 │≥│λ3 │≥…≥│λn │特征值λi 对应的特征向量为x i 。
乘幂法的做法是：①取n维非零向量v0 作为初始向量；②对于k=1,2, …,做如下迭代：直至 uk +1 ∞- u k u k =Av k -1 v k = u k /║u k ║∞ < ε为止，这时v k +1 就是A 的绝对值最大的特征值λ1 所对应的特征向∞量x 1 。
若v k -1 与v k 的各个分量同号且成比例，则λ1 =║u k ║∞ ；若v k -1 与v k 的各个分量异号且成比 例，则λ1 = -║u k ║∞ 。
若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时 间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格 化操作，所以下述乘幂法串行算法 21.1 的一轮计算时间为n 2+2n =O(n 2 )。
一种厄米特矩阵特征值并行求解方法
一种厄米特矩阵特征值并行求解方法作者：严新文彭永健张续莹来源：《科技资讯》2018年第11期摘要：厄米特矩阵特征值特征向量快速求解是在实际工程中常常会碰到的问题，但大部分基于计算机实现的算法都是利用串行思想进行编程求解，不能满足实时求解需求。
本文提出基于FPGA的并行算法来求解本类问题。
首先，将厄米特矩阵实数化。
其次，通过并行雅克比分解和特征值快速排序来恢复厄米特矩阵的特征值特征向量。
最后，通过MATLAB仿真验证了算法的有效性，并通过蒙特卡洛仿真确定了雅克比sweep的迭代次数，同时通过FPGA实现验证了算法的实时性。
关键词：厄米特矩阵特征值并行雅克比分解蒙特卡洛仿真中图分类号：O129 文献标识码：A 文章编号：1672-3791（2018）04（b）-0108-02矩阵特征值特征向量求解一直是线性代数中最重要的操作之一，在科学计算中有着非常广泛的应用[1]。
工程技术和物理等领域中的很多问题在数学上都归结为求矩阵的特征值特征向量问题。
例如振动问题和物理学中某些临界值的确定，均归结为求矩阵的特征值问题。
矩阵特征值和特征向量的求解方法有很多，幂法可以求解一般矩阵绝对值最大的特征值，雅可比方法可以求解对称矩阵的全部特征值，QR方法可以求解一般实矩阵的全部特征值[2]。
在一般的工程当中，常用的矩阵特征分解方法大部分都是针对实矩阵的，而在实际工程中，我们常需要处理复矩阵，而且大部分情况下为一个厄米特矩阵，由于厄米特矩阵具有对称性，能简化其特征值特征向量求解过程。
但是上述算法的计算机实现大部分都是利用串行思想进行编程求解，并不适合用于并行计算，本文详细介绍一种可用于FPGA并行高效实现的厄米特矩阵特征值特征向量求解方法[3]。
1 厄米特矩阵特征值求解方法1.1 实数化如果直接对厄米特矩阵进行特征分解显然是不合适的，也是没有效率的，因此我们在对其进行特征分解前先进行一项预处理过程，即将它实数化。
1.2 并行雅克比分解对于任意的一个实对称矩阵A，只要能够求得一个正交矩阵U，使得成为一个对角阵D，就得到了A的所有特征值和对应的特征向量。
实对称矩阵特征值分解高速并行算法的fpga实现
实对称矩阵特征值分解高速并行算法的fpga实现高速并行算法在计算机视觉中有着重要作用，不仅可提高运算速度，而且可以节省计算资源。
随着硬件技术的发展，fpga的优势在数值运算领域也越来越明显。
本文研究fpga实现对称矩阵特征值分解高速并行算法。
本文根据传统的特征值分解算法构造特征值分解算法流水线，利用fpga模块化结构实现模块的单周期执行，实现对称矩阵特征值分解高速并行算法。
首先，给出算法原理和流程：（1）特征值分解算法以矩阵A的n阶为输入，求解矩阵A的n个主特征值和对应的特征向量。
（2）算法在单指令多数据模式下并行执行，把工作分解入多个分组，在每个分组内并行计算，每次迭代依次更新下阶段分组所需要的数据，从而充分利用多核计算资源。
（3）算法利用 fpga 的优势，针对每一阶分组构造一个模块，实现模块的单周期执行，以实现高速并行的特征值分解算法。
其次，算法的逻辑设计和实现，主要有以下几步：（1）设计特征值分解模块的框架，完成模块的调用连接方式。
（2）根据特征值分解算法构造特征值分解算法流水线，根据不同的阶数应用相应的指令，实现特征值分解模块，形成模块化结构。
（3）将计算模块输出的结果和输入线、控制线以及控制信号联系起来，以实现算法流水线控制。
（4）实现计算结果存储方式，支持特征值和特征向量的对称矩阵特征值置换保存。
（5）根据需要的高性能运算，设计适当的特性模块，针对特征值分解算法构造fpga 实现模块。
最后，将特征值分解算法集成于定制的fpga硬件设备中，最终实现fpga实现对称矩阵特征值分解算法，提高算法的计算速度，节省能源。
以上就是fpga实现对称矩阵特征值分解高速并行算法，结合fpga 最大限度发挥多核计算资源，通过模块化结构实现单周期执行，从而达到高速并行算法、节省能源的目的。
矩阵特征值分解算法实现的并行化与优化
矩阵特征值分解算法实现的并行化与优化矩阵特征值分解是一种重要的数值计算方法，可以将一个矩阵拆解为特征向量和对应的特征值。
这一算法的应用广泛，涵盖了数学、物理、工程等多个领域。
然而，矩阵特征值分解的计算复杂度较高，限制了其在大规模数据和实时计算场景中的应用。
因此，实现并行化与优化的研究成为目前的热点之一。
在本文中，我们将探讨矩阵特征值分解算法的并行化与优化方法。
首先，我们介绍传统的特征值分解算法，然后讨论并行化的思路和方法。
最后，我们将针对不同的优化需求，提出相应的优化策略。
一、传统的特征值分解算法传统的特征值分解算法主要有幂法、QR方法和雅可比方法等。
这些方法在串行计算环境下运行，计算复杂度较高。
幂法是最简单的特征值分解算法，但它在实际应用中存在精度不高和收敛速度慢的问题。
QR方法通过将矩阵分解为正交矩阵和上三角矩阵的乘积，然后迭代求得特征值。
雅可比方法将矩阵通过相似变换转化为对角矩阵，再求得其特征值。
然而，这些方法在大规模矩阵计算时，计算时间较长，效率较低。
二、并行化思路和方法为了提高特征值分解算法的计算效率，研究者们开始将并行计算引入其中。
并行化可以通过使用多个处理器或者分布式计算系统来实现。
其中，多核并行化和GPU并行化是两种主要的方向。
1. 多核并行化多核并行化的思路是利用多个处理器并行计算。
可以将矩阵分成多个子矩阵，然后分别在不同的处理器上计算特征向量和特征值。
这样可以大大缩短计算时间。
此外，还可以采用线程池和任务调度等技术来合理管理计算资源，提高并行计算效率。
2. GPU并行化GPU并行化是利用图形处理器的并行计算能力来加速特征值分解算法。
GPU具有大量的计算单元和高带宽的内存，适合于并行计算。
可以将矩阵的计算任务分配给不同的GPU核心，并采用适当的并行算法来提高计算效率。
此外，还可以通过利用GPU的共享内存和纹理内存等特性来加速计算过程。
三、优化策略除了并行化计算，还有一些优化策略可以提高矩阵特征值分解算法的计算效率。
矩阵特征值分解算法的并行实现优化策略
矩阵特征值分解算法的并行实现优化策略矩阵特征值分解是线性代数中一个重要的问题，它在众多科学与工程领域中都有广泛的应用。
然而，传统的特征值分解算法在处理大规模矩阵时，往往需要消耗大量的计算资源和时间。
为了提高特征值分解算法的效率和性能，研究人员致力于开发并行实现优化策略。
本文将讨论一些常用的并行实现优化策略，并探讨它们的优势和局限性。
一、基于任务并行的策略任务并行是一种将计算任务分解为多个子任务，每个子任务在独立的处理器上并行执行的策略。
在矩阵特征值分解算法中，任务并行可以应用于矩阵的乘法和特征向量的计算。
例如，可以将矩阵分割成多个子矩阵，并将每个子矩阵分配给不同的处理器进行计算。
这样可以极大地减少计算时间，提高算法的效率。
然而，任务并行也存在一些问题，比如任务分配、负载均衡等，需要通过合理的调度算法进行解决。
二、基于数据并行的策略数据并行是一种将数据划分为多个部分，每个部分在不同的处理器上并行处理的策略。
在矩阵特征值分解算法中，可以将矩阵的行或列划分为多个部分，并分配给不同的处理器进行计算。
这种策略可以充分利用计算资源，提高算法的并行度，从而加速特征值分解的过程。
然而，数据并行也会引入通信开销和数据重组的问题，需要采取相应的优化措施来减少这些开销。
三、混合并行策略混合并行是一种将任务并行和数据并行相结合的策略，旨在充分发挥两者的优势并弥补各自的劣势。
在矩阵特征值分解算法中，可以将矩阵的乘法操作使用任务并行实现，而将特征向量的计算使用数据并行实现。
这样可以在保证算法的高效性的同时，最大限度地利用计算资源，提高算法的性能。
四、硬件优化策略除了并行实现策略外，还可以通过硬件优化来进一步提高矩阵特征值分解算法的性能。
例如，可以使用并行计算硬件，如GPU（图形处理器）进行计算加速。
GPU具有大量的并行计算单元，适合并行计算密集型任务，可以显著提高特征值分解算法的计算速度。
此外，还可以使用优化的存储器结构、算法优化等手段来提高算法的效率。
Jacobi矩阵特征值的并行算法
Jacobi矩阵特征值的并行算法
刘艳红;吕全义
【期刊名称】《纺织高校基础科学学报》
【年(卷),期】2011(024)001
【摘要】提出了并行求解实三对角矩阵特征值方法,该方法主要针对Jacobi矩阵.应用求多项式根的Sturm法,将矩阵特征多项式的求根区间隔离成单根区间;对已隔离出的单根区间先用二分法求解,达到一定精度后再用牛顿法精确求解.考虑到处理机负载平衡问题,将求根区间分成若干等分,然后按区间循环地将其分给各个处理机.各处理机并行地进行求根计算,它们之间无通信.通过此方法实现了处理机负载平衡,算法并行效率达0.85以上.数值算例表明了此并行算法的高效性.
【总页数】5页(P21-25)
【作者】刘艳红;吕全义
【作者单位】西北工业大学,应用数学系,陕西,西安,710129;西北工业大学,应用数学系,陕西,西安,710129
【正文语种】中文
【中图分类】O246
【相关文献】
1.由部分特征值和顺序主子阵构造广义Jacobi矩阵的逆特征值问题 [J], 潘云兰;秦立
2.由特征值和顺序主子阵构造广义Jacobi矩阵的逆特征值问题 [J], 徐秀斌;秦立
3.一类广义Jacobi矩阵的逆特征值问题 [J], 孟纯军;姜婷婷
4.由混合数构造的伪Jacobi矩阵的广义逆特征值问题 [J], 薛昕; 雷英杰; 郑志勇
5.由混合数构造的伪Jacobi矩阵的广义逆特征值问题 [J], 薛昕; 雷英杰; 郑志勇因版权原因，仅展示原文概要，查看原文内容请购买。
矩阵特征值分解算法的优化与并行实现
矩阵特征值分解算法的优化与并行实现矩阵特征值分解是一种重要的数学问题，在科学计算、数据处理、图像处理等领域都有广泛的应用。
特征值分解将一个矩阵分解为特征向量矩阵和特征值矩阵的乘积，它能够揭示矩阵的结构和特性，对于理解和分析矩阵的意义具有重要意义。
然而，由于特征值分解算法的高计算复杂度，效率成为制约其应用的一个瓶颈。
因此，研究和优化矩阵特征值分解算法，以提高其计算效率，对于推动相关领域的发展具有重要意义。
一、传统特征值分解算法的优化传统的特征值分解算法通常采用迭代方法，如幂迭代、反迭代、QR方法等。
这些算法的主要问题在于计算复杂度较高，尤其是在处理大规模矩阵时，计算时间会非常长。
因此，研究和优化这些传统算法成为降低计算复杂度的重要途径。
一种优化方法是加速收敛速度。
传统特征值分解算法通常需进行多次迭代才能得到精确的特征值和特征向量。
为了加速迭代过程，可以采用加速技术，如Krylov子空间方法和快速QR分解等。
这些方法能够通过寻找适当的迭代初始向量和加速技术，在较少的迭代次数内得到较好的近似解。
另一种优化方法是降低计算复杂度。
传统特征值分解算法通常需要进行大量的矩阵乘法和矩阵分解操作，这些操作的计算复杂度较高。
为了降低计算复杂度，可以采用矩阵分块技术、稠密矩阵的低秩近似和快速傅里叶变换等方法。
这些方法能够有效地减少计算量，并提高特征值分解算法的效率。
二、矩阵特征值分解算法的并行实现并行计算是提高特征值分解算法效率的重要手段。
由于特征值分解算法中的大量矩阵运算是相互独立的，因此可以通过并行计算来加速算法的执行。
并行实现可以利用多核处理器、GPU和分布式计算平台等并行计算设备来实现。
在并行计算中，可以采用各种并行策略来实现特征值分解算法。
例如，可以将矩阵分成多个小块，分配给不同的计算节点，并行地进行特征值分解计算。
另外，可以采用数据并行的方式，在不同的计算节点上并行地执行矩阵乘法和矩阵分解操作。
还可以采用任务并行的方式，在不同的计算节点上并行地执行不同的计算任务。
求解大型矩阵特征值问题的并行精化Davidson方法
第26卷第5期2009年10月工程数学学报CHINESE JOURNAL OF ENGINEERING MATHEMATICSVol.26No.5Oct.2009文章编号:1005-3085(2009)05-0922-07求解大型矩阵特征值问题的并行精化Davidson方法王顺绪1,2,戴华1(1-南京航空航天大学理学院，南京210016;2-淮海工学院数理系，连云港222005)摘要:针对共享主存的并行计算环境和微机网络并行计算环境，本文给出了求解大型稀疏对称矩阵的部分极端特征对的并行精化Davidson方法，分析了该法的内在并行性。
各处理器利用矩阵的行块和投影子空间的正交基所组成矩阵的行块进行运算，结合重新启动策略求解矩阵多个特征对的近似值，并用以计算某型号机翼的固有频率，在微机网络并行计算环境和拥有共享主存并行计算环境IBM-P650上进行了数值试验。
关键词:并行计算；特征值问题；Davidson方法；精化方法分类号:AMS(2000)65F15;65Y05中图分类号:O241;O246文献标识码:A1引言在科学和工程技术领域中，经常需要计算大型稀疏对称矩阵的若干个特征对。
解决这一问题的有效方法是正交投影方法，包括子空间迭代法，Lanczos方法，Arnoldi方法，Davidson方法和Jacobi-Davidson方法[1,2]。
为进一步提高Lanczos方法和Davidson方法的可靠性，Under-wood提出了块Lanczos方法[3]，戴等提出了预处理块Lanczos方法[4]，Crouzeix等人提出了块Davidson方法[5,6]。
用正交投影方法计算特征值问题Ax=λx(1)的部分特征对时，Ritz向量˜x i可能收敛很慢甚至发散，针对这一问题，贾提出了保留Ritz值˜λi，用满足(A−˜λi I)ˆx i 2=minx∈V i, x 2=1(A−˜λi I)x 2(2)的精化向量ˆx i代替Ritz向量˜x i的精化思想[7-12]，并继而指出(A−˜λi I)ˆx i 2=σmin(AV i−˜λi V i),(3)其中ˆx i=V iˆz i,ˆz i是矩阵(A−˜λi I)V i的最小奇异值σmin对应的右奇异向量，且 ˆz i 2=1。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
9 矩阵特征值计算在实际的工程计算中，经常会遇到求n 阶方阵 A 的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。
对于一个方阵A，如果数值λ使方程组Ax=λx即(A-λI n )x=0 有非零解向量(Solution Vector)x，则称λ为方阵A的特征值，而非零向量x为特征值λ所对应的特征向量，其中I n 为n阶单位矩阵。
由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。
本章主要介绍求一般方阵绝对值最大的特征值的乘幂(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-side Rotation)法以及求一般矩阵全部特征值的QR 方法及一些相关的并行算法。
1.1 求解矩阵最大特征值的乘幂法1.1.1 乘幂法及其串行算法在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。
乘幂法是一种求矩阵绝对值最大的特征值的方法。
记实方阵A的n个特征值为λi i=(1,2, …,n)，且满足：│λ1 │≥│λ2 │≥│λ3 │≥…≥│λn │特征值λi 对应的特征向量为x i 。
乘幂法的做法是：①取n维非零向量v0 作为初始向量；②对于k=1,2, …,做如下迭代：直至uk+1 ∞-uku k =Av k-1 v k = u k /║u k ║∞<ε为止，这时v k+1 就是A的绝对值最大的特征值λ1 所对应的特征向∞量x1 。
若v k-1 与v k 的各个分量同号且成比例，则λ1 =║u k ║∞；若v k-1 与v k 的各个分量异号且成比例，则λ1 = -║u k ║∞。
若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格化操作，所以下述乘幂法串行算法21.1 的一轮计算时间为n2+2n=O(n2 )。
算法21.1 单处理器上乘幂法求解矩阵最大特征值的算法输入：系数矩阵A n×n ，初始向量v n×1 ，ε输出：最大的特征值m axBeginwhile (│diff│>ε) do(1)for i=1 to n do(1.1)sum=0(1.2)for j= 1 to n dosum=sum+a[i,j]*x[j]end for(1.3)b[i]= sumend for(2)max=│b[1]│(3)for i=2 to n doif (│b[i]│>max) then max=│b[i]│ end ifend for(4)for i=1 to n dox[i] =b[i]/maxend for(5)diff=max-oldmax , oldmax=maxend whileEnd1.1.2 乘幂法并行算法乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。
设处理器个数为p，对矩阵A 按行划分为p 块，每块含有连续的m 行向量，这里m=⎡n/ p⎤，编号为i的处理器含有A的第i m 至第(i+1)m-1 行数据，(i=0,1, …,p-1)，初始向量v被广播给所有处理器。
各处理器并行地对存于局部存储器中a 的行块和向量v 做乘积操作，并求其m 维结果向量中的最大值localmax，然后通过归约操作的求最大值运算得到整个n 维结果向量中的最大值max 并广播给所有处理器，各处理器利用max 进行规格化操作。
最后通过扩展收集操作将各处理器中的m维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。
具体算法框架描述如下：算法21.2 乘幂法求解矩阵最大特征值的并行算法输入：系数矩阵A n×n ，初始向量v n×1 ，ε输出：最大的特征值m axBegin对所有处理器m y_rank(my_rank=0,…, p-1)同时执行如下的算法:while (│diff│>ε) do /* diff 为特征向量的各个分量新旧值之差的最大值*/ (1)for i=0 to m-1 do /*对所存的行计算特征向量的相应分量*/(1.1)sum=0(1.2)for j= 0 to n-1 dosum=sum+a[i,j]*x[j]end for(1.3)b[i]= sumend for(2)localmax=│b[0]│/*对所计算的特征向量的相应分量求新旧值之差的最大值l ocalmax */(3)for i=1 to m-1 doif (│b[i]│>localmax) then localmax=│b[i]│ end ifend for(4)用A llreduce 操作求出所有处理器中l ocalmax 值的最大值m ax并广播到所有处理器中(5)for i=0to m-1 do /*对所计算的特征向量归一化*/b [i ] =b [i ]/maxend for (6)用 A llgather 操作将各个处理器中计算出的特征向量的分量的新值组合并广播到 所有处理器中(7)diff=max -oldmax, oldmax =max end while End若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为 m n +2m ；一轮迭代中，各处理器做一次归约操作，通信量为 1， 一次扩展收集操作，通信量为 m ，则通信时间为 4t s (p - 1) + (m + 1)t w ( p - 1) 。
因此乘幂法的一轮并行计算时间为T p = mn + 2m + 4t s (MPI 源程序请参见所附光盘。
p - 1) + (m + 1)t w ( p - 1) 。
1.2 求对称矩阵特征值的雅可比法1.2.1 雅可比法求对称矩阵特征值的串行算法若矩阵A =[a ij ]是n 阶实对称矩阵，则存在一个正交矩阵U ，使得U T AU =D 其中 D 是一个对角矩阵，即⎡λ1⎢ D = ⎢0⎢ ⎢ ⎢⎣0λ20 ⎤⎥ 0 ⎥ ⎥ ⎥ λn ⎥⎦这里λi (i =1,2,…,n )是A 的特征值，U 的第i 列是与λi 对应的特征向量。
雅可比算法主要是通过 正交相似变换将一个实对称矩阵对角化，从而求出该矩阵的全部特征值和对应的特征向量。
因此可以用一系列的初等正交变换逐步消去A 的非对角线元素，从而使矩阵A 对角化。
设初 等正交矩阵为R (p ,q ,θ)，其(p ,p )( q ,q )位置的数据是cos θ，(p , q )( q , p )位置的数据分别是-sin θ和 sin θ(p ≠ q )，其它位置的数据和一个同阶数的单位矩阵相同。
显然可以得到：R (p ,q ,θ) T R (p ,q ,θ)=I n不妨记B= R (p ,q ,θ)T AR (p ,q ,θ)，如果将右端展开，则可知矩阵B 的元素与矩阵A 的元素之 间有如下关系：b pp = a pp cos 2θ+a qq sin 2θ+a pq sin2θ b qq = a pp sin 2θ+a qq cos 2θ-a pq sin2θb pq = ((a qq -a pp )sin2θ)/2+a pq cos2θ b qp = b pq b pj = a pj cos θ+a qj sin θ b qj = -a pj sin θ +a qj cos θ b ip = a ip cos θ+a iq sin θ b iq = -a ip sin θ +a iq cos θ b ij = a ij其中 1 ≤ i , j ≤ n 且i ,j ≠ p ,q 。
因为A 为对称矩阵，R 为正交矩阵，所以B 亦为对称矩阵。
若 要求矩阵B 的元素b pq =0，则只需令 ((a qq -a pp )sin2θ)/2+a pq cos2θ=0，即：pqqqtg 2θ =-a pq (a q q - a pp ) 2在实际应用时，考虑到并不需要解出 θ，而只需要求出 s in2θ，sin θ 和 c os θ 就可以了。
如果限制 θ 值在-π/2 < 2θ ≤ π/2，则可令m = - a pq , n = 1( a- app ) , w = sgn( n ) m容易推出： sin 2θ = w , 2sin θ = 2(1 + w ,1 - w2 )cos θ = m 21 - sin2 θ+ n 2利用sin2θ，sin θ和cos θ的值，即得矩阵B 的各元素。
矩阵A 经过旋转变换，选定的非主对角元素a pq 及a qp （一般是绝对值最大的）就被消去，且其主对角元素的平方和增加了 2a2，而非主对角元素的平方和减少了 2a 2 pq ，矩阵元素总的平方和不变。
通过反复选取主元素a pq ， 并作旋转变换，就逐步将矩阵A 变为对角矩阵。
在对称矩阵中共有(n 2-n )/2 个非主对角元素要 被消去, 而每消去一个非主对角元素需要对 2n 个元素进行旋转变换, 对一个元素进行旋转 变换需要 2 次乘法和 1 次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时 间，则消去一个非主对角元素需要 3 个单位运算时间，所以下述算法 21.3 的一轮计算时间 复杂度为 (n 2-n )/2*2n *3=3n 2(n-1)=O(n 3)。
算法 21.3 单处理器上雅可比法求对称矩阵特征值的算法输入：矩阵A n ×n ，ε 输出：矩阵特征值 E igenvalueBegin(1)while (max >ε) do(1.1) max =a[1,2] (1.2)for i =1 to n dofor j = i +1 to n doif (│a [i ,j ]) │>max ) then max =│a [i ,j ] │，p =i ，q =j end if end for end for (1.3)Compute:m =- a [p ,q ]，n=(a [q ,q ]- a [p ,p ])/2，w =sgn(n )*m/sqrt(m *m +n *n )， si 2=w ，si 1=w /sqrt(2(1+ sqrt(1-w *w )))，co1=sqrt(1-si 1*si 1)， b [p ,p ]= a [p ,p ]*co1*co1+ a [q ,q ]*si 1*si 1+ a [p ,q ]*si 2， b [q ,q ]= a [p ,p ]*si 1*si 1+ a [q ,q ]*co1*co1- a [p ,q ]*si 2， b [q , p ]=0，b [p ,q ]=0(1.4)for j =1 to n doif ((j ≠ p ) and ( j ≠ q )) then (i)b [p ,j ]= a [p ,j ]*co1+ a [q ,j ]*si 1 (ii)b [q ,j ]= -a [p ,j ]*si 1 + a [q ,j ]*co1 end if end for(1.5)for i =1 to n doEndif ((i ≠ p ) and (i ≠ q )) then(i)b [i , p ]= a [i ,p ]*co1+ a [i ,q ]*si 1 (ii)b [i , q ]= - a [i ,p ]*si 1+ a [i ,q ]*co1 end if end for(1.6)for i=1 to n dofor j=1 to n do a [i,j ]=b [i,j ] end for end for end while (2)for i =1 to n doEigenvalue [i ]= a [i , i ] end for1.2.2 雅可比法求对称矩阵特征值的并行算法串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。