邻接矩阵与阶层计算(完整版)实用资料

邻接矩阵与阶层计算（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）

邻接矩阵与阶层计算邻接矩阵与阶层计算（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）针对下图给出其项目的邻接矩阵，通过计算得到各阶层的调整信息。邻接矩阵如下：A=0100110000000000000000001000000000000000000000010010000000000000000000000000000001000000000000010000000000000000000000010000000000000000000000111100000000000000001000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000000100000000000000000000010000000000000000000001000000000000000000000011011000000000000000000000000000000000000000001000000000000000000000000000000000000000000000010000000000000000000001000000000000000000000010000000000000000000000经计算得到：可达矩阵为：1111111111111111111111011100100011111111111100110010001111111111110001000000000001000000000010011000000000000000000100010000000000000000001000111111111111000000011000000000000000000000100000000000000000000001000000000000000000000010001111111100000000000100111111110000000000001011111111000000000000011111111100000000000000111111110000000000000001000000000000000000000011000000000000000000000100000000000000000000001011000000000000000000011100000000000000000000110000000000000000000001层次信息如下：层：1层：2,5,6层：3,8,10层：4,7,9141211层：11,12,13,14141211层：15222010681532层：16,17,19,2022201068153213层：18,211315层：221521181997421181997416171617#ifndef

_MATRIX_H

#define

_MATRIX_H

class

Matrix

{

private:

int

row;

//

矩阵的行数

int

col;

//

矩阵的列数

int

n;

//

矩阵元素个数

double*

mtx;

//

动态分配用来存放数组的空间

public:

Matrix(int

row=1,

int

col=1);

//

带默认参数的构造函数

Matrix(int

row,

int

col,

double

mtx[]);

//

用数组创建一个矩阵

Matrix(const

Matrix

&obj);

//

copy构造函数

~Matrix()

{

delete[]

this->mtx;

}

void

print()const;

//

格式化输出矩阵

int

getRow()const

{

return

this->row;

}//

访问矩阵行数

int

getCol()const

{

return

this->col;

}//

访问矩阵列数

int

getN()const

{

return

this->n;

}//

访问矩阵元素个数

double*

getMtx()const

{

return

this->mtx;

}//

获取该矩阵的数组

//

用下标访问矩阵元素

double

get(const

int

i,

const

int

j)const;

//

用下标修改矩阵元素值

void

set(const

int

i,

const

int

j,

const

double

e);

//

重载了一些常用操作符，包括

+，-，x，=，负号，正号，

//

A

=

B

Matrix

&operator=

(const

Matrix

&obj);

//

+A

Matrix

operator+

()const

{

return

*this;

}

//

-A

Matrix

operator-

()const;

//

A

+

B

friend

Matrix

operator+

(const

Matrix

&A,

const

Matrix

&B);

//

A

-

B

friend

Matrix

operator-

(const

Matrix

&A,

const

Matrix

&B);

//

A

*

B

两矩阵相乘

friend

Matrix

operator*

(const

Matrix

&A,

const

Matrix

&B);

//

a

*

B

实数与矩阵相乘

friend

Matrix

operator*

(const

double

&a,

const

Matrix

&B);

//

A

的转置

friend

Matrix

trv(const

Matrix

&A);

//

A

的行列式值，采用列主元消去法

//

求行列式须将矩阵化为三角阵,此处为了防止修改原矩阵，采用传值调用

friend

double

det(Matrix

A);

//

A

的逆矩阵，采用高斯-若当列主元消去法

friend

Matrix

inv(Matrix

A);

};

#endif

实现

#include<iostream.h>

#include<math.h>

#include<stdlib.h>

#include<iomanip.h>

#include"matrix.h"

//

带默认参数值的构造函数

//

构造一个row行col列的零矩阵

Matrix::Matrix(int

row,

int

col){

this->row

=

row;

this->col

=

col;

this->n

=

row

*

col;

this->mtx

=

new

double[n];

for(int

i=0;

i<n;

i++)

this->mtx[i]

=

0.0;

}

//

用一个数组初始化矩阵

Matrix::Matrix(int

row,

int

col,

double

mtx[])

{

this->row

=

row;

this->col

=

col;

this->n

=

row

*

col;

this->mtx

=

new

double[n];

for(int

i=0;

i<n;

i++)

this->mtx[i]

=

mtx[i];

}

//

拷贝构造函数，因为成员变量含有动态空间，防止传递参数

//

等操作发生错误

Matrix::Matrix(const

Matrix

&obj){

this->row

=

obj.getRow();

this->col

=

obj.getCol();

this->n

=

obj.getN();

this->mtx

=

new

double[n];

for(int

i=0;

i<n;

i++)

this->mtx[i]

=

obj.getMtx()[i];

}

//

格式化输出矩阵所有元素

void

Matrix::print()const{

for(int

i=0;

i<this->row;

i++){

for(int

j=0;

j<this->col;

j++)

if(fabs(this->get(i,j))

<=

1.0e-10)

cout

<<

setiosflags(ios::left)

<<

setw(12)

<<

0.0

<<

'

';

else

cout

<<

setiosflags(ios::left)

<<

setw(12)

<<

this->get(i,j)

<<

'

';

cout

<<endl;

}

}

//

获取矩阵元素

//

注意这里矩阵下标从(0,0)开始

double

Matrix::get(const

int

i,

const

int

j)const{

return

this->mtx[i*this->col

+

j];

}

//

修改矩阵元素

void

Matrix::set(const

int

i,

const

int

j,

const

double

e){

this->mtx[i*this->col

+

j]

=

e;

}

//

重载赋值操作符，由于成员变量中含有动态分配

Matrix

&Matrix::operator=

(const

Matrix

&obj){

if(this

==

&obj)

//

将一个矩阵赋给它自己时简单做返回即可

return

*this;

delete[]

this->mtx;

//

首先删除目的操作数的动态空间

this->row

=

obj.getRow();

this->col

=

obj.getCol();

this->n

=

obj.getN();

this->mtx

=

new

double[n];

//

重新分配空间保存obj的数组

for(int

i=0;

i<n;

i++)

this->mtx[i]

=

obj.getMtx()[i];

return

*this;

}

//

负号操作符，返回值为该矩阵的负矩阵，原矩阵不变

Matrix

Matrix::operator-

()const{

//

为了不改变原来的矩阵，此处从新构造一个矩阵

Matrix

_A(this->row,

this->col);

for(int

i=0;

i<_A.n;

i++)

_A.mtx[i]

=

-(this->mtx[i]);

return

_A;

}

//

矩阵求和，对应位置元素相加

Matrix

operator+

(const

Matrix

&A,

const

Matrix

&B){

Matrix

AB(A.row,

A.col);

if(A.row!=B.row

||

A.col!=B.col){

cout

<<

"Can't

do

A+B\n";

//

如果矩阵A和B行列数不一致则不可相加

exit(0);

}

for(int

i=0;

i<AB.n;

i++)

AB.mtx[i]

=

A.mtx[i]

+

B.mtx[i];

return

AB;

}

//

矩阵减法，用加上一个负矩阵来实现

Matrix

operator-

(const

Matrix

&A,

const

Matrix

&B){

return

(A

+

(-B));

}

//

矩阵乘法

Matrix

operator*

(const

Matrix

&A,

const

Matrix

&B){

if(A.col

!=

B.row){

//

A的列数必须和B的行数一致

cout

<<

"Can't

multiply\n";

exit(0);

}

Matrix

AB(A.row,

B.col);

//

AB用于保存乘积

for(int

i=0;

i<AB.row;

i++)

for(int

j=0;

j<AB.col;

j++)

for(int

k=0;

k<A.col;

k++)

AB.set(i,

j,

AB.get(i,j)

+

A.get(i,k)*B.get(k,j));

return

AB;

}

//

矩阵与实数相乘

Matrix

operator*

(const

double

&a,

const

Matrix

&B){

Matrix

aB(B);

for(int

i=0;

i<aB.row;

i++)

for(int

j=0;

j<aB.col;

j++)

aB.set(i,j,

a*B.get(i,j));

return

aB;

}

//

矩阵的转置

将(i,j)与(j,i)互换

//

此函数返回一个矩阵的转置矩阵，并不改变原来的矩阵

Matrix

trv(const

Matrix

&A){

Matrix

AT(A.col,

A.row);

for(int

i=0;

i<AT.row;

i++)

for(int

j=0;

j<AT.col;

j++)

AT.set(i,

j,

A.get(j,i));

return

AT;

}

//

矩阵行列式值，采用列主元消去法

double

det(Matrix

A){

if(A.row

!=

A.col)

{

//

矩阵必须为n*n的才可进行行列式求值

cout

<<

"error"

<<

endl;

return

0.0;

//如果不满足行列数相等返回0.0

}

double

detValue

=

1.0;

//

用于保存行列式值

for(int

i=0;

i<A.getRow()-1;

i++){

//

需要n-1步列化零操作

//------------------

选主元

---------------------------------

double

max

=

fabs(A.get(i,i));

//

主元初始默认为右下方矩阵首个元素

int

ind

=

i;

//

主元行号默认为右下方矩阵首行

for(int

j=i+1;

j<A.getRow();

j++){

//

选择列主元

if(fabs(A.get(j,i))

>

max){

//

遇到绝对值更大的元素

max

=

fabs(A.get(j,i));

//

更新主元值

ind

=

j;

//

更新主元行号

}

}//loop

j

//-------------------

移动主元行

-----------------------------

if(max

<=

1.0e-10)

return

0.0;

//

右下方矩阵首行为零，显然行列式值为零

if(ind

!=

i){//

主元行非右下方矩阵首行

for(int

k=i;

k<A.getRow();

k++){

//

将主元行与右下方矩阵首行互换

double

temp

=

A.get(i,k);

A.set(i,k,A.get(ind,k));

A.set(ind,k,temp);

}

detValue

=

-detValue;

//

互换行列式两行，行列式值反号

}

//-------------------

消元

----------------------------------

for(j=i+1;

j<A.getRow();

j++){

//

遍历行

double

temp

=

A.get(j,i)/A.get(i,i);

for(int

k=i;

k<A.getRow();

k++)

//

遍历行中每个元素，行首置0

A.set(j,

k,

A.get(j,k)-A.get(i,k)*temp);

}

detValue

*=

A.get(i,i);

//

每步消元都会产生一个对角线上元素，将其累乘

}//

loop

i

//

注意矩阵最后一个元素在消元的过程中没有被累乘到

return

detValue

*

A.get(A.getRow()-1,

A.getRow()-1);}//det()

//

A的逆矩阵

高斯-若当消去法，按列选主元

Matrix

inv(Matrix

A){

if(A.row

!=

A.col){

//

只可求狭义逆矩阵，即行列数相同

cout

<<

"Matrix

should

be

N

x

N\n";

exit(0);

}

//

构造一个与A行列相同的单位阵B

Matrix

B(A.row,A.col);

for(int

r=0;

r<A.row;

r++)

for(int

c=0;

c<A.col;

c++)

if(r

==

c)

B.set(r,c,1.0);

//

对矩阵A进行A.row次消元运算，每次保证第K列只有对角线上非零

//

同时以同样的操作施与矩阵B，结果A变为单位阵B为所求逆阵

for(int

k=0;

k<A.row;

k++){

//------------------

选主元

--------------------------------------

double

max

=

fabs(A.get(k,k));

//

主元初始默认为右下方矩阵首个元素

int

ind

=

k;

//

主元行号默认为右下方矩阵首行

//

结果第ind行为列主元行

for(int

n=k+1;

n<A.getRow();

n++){

if(fabs(A.get(n,k))

>

max){//

遇到绝对值更大的元素

max

=

fabs(A.get(n,k));

//

更新主元值

ind

=

n;//

更新主元行号

}

}

//-------------------

移动主元行

--------------------------------

if(ind

!=

k){//

主元行不是右下方矩阵首行

for(int

m=k;

m<A.row;

m++){//

将主元行与右下方矩阵首行互换

double

tempa

=

A.get(k,m);

A.set(k,

m,

A.get(ind,m));

A.set(ind,

m,

tempa);

}

for(m=0;

m<B.row;

m++){

double

tempb

=

B.get(k,m);

//

对矩阵B施以相同操作

B.set(k,

m,

B.get(ind,m));

//

B与A阶数相同，可在一个循环中

B.set(ind,

m,

tempb);

}

}

//---------------------

消元

-----------------------------------

//

第k次消元操作，以第k行作为主元行，将其上下各行的第k列元素化为零

//

同时以同样的参数对B施以同样的操作，此时可以将B看作A矩阵的一部分

for(int

i=0;

i<A.col;

i++){

if(i

!=

k){

double

Mik

=

-A.get(i,k)/A.get(k,k);

for(int

j=k+1;

j<A.row;

j++)

A.set(i,

j,

A.get(i,j)

+

Mik*A.get(k,j));

for(j=0;

j<B.row;

j++)

B.set(i,

j,

B.get(i,j)

+

Mik*B.get(k,j));

}//end

if

}//loop

i

double

Mkk

=

1.0/A.get(k,k);

for(int

j=0;

j<A.row;

j++)

A.set(k,

j,

A.get(k,j)

*

Mkk);

for(j=0;

j<B.row;

j++)

B.set(k,

j,

B.get(k,j)

*

Mkk);

}//loop

k

return

B;

}//inv()9矩阵特征值计算在实际的工程计算中，经常会遇到求n阶方阵A的特征值(Eigenvalue)与特征向量(Eigenvector)的问题。对于一个方阵A，如果数值λ使方程组Ax=λx即(A-λIn)x=0有非零解向量(SolutionVector)x，则称λ为方阵A的特征值，而非零向量x为特征值λ所对应的特征向量，其中In为n阶单位矩阵。由于根据定义直接求矩阵特征值的过程比较复杂，因此在实际计算中，往往采取一些数值方法。本章主要介绍求一般方阵绝对值最大的特征值的乘幂(Power)法、求对称方阵特征值的雅可比法和单侧旋转(One-sideRotation)法以及求一般矩阵全部特征值的QR方法及一些相关的并行算法。求解矩阵最大特征值的乘幂法乘幂法及其串行算法在许多实际问题中，只需要计算绝对值最大的特征值，而并不需要求矩阵的全部特征值。乘幂法是一种求矩阵绝对值最大的特征值的方法。记实方阵A的n个特征值为λii=(1,2,…,n)，且满足：│λ1│≥│λ2│≥│λ3│≥…≥│λn│特征值λi对应的特征向量为xi。乘幂法的做法是：①取n维非零向量v0作为初始向量；②对于k=1,2,…,做如下迭代：uk=Avk-1vk=uk/║uk║∞直至ε为止，这时vk+1就是A的绝对值最大的特征值λ1所对应的特征向量x1。若vk-1与vk的各个分量同号且成比例，则λ1=║uk║∞；若vk-1与vk的各个分量异号且成比例，则λ1=-║uk║∞。若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，则因为一轮计算要做一次矩阵向量相乘、一次求最大元操作和一次规格化操作，所以下述乘幂法串行算法21.1的一轮计算时间为n2+2n=O(n2)。算法21.1单处理器上乘幂法求解矩阵最大特征值的算法输入：系数矩阵An×n，初始向量vn×1，ε输出：最大的特征值maxBegin while(│diff│>ε)do(1)fori=1tondo(1.1)sum=0(1.2)forj=1tondosum=sum+a[i,j]*x[j]endfor(1.3)b[i]=sumendfor(2)max=│b[1]│(3)fori=2tondoif(│b[i]│>max)thenmax=│b[i]│endifendfor(4)fori=1tondox[i]=b[i]/maxendfor(5)diff=max-oldmax,oldmax=maxendwhileEnd乘幂法并行算法乘幂法求矩阵特征值由反复进行矩阵向量相乘来实现，因而可以采用矩阵向量相乘的数据划分方法。设处理器个数为p，对矩阵A按行划分为p块，每块含有连续的m行向量，这里，编号为i的处理器含有A的第im至第(i+1)m-1行数据，(i=0,1,…,p-1)，初始向量v被广播给所有处理器。各处理器并行地对存于局部存储器中a的行块和向量v做乘积操作，并求其m维结果向量中的最大值localmax，然后通过归约操作的求最大值运算得到整个n维结果向量中的最大值max并广播给所有处理器，各处理器利用max进行规格化操作。最后通过扩展收集操作将各处理器中的维结果向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下：算法21.2乘幂法求解矩阵最大特征值的并行算法输入：系数矩阵An×n，初始向量vn×1，ε输出：最大的特征值maxBegin 对所有处理器my_rank(my_rank=0,…,p-1)同时执行如下的算法: while(│diff│>ε)do/*diff为特征向量的各个分量新旧值之差的最大值*/(1)fori=0tom-1do/*对所存的行计算特征向量的相应分量*/(1.1)sum=0(1.2)forj=0ton-1dosum=sum+a[i,j]*x[j]endfor(1.3)b[i]=sumendfor(2)localmax=│b[0]│/*对所计算的特征向量的相应分量求新旧值之差的最大值localmax*/(3)fori=1tom-1doif(│b[i]│>localmax)thenlocalmax=│b[i]│endifendfor(4)用Allreduce操作求出所有处理器中localmax值的最大值max并广播到所有处理器中 (5)fori=0tom-1do/*对所计算的特征向量归一化*/b[i]=b[i]/maxendfor(6)用Allgather操作将各个处理器中计算出的特征向量的分量的新值组合并广播到所有处理器中(7)diff=max-oldmax,oldmax=maxendwhileEnd若各取一次乘法和加法运算时间、一次除法运算时间、一次比较运算时间为一个单位时间，一轮迭代的计算时间为mn+2m；一轮迭代中，各处理器做一次归约操作，通信量为1，一次扩展收集操作，通信量为m，则通信时间为。因此乘幂法的一轮并行计算时间为。MPI源程序请参见所附光盘。求对称矩阵特征值的雅可比法雅可比法求对称矩阵特征值的串行算法若矩阵A=[aij]是n阶实对称矩阵，则存在一个正交矩阵U，使得UTAU=D其中D是一个对角矩阵，即这里λ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,θ)TR(p,q,θ)=In不妨记B=R(p,q,θ)TAR(p,q,θ)，如果将右端展开，则可知矩阵B的元素与矩阵A的元素之间有如下关系：bpp=appcos2θ+aqqsin2θ+apqsin2θbqq=appsin2θ+aqqcos2θ-apqsin2θbpq=((aqq-app)sin2θ)/2+apqcos2θbqp=bpqbpj=apjcosθ+aqjsinθbqj=-apjsinθ+aqjcosθbip=aipcosθ+aiqsinθbiq=-aipsinθ+aiqcosθbij=aij其中1≤i,j≤n且i,j≠p,q。因为A为对称矩阵，R为正交矩阵，所以B亦为对称矩阵。若要求矩阵B的元素bpq=0，则只需令((aqq-app)sin2θ)/2+apqcos2θ=0，即：在实际应用时，考虑到并不需要解出θ，而只需要求出sin2θ，sinθ和cosθ就可以了。如果限制θ值在-π/2<2θ≤π/2，则可令容易推出：,,利用sin2θ，sinθ和cosθ的值，即得矩阵B的各元素。矩阵A经过旋转变换，选定的非主对角元素apq及aqp（一般是绝对值最大的）就被消去，且其主对角元素的平方和增加了，而非主对角元素的平方和减少了2a2pq，矩阵元素总的平方和不变。通过反复选取主元素apq，并作旋转变换，就逐步将矩阵A变为对角矩阵。在对称矩阵中共有(n2-n)/2个非主对角元素要被消去,而每消去一个非主对角元素需要对2n个元素进行旋转变换,对一个元素进行旋转变换需要2次乘法和1次加法，若各取一次乘法运算时间或一次加法运算时间为一个单位时间，则消去一个非主对角元素需要3个单位运算时间，所以下述算法21.3的一轮计算时间复杂度为(n2-n)/2*2n*3=3n2(n-1)=O(n3)。算法21.3单处理器上雅可比法求对称矩阵特征值的算法输入：矩阵An×n，ε输出：矩阵特征值EigenvalueBegin (1)while(max>ε)do (1.1)max=a[1,2] (1.2)fori=1tondo forj=i+1tondo if(│a[i,j])│>max)thenmax=│a[i,j]│，p=i，q=jendif endfor endfor(1.3)Compute:m=-a[p,q]，n=(a[q,q]-a[p,p])/2，w=sgn(n)*m/sqrt(m*m+n*n)，si2=w，si1=w/sqrt(2(1+sqrt(1-w*w)))，co1=sqrt(1-si1*si1)，b[p,p]=a[p,p]*co1*co1+a[q,q]*si1*si1+a[p,q]*si2，b[q,q]=a[p,p]*si1*si1+a[q,q]*co1*co1-a[p,q]*si2，b[q,p]=0，b[p,q]=0(1.4)forj=1tondoif((j≠p)and(j≠q))then (i)b[p,j]=a[p,j]*co1+a[q,j]*si1 (ii)b[q,j]=-a[p,j]*si1+a[q,j]*co1endifendfor(1.5)fori=1tondoif((i≠p)and(i≠q))then(i)b[i,p]=a[i,p]*co1+a[i,q]*si1 (ii)b[i,q]=-a[i,p]*si1+a[i,q]*co1endifendfor(1.6)fori=1tondoforj=1tondoa[i,j]=b[i,j]endforendforendwhile(2)fori=1tondo Eigenvalue[i]=a[i,i]endforEnd雅可比法求对称矩阵特征值的并行算法串行雅可比算法逐次寻找非主对角元绝对值最大的元素的方法并不适合于并行计算。因此，在并行处理中，我们每次并不寻找绝对值最大的非主对角元消去，而是按一定的顺序将A中的所有上三角元素全部消去一遍，这样的过程称为一轮。由于对称性，在一轮中，A的下三角元素也同时被消去一遍。经过若干轮，可使A的非主对角元的绝对值减少，收敛于一个对角方阵。具体算法如下：设处理器个数为p，对矩阵A按行划分为2p块，每块含有连续的m/2行向量，记，这些行块依次记为A0,A1,…,A2p-1，并将A2i与A2i+1存放与标号为i的处理器中。每轮计算开始，各处理器首先对其局部存储器中所存的两个行块的所有行两两配对进行旋转变换，消去相应的非对角线元素。然后按图21.1所示规律将数据块在不同处理器之间传送，以消去其它非主对角元素。图STYLEREF1\s21.SEQ图\*ARABIC\s11p=4时的雅可比算法求对称矩阵特征值的数据交换图这里长方形框中两个方格内的整数被看成是所移动的行块的编号。在要构成新的行块配对时，只要将每个处理器所对应的行块按箭头方向移至相邻的处理器即可，这样的传送可以在行块之间实现完全配对。当编号为i和j的两个行块被送至同一处理器时，令编号为i的行块中的每行元素和编号为j的行块中的每行元素配对，以消去相应的非主对角元素，这样在所有的行块都两两配对之后，可以将所有的非主对角元素消去一遍，从而完成一轮计算。由REF_Ref37669518\h图21.1可以看出，在一轮计算中，处理器之间要2p-1次交换行块。为保证不同行块配对时可以将原矩阵A的非主对角元素消去，引入变量b记录每个处理器中的行块数据在原矩阵A中的实际行号。在行块交换时，变量b也跟着相应的行块在各处理器之间传送。处理器之间的数据块交换存在如下规律：0号处理器前一个行块(简称前数据块，后一个行块简称后数据块)始终保持不变，将后数据块发送给1号处理器，作为1号处理器的前数据块。同时接收1号处理器发送的后数据块作为自己的后数据块。p-1号处理器首先将后数据块发送给编号为p-2的处理器，作为p-2号处理器的后数据块。然后将前数据块移至后数据块的位置上，最后接收p-2号处理器发送的前数据块作为自己的前数据块。编号为my_rank的其余处理器将前数据块发送给编号为my_rank+1的处理器，作为my_rank+1号处理器的前数据块。将后数据块发送给编号为my_rank-1的处理器，作为my_rank-1号处理器的后数据块。为了避免在通信过程中发生死锁，奇数号处理器和偶数号处理器的收发顺序被错开。使偶数号处理器先发送后接收；而奇数号处理器先将数据保存在缓冲区中，然后接收偶数号处理器发送的数据，最后再将缓冲区中的数据发送给偶数号处理器。数据块传送的具体过程描述如下：算法21.4雅可比法求对称矩阵特征值的并行过程中处理器之间的数据块交换算法输入：矩阵A的行数据块和向量b的数据块分布于各个处理器中输出：在处理器阵列中传送后的矩阵A的行数据块和向量b的数据块Begin 对所有处理器my_rank(my_rank=0,…,p-1)同时执行如下的算法: /*矩阵A和向量b为要传送的数据块*/ (1)if(my-rank=0)then/*0号处理器*/ (1.1)将后数据块发送给1号处理器 (1.2)接收1号处理器发送来的后数据块作为自己的后数据块endif(2)if((my-rank=p-1)and(my-rankmod2≠0))then/*处理器p-1且其为奇数*/ (2.1)fori=m/2tom-1do/*将后数据块保存在缓冲区buffer中*/ forj=0ton-1do buffer[i-m/2,j]=a[i,j] endfor endfor (2.2)fori=m/2tom-1do buf[i-m/2]=b[i] endfor (2.3)fori=0tom/2-1do/*将前数据块移至后数据块的位置上*/ forj=0ton-1do a[i+m/2,j]=a[i,j] endfor endfor (2.4)fori=0tom/2-1do b[i+m/2]=b[i]endfor (2.5)接收p-2号处理器发送的数据块作为自己的前数据块 (2.6)将buffer中的后数据块发送给编号为p-2的处理器endif(3)if((my-rank=p-1)and(my-rankmod2=0))then/*处理器p-1且其为偶数*/ (3.1)将后数据块发送给编号为p-2的处理器 (3.2)fori=0tom/2-1do/*将前数据块移至后数据块的位置上*/ forj=0ton-1do a[i+m/2,j]=a[i,j] endfor endfor (3.3)fori=0tom/2-1do b[i+m/2]=b[i] endfor (3.4)接收p-2号处理器发送的数据块作为自己的前数据块endif(4)if((my-rank≠p-1)and(my-rank≠0))then/*其它的处理器*/ (4.1)if(my-rankmod2=0)then/*偶数号处理器*/(i)将前数据块发送给编号为my_rank+1的处理器(ii)将后数据块发送给编号为my_rank-1的处理器(ii)接收编号为my_rank-1的处理器发送的数据块作为自己的前数据块(iv)接收编号为my_rank+1的处理器发送的数据块作为自己的后数据块else/*奇数号处理器*/(v)fori=0tom-1do/*将前后数据块保存在缓冲区buffer中*/ forj=0ton-1do buffer[i,j]=a[i,j] endforendfor(vi)fori=0tom-1do buf[i]=b[i]endfor(vii)接收编号为my_rank-1的处理器发送的数据块作为自己的前数据块(viii)接收编号为my_rank+1的处理器发送的数据块作为自己的后数据块(ix)将存于buffer中的前数据块发送给编号为my_rank+1的处理器(x)将存于buffer中的后数据块发送给编号为my_rank-1的处理器endifendifEnd各处理器并行地对其局部存储器中的非主对角元素aij进行消去，首先计算旋转参数并对第i行和第j行两行元素进行旋转行变换。然后通过扩展收集操作将相应的旋转参数及第i列和第j列的列号按处理器编号连接起来并广播给所有处理器。各处理器在收到这些旋转参数和列号之后，按0,1,…,p-1的顺序依次读取旋转参数及列号并对其m行中的第i列和第j列元素进行旋转列变换。经过一轮计算的2p-1次数据交换之后，原矩阵A的所有非主对角元素都被消去一次。此时，各处理器求其局部存储器中的非主对角元素的最大元localmax，然后通过归约操作的求最大值运算求得将整个n阶矩阵非主对角元素的最大元max，并广播给所有处理器以决定是否进行下一轮迭代。具体算法框架描述如下：算法21.5雅可比法求对称矩阵特征值的并行算法输入：矩阵An×n，ε输出：矩阵A的主对角元素即为A的特征值Begin 对所有处理器my_rank(my_rank=0,…,p-1)同时执行如下的算法: (a)fori=0tom-1do b[i]=myrank*m+i/*b记录处理器中的行块数据在原矩阵A中的实际行号*/endfor (b)while(│max│>ε)do/*max为A中所有非对角元最大的绝对值*/ (1)fori=my_rank*mto(my_rank+1)*m-2do/*对本处理器内部所有行两两配对进行旋转变换*/ forj=i+1to(my_rank+1)*m-1do(1.1)r=imodm,t=jmodm/*i,j为进行旋转变换行(称为主行)的实际行号,r,t为它们在块内的相对行号*/ (1.2)if(a[r,j]≠0)then/*对四个主元素的旋转变换*/ (i)Compute:f=-a[r,j],g=(a[t,j]-a[r,i])/2,h=sgn(g)*f/sqrt(f*f+g*g),sin2=h,sin1=h/sqrt(2*(1+sqrt(1-h*h))),cos1=sqrt(1-sin1*sin1)，bpp=a[r,i]*cos1*cos1+a[t,j]*sin1*sin1+a[r,j]*sin2，bqq=a[r,i]*sin1*sin1+a[t,j]*cos1*cos1-a[r,j]*sin2，bpq=0,bqp=0 (ii)forv=0ton-1do/*对两个主行其余元素的旋转变换*/ if((v≠i)and(v≠j))then br[v]=a[r,v]*cos1+a[t,v]*sin1 a[t,v]=-a[r,v]*sin1+a[t,v]*cos1 endifendfor(iii)forv=0ton-1do if((v≠i)and(v≠j))then a[r,v]=br[v] endifendfor (iv)forv=0tom-1do/*对两个主列在本处理器内的其余元素的旋转变换*/if((v≠r)and(v≠t))then bi[v]=a[v,i]*cos1+a[v,j]*sin1 a[v,j]=-a[v,i]*sin1+a[v,j]*cos1endifendfor (v)forv=0tom-1doif((v≠r)and(v≠t))then a[v,i]=bi[v]endifendfor (vi)Compute: a[r,i]=bpp,a[t,j]=bqq,a[r,j]=bpq,a[t,i]=bqp，/*用temp1保存本处理器主行的行号和旋转参数*/temp1[0]=sin1,temp1[1]=cos1，temp1[2]=(float)i,temp1[3]=(float)jelse(vii)Compute: temp1[0]=0,temp1[1]=0，temp1[2]=0,temp1[3]=0endif (1.3)将所有处理器temp1中的旋转参数及主行的行号 按处理器编号连接起来并广播给所有处理器,存于temp2中 (1.4)current=0 (1.5) forv=1topdo/*根据temp2中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/ (i)Compute: s1=temp2[(v-1)*4+0],c1=temp2[(v-1)*4+1]， i1=(int)temp2[(v-1)*4+2],j1=(int)temp2[(v-1)*4+3] (ii)if(s1、c1、i1、j1中有一不为0)then if(my-rank≠current)then forz=0tom-1do zi[z]=a[z,i1]*c1+a[z,j1]*s1 a[z,j1]=-a[z,i1]*s1+a[z,j1]*c1 endfor forz=0tom-1do a[z,i1]=zi[z] endfor endifendif (iii)current=current+1endfor endforendfor (2)forcounter=1to2p-2do/*进行2p-2次处理器间的数据交换,并对交换后处理器中所有行两两配对进行旋转变换*/ (2.1)Data_exchange()/*处理器间的数据交换*/ (2.2)fori=0tom/2-1doforj=m/2tom-1do (i)if(a[i,b[j]]≠0)then/*对四个主元素的旋转变换*/ ①Compute:f=-a[i,b[j]]，g=(a[j,b[j]]-a[i,b[i]])/2，h=sgn(g)*f/sqrt(f*f+g*g)，sin2=h,sin1=h/sqrt(2*(1+sqrt(1-h*h)))，cos1=sqrt(1-sin1*sin1)，bpp=a[i,b[i]]*cos1*cos1+a[j,b[j]]*sin1*sin1+a[i,b[j]]*sin2，bqq=a[i,b[i]]*sin1*sin1+a[j,b[j]]*cos1*cos1-a[i,b[j]]*sin2，bpq=0,bqp=0 ②forv=0ton-1do/*对两个主行其余元素的旋转变换*/ if((v≠b[i])and(v≠b[j]))then br[v]=a[i,v]*cos1+a[j,v]*sin1 a[j,v]=-a[i,v]*sin1+a[j,v]*cos1 endifendfor ③forv=0ton-1do if((v≠b[i])and(v≠b[j]))then a[i,v]=br[v] endifendfor ④forv=0tom-1do/*对本处理器内两个主列的其余元素旋转变换*/ if((v≠i)and(v≠j))thenbi[v]=a[v,b[i]]*cos1+a[v,b[j]]*sin1a[v,b[j]]=-a[v,b[i]]*sin1+a[v,b[j]]*cos1 endifendfor ⑤forv=0tom-1do if((v≠i)and(v≠j))then a[v,b[i]]=bi[v] endifendfor ⑥Compute:a[i,b[i]]=bpp,a[j,b[j]]=bqq,a[i,b[j]]=bpq,a[j,b[i]]=bqp/*用temp1保存本处理器主行的行号和旋转参数*/temp1[0]=sin1,temp1[1]=cos1,temp1[2]=(float)b[i],temp1[3]=(float)b[j] else ⑦Compute:temp1[0]=0,temp1[1]=0,temp1[2]=0,temp1[3]=0 endif(ii)将所有处理器temp1中的旋转参数及主行的行号按处理器编号连接起来并广播给所有处理器,存于temp2中 (iii)current=0 (iv)forv=1topdo/*根据temp2中的其它处理器的旋转参数及主行的行号对相关的列在本处理器的部分进行旋转变换*/ ①Compute: s1=temp2[(v-1)*4+0],c1=temp2[(v-1)*4+1], i1=(int)temp2[(v-1)*4+2],j1=(int)temp2[(v-1)*4+3] ②if(s1、c1、i1、j1中有一不为0)then if(my-rank≠current)then forz=0tom-1do zi[z]=a[z,i1]*c1+a[z,j1]*s1 a[z,j1]=-a[z,i1]s1+a[z,j1]*c1 endfor forz=0tom-1do a[z,i1]=zi[z] endfor endif endif③current=current+1endforendforendforendfor (3)Data-exchange()/*进行一轮中的最后一次处理器间的数据交换,使数据回到原来的位置*/ (4)localmax=max/*localmax为本处理器中非对角元最大的绝对值*/ (5)fori=0tom-1do forj=0ton-1do if((m*my-rank+i)≠j)then if(│a[i,j]│>localmax)thenlocalmax=│a[i,j]│endif endif endforendfor (6)通过Allreduce操作求出所有处理器中localmax的最大值为新的max值endwhileEnd在上述算法中,每个处理器在一轮迭代中要处理2p个行块对,由于每一个行块对含有m/2行,因而对每一个行块对的处理要有(m/2)2=m2/4个行的配对,即消去m2/4个非主对角元素.每消去一个非主对角元素要对同行n个元素和同列n个元素进行旋转变换.由于按行划分数据块,对同行n个元素进行旋转变换的过程是在本处理器中进行的.而对同列n个元素进行旋转变换的过程则分布在所有处理器中进行.但由于所有处理器同时在为其它处理器的元素消去过程进行列的旋转变换，对每个处理器而言，对列元素进行旋转变换的处理总量仍然是n个元素。因此，一个处理器每消去一个非主对角元素共要对2n个元素进行旋转变换。而对一个元素进行旋转变换需要2次乘法和1次加法，若取一次乘法运算时间或一次加法运算时间为一个单位时间，则其需要3个单位运算时间，即一轮迭代的计算时间为T1＝3×2p×2n×m2/4＝3n3/p；在每轮迭代中，各个处理器之间以点对点的通信方式相互错开交换数据2p-1次，通信量为mn+m，扩展收集操作n(n-1)/(2p)次，通信量为4，另外有归约操作1次，通信量为1，从而不难得出上述算法求解过程中的总通信时间为：因此雅可比算法求对矩阵特征值的一轮并行计算时间为。我们的大量实验结果说明，对于n阶方阵，用上述算法进行并行计算，一般需要不超过O(logn)轮就可以收敛。 MPI源程序请参见章末附录。求对称矩阵特征值的单侧旋转法单侧旋转法的算法描述求解对称方阵特征值及特征向量的雅可比算法在每次消去计算前都要对非主对角元素选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变换，这称之为双侧旋转(Two-sideRotation)变换。在双侧旋转变换中，方阵的行、列方向都有数据相关关系，使得整个计算中的相关关系复杂，特别是在对高阶矩阵进行特征值求解时，增加了处理器间的通信开销。针对传统雅可比算法的缺点，可以将双侧旋转改为单侧旋转，得出一种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换，使得数据相关关系仅在同列之间，因此方便数据划分，适合并行计算，称为单侧旋转法（One-sideRotation）。若A为一对称方阵，λ是A的特征值，x是A的特征向量，则有：Ax=λx。又A=AT，所以ATAx=Aλx=λλx，这说明了λ2是ATA的特征值，x是ATA的特征向量，即ATA的特征值是A的特征值的平方，并且它们的特征向量相同。我们使用节中所介绍的Givens旋转变换对A进行一系列的列变换，得到方阵Q使其各列两两正交，即AV＝Q，这里V为表示正交化过程的变换方阵。因QTQ=(AV)TAV=VTATAV=diag(δ1,δ2,…,δn)为n阶对角方阵，可见这里δ1,δ2,…,δn是矩阵ATA的特征值，V的各列是ATA的特征向量。由此可得：(i=1,2,…,n)。设Q的第i列为qi,则qiTqi=δi,||qi||2=。因此在将A进行列变换得到各列两两正交的方阵Q后，其各列的谱范数||qi||2即为A的特征值的绝对值。记V的第i列为vi，AV＝Q，所以Avi=qi。又因为vi为A的特征向量，所以Avi=λivi，即推得λivi=qi。因此λi的符号可由向量qi及vi的对应分量是否同号判别，实际上在具体算法中我们只要判断它们的第一个分量是否同号即可。若相同，则λi取正值，否则取负值。求对称矩阵特征值的单侧旋转法的串行算法如下：算法21.6求对称矩阵特征值的单侧旋转法输入：对称矩阵A，精度e输出：矩阵A的特征值存于向量B中Begin(1)while(p>ε)p=0fori=1tondoforj=i+1tondo(1.1)sum[0]=0,sum[1]=0,sum[2]=0(1.2)fork=1tondosum[0]=sum[0]+a[k,i]*a[k,j]sum[1]=sum[1]+a[k,i]*a[k,i]sum[2]=sum[2]+a[k,j]*a[k,j]endfor(1.3)if(│sum[0]│>ε)then(i)aa=2*sum[0]bb=sum[1]-sum[2]rr=sqrt(aa*aa+bb*bb)(ii)if(bb≥0)thenc=sqrt((bb+rr)/(2*rr))s=aa/(2*rr*c)elses=sqrt((rr-bb)/(2*rr))c=aa/(2*rr*s)endif(iii)fork=1tondotemp[k]=c*a[k,i]+s*a[k,j]a[k,j]=[-s]*a[k,i]+c*a[k,j] endfor(iv)fork=1tondotemp1[k]=c*e[k,i]+s*e[k,j]e[k,j]=[-s]*e[k,i]+c*e[k,j]endfor(v)fork=1tondoa[k,i]=temp[k]endfor(vi)fork=1tondoe[k,i]=temp1[k]endfor(vii)if(│sum[0]│>p)thenp=│sum[0]│endifendifendforendforendwhile(2)fori=1tondosu=