北航并行计算矩阵转置作业
北航 数值分析第二次大作业(带双步位移的QR方法)
一、算法设计方案:按题目要求,本程序运用带双步位移的QR方法求解给定矩阵的特征值,并对每一实特征值,求解其相应的特征向量。
总体思路:1)初始化矩阵首先需要将需要求解的矩阵输入程序。
为了防止矩阵在后面的计算中被破坏保存A[][]。
2)对给定的矩阵进行拟上三角化为了尽量减少计算量,提高程序的运行效率,在对矩阵进行QR分解之前,先进行拟上三角化。
由于矩阵的QR 分解不改变矩阵的结构,所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。
这里用函数void QuasiTriangularization()来实现,函数形参为double型N维方阵double a[][N]。
3)对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。
用子程序void QR_decomposition()来实现,将Q、R设为形参,然后将计算出来的结果传入Q和R,然后求出RQ乘积。
4)对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛,对QR分解引入双步位移,适当选取位移量,可以避免进行复数运算。
为了进一步减少计算量,在每次进行QR分解之前,先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。
若可以,则得到特征值,同时对矩阵进行降阶处理;若不可以,则进行QR分解。
这里用函数intTwoStepDisplacement_QR()来实现。
这是用来存储计算得到的特征值的二维数组。
考虑到特征值可能为复数,因此将所有特征值均当成复数处理。
此函数中,QR分解部分用子函数void QR_decompositionMk()实现。
这里形参有三个,分别用来传递引入双步位移后的Mk阵,A矩阵,以及当前目标矩阵的维数m。
5)计算特征向量得到特征值后,计算实特征值相应的特征向量。
这里判断所得特征值的虚数部分是否为零。
求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。
因此先初始化矩阵Array,计算(A-λI),再求解方程组。
北航数值分析计算实习题目二 矩阵QR分解
数值分析实习二院(系)名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值,并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明:1、求矩阵特征值时,要求迭代的精度水平为1210ε-=。
2、打印以下内容:算法的设计方案;全部源程序(要求注明主程序和每个子程序的功能); 矩阵A 经过拟上三角话之后所得的矩阵()1n A -;对矩阵()1n A-进行QR 分解方法结束后所得的矩阵;矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=,和A 的相应于实特征值的特征向量;其中()(),.i e i m i R R I I λλ==如果i λ是实数,则令0.i I =3、采用e 型输出数据,并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁,自定义类Xmatrix ,其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数,即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销,使用引用(&)代替值传递,以节省内存空间,避免溢出.(1)此程序的主要部分为Xmatrix 中的doubleQR()方法,具体如下:Step1:使用矩阵拟上三角化的算法将A 化为拟上三角阵A (n-1)(此处调用Xmatrix 中的preQR()方法)Step2:令121,,10k m n ε-===, 其中k 为迭代次数。
Step3:如果,1m m a ε-≤,则得到A 的一个特征值,m m a ,令1m m =-,goto Step4;否则goto Step5.Step4: 如果1m =,则得到A 的一个特征值11a ,goto Step11;如果0m =,则goto Step11;如果1m >,则goto Step3;Step5(Step6):如果2m =,则得到A 的两个特征值12s s 和(12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
矩阵快速转置算法详解
矩阵快速转置算法详解好啦,今天我们来聊聊矩阵快速转置算法,这可是个让人觉得又爱又恨的话题。
矩阵嘛,想必大家都不陌生,毕竟它就像我们生活中的一张大表格,记录着各种数据。
可这转置,听起来就像是魔法一样,把行变成列,列又变回行,这其中可有不少讲究。
想象一下,一个矩阵就像是你的冰箱,里面装着各种各样的食材。
平常你打开冰箱,看到一层层整齐的食物,心里就觉得舒服。
不过,有时候你也想换个角度,看看这些食材在不同位置的组合。
这个时候,转置就登场啦!转置后的矩阵,就像是你重新整理过的冰箱,虽然东西没变,但看上去就是更有条理,更让人一目了然。
转置到底是怎么个操作呢?其实就是把矩阵里的元素交换位置。
比如,原本在第一行第二列的元素,现在要跑到第二行第一列。
听起来简单对吧?但当你的矩阵大得像个足球场,里面的元素多得跟星星一样,那可就得小心翼翼,别让它们跑丢了。
这里就需要一个“快速”的方法,避免慢吞吞的操作,简直就是搬家时的“火速搬迁”,让你几乎感觉不到麻烦。
说到快速转置算法,最经典的就是“分块转置”。
这就像你在打扫房间的时候,通常不会一次性把所有东西都搬出去,而是先把几个小块整理好。
这个算法也是如此,把大矩阵分成小块,分别转置,然后再合并。
这种方法就像是在料理过程中,把大菜分成几道小菜,最后汇聚成一桌美食,既高效又省力。
分块的方法还有个好处,就是它能更好地利用缓存,提高处理速度,简直是个聪明的“小调皮”。
转置过程中,你可能会碰到一些小坑,比如说对称矩阵。
它们就像是一对双胞胎,无论怎么交换,最后看上去都是一样的。
对此,快速转置算法能巧妙地利用这一点,不用额外的空间,就能轻松搞定,简直让人拍手叫好。
就像生活中遇到的那些小聪明,总是能在关键时刻帮你一把,真是妙不可言。
转置算法的实现可不仅仅是写代码那么简单。
想象一下,代码里的每一行都像是一位舞者,它们需要配合默契,才能让整个矩阵的转置变得流畅自如。
这时候,算法的优化就显得尤为重要。
北航数值分析大作业第二题(fortran)
!计算A(r+1) DO I=1,N DO J=1,N A(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J) ENDDO ENDDO ENDIF ENDDO RETURN END
!***************符号函数子程序*****************! FUNCTION SGN(X) REAL(8) X IF(X>0) THEN SGN=1 ELSE IF(X<0) THEN SGN=-1 ELSE IF(X==0) THEN SGN=0 ENDIF END
DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数 组,1为实部,为虚部 REAL(8) A,A1,A2,C,Q,R,CM E=1E-12 L=1000 !精度水平 !迭代最大次数
OPEN(1,FILE='数值分析大作业第二题计算结果.TXT') DO I=1,N DO J=1,N IF(I==J) THEN A(I,J)=1.52*COS(I+1.2*J) ELSE A(I,J)=SIN(0.5*I+0.2*J) ENDIF ENDDO ENDDO A1=A WRITE(*,"('矩阵A为:')") WRITE(1,"('矩阵A为:')") DO I=1,N DO J=1,N WRITE(*,"(2X,E20.13,2X,\)") A(I,J) WRITE(1,"(2X,E20.13,2X,\)") A(I,J) ENDDO WRITE(*,"(' ')") WRITE(1,"(' ')") ENDDO !使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A(n-1) CALL HESSENBERG(A,N) WRITE(*,"('拟上三角化后矩阵A(n-1)为:')") WRITE(1,"('拟上三角化后矩阵A(n-1)为:')") DO I=1,N DO J=1,N WRITE(*,"(2X,E20.13,2X,\)") A(I,J) WRITE(1,"(2X,E20.13,2X,\)") A(I,J) ENDDO WRITE(*,"('')") WRITE(1,"('')") ENDDO !计算对矩阵A(n-1)实行QR方法迭代结束后所得矩阵 A2=A CALL QRD(A2,N,Q,R)
北航矩阵考题A
an 1 2an 2 (a , a , 1 2 nan n
, an ) (分解不唯一) ( 4 分)
a1 2a2
nan trA
A2 (trA) A 0
A100 ( )99 A (a1 nan )99 A.
(2)可知 f ( A) 的谱公式为 f ( A) f (1 )G1 f (2 )G2 ,(且 G1 , G2 同上), 令 f ( x) sin x f (7) sin 7, f (2) sin 2 得 sin A 的谱分解为
4 4 5 4 1 sin A (sin 7)G1 (sin 2)G2 sin 7 1 sin 2 9 9 5 5 5 4
(3) 由 (e A ) {e7 , e 2 } 谱半径 (e A ) e7 由 (esin A ) {esin 7 , e sin 2 } 行列式 det(esin A ) esin 7 sin 2
2
1 1 1 1 四、 (15 分)设 A 1 2 3 , b 2 , 0 0 0 0
1 n 1 n 2
n (A 2 ) 为收敛
( 4 分)
以下的六、七题中只需任选一题:
4
六、(15 分)(1) 设矩阵 A 最小式 m( x) ( x 2) 2 且 f ( A) 收敛,推导 f ( A) 的广谱 计算公式。解(1)由 ( A 2)2 0 与台乐公式 f ( A) 0 可得公式
(3)求谱半径 (e A ) 与行列式 det(esin A ) .(4):求 ln( I A ) ?
北航并行计算矩阵转置作业
矩阵转置并行实现1、算法描述:若处理器个数为P,且它们的编号依次是0,1, 则将n阶矩阵A分成p 个大小为mxm的子块,m=[n/p]o p个子块组成一个子块阵列,记其中第i行第j列的子块为Aij,它含有第(i-l)m+l至第im行中的第(j-l)m+l至第jm列的所有元素。
对每一处理器按行主方式赋以二维下标,记编号为i的处理器的二维下标为(v,u),其中u=imod A/p ,将A的子块存入下标为(v,u)表示的对应处理器中。
转置分为两步进行:第一步,子块转置;第二步,处理器内部转置。
为了避免对应子块交换数据是处理器发生死锁,可令下三角块先向与之对应的上三角子块发送数据,然后从上三角子块接收数据;上三角子块数据先存放在缓冲区buffe「中,然后从与之对应的下三角子块接收数据,最后再将缓冲区中的数据发送给下三角子块,流程图如下所示:2、程序代码:#include "stdio.h"#include "stdlib.h"#include ”mpi.h”#include "math.h"#define E 0.0001#define a(x,y) a[x*m+y]#define b(x,y) b[x*m+y]#define A(x,y) A[x*size+y]#define B(x,y) B[x*size+y]#define intsize sizeof(int)#define floatsize sizeof(float)#define charsize sizeof(char)int size,N;int m;int t;float *A, *B;double starttime;double timel;double time2;int my_rank;int p;MPI_Status status;FILE *fdA;void Environment_Finalize(float *ajloat *b){free(a);free(b);}int main(int argc, char **argv){int ij/kmy^rank,group_size;float *a,*b;int u,v;float temp;MPIJ nit(&argc,&argv);MPI_Comm_size(MPI_COMM_WORLD/&group_size);M P l_Co m m_ra nk(MPI_COM M_W0 RLD,&my_ra nk); p=group_size;if(my_ra nk==0) starttime=MPI_Wtime(); fdA=fopen("dataln.txt,,;,r H);fscanf(fdA z"%d %d", &size, &N);if(size != N){puts(H The in put is error!"); exit(O);}A=(float*)malloc(floatsize*size*size);B=(float*)malloc(floatsize*size*size);for(i = 0; i < size; i ++){for(j = 0; j < size; j ++) fscanf(fdA, ”%化A+i*size+j);}fclose(fdA);MPI_Bcast(&size,l,MPI JN TAMPI_COMM_WORLD); t=(int)sqrt(p);if (t>size)t=size;i#si ze%t!=0)for(;;){t-;if(size%t==0)break;p=t*t;m=size/t;a=(float *)malloc(floatsize*m*m); b=(float *)malloc(floatsize*m*m); if(a==NULL||b==NULL)printf("allocate space fail!11);if (my_ra nk==0) for(i=0;i<m;i++)for(j=0;j<m;j++)a(i,j)=A(i,j);}if (my_rank==O){for(i=l;i<p;i++){v=i/t;u=i%t;for(j=v*m;j<(v+l)*m;j++)for(k=u*m;k<(u+l)*m;k++)b((j%m),(k%m))=A(j,k);MPI_Send(bm*m,l\/IPI_FLOAl;ijMPI_COIVIM_WORLD);}}else if (my_rank<p)MPI_Recv(am*m,MPI_FLOATQmy_rank,MPI_COI\/IM_WORLD&status);timel=MPI_Wtime();if ((my_rank/t)>(my_rank%t)&&my_rank<p){v= my_rank/t;u= my_ra nk%t;MPLSend(a/m*m/MPLFLOAT, (u*t+v),(u*t+v)/MPLCOMM_WORLD);MPLRecvfa^^^PLFLOATfu^+vhmy^ran^MPLCOMM^WORL^&status);}if ((my_rank/t)<(my_rank%t)&&my_rank<p){v= my_rank/t;u= my_ra nk%t;for(i=0;i<m;i++)for(j=0;j<m;j++)b(ij)=a(ij);MPI_Recv(am*rrLMPI_FLOA7;(屮t+v), my_rank,MPI_COMM_WORLD&status);MPLSend(b/m*m/MPLFLOAl;(u*t+v)/(u*t+v)/MPLCOMM_WORLD);for(i=l;i<m;i++)for(j=0;j<i;j++){temp=a(ij);a(i,j)=a(j,i);a(jj)=temp;}if (my_rank==O){for(i=0;i<m;i++)for(j=0;j<m;j++)B(i,j)=a(i,j);}if (my_rank==O){for(i=l;i<p;i++){MPLRecv(a z m*m z MPLFLOATJJ z MPLCOMM_WORLD,&status); V= i/t; u= i%t;for(j=v*m;j<(v+l)*m;j++) for(k=u*m;k<(u+l)*m;k++)B(j,k)=a((j%m),(k%m));}}else if(my_rank<p)MPLSend(a/m*m/MPLFLOAT,0,my-rank/MPLCOMM_WORLD);if (my_rank==O){printf(H lnput of file V'dataln.txtV'Xn"); printf(,,%d\t%d\n,,/ size, size);for(i=0;i<size;i++)for(j=0;j<size;j++) printf(,,%f\t,,,A(i/j)); printfCV);printf(M\nOutput of Matrix AT\n M);for(i=0;i<size;i++){for(j=0;j<size;j++) printf(l,%f\t,,,B(i/j)); printfCV);}}time2=MPI_Wtime();if (my_rank==0){printf(”\n”);printf(n processor number: %d\n H z p);printf(H Whole running time = %f seconds\n,,/time2-starttime);printf(H Distribute data time = %f sec onds\n'\timel-starttime); prin tf(H Parallel computetime = %f sec on ds\n':time2・time:l);}MPI_Barrier(MPLCOMM_WORLD); MPI_Finalize();Environme nt_Fin alize(a,b); return(O);。
北航数值分析大作业第二题
北航数值分析大作业第二题本页仅作为文档封面,使用时可以删除This document is for reference only-rar21year.March数值分析第二次大作业史立峰SY1505327一、 方案(1)利用循环结构将sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值,得到需要变换的矩阵A ;(2)然后,对矩阵A 利用Householder 矩阵进行相似变换,把A 化为上三角矩阵A (n-1)。
对A 拟上三角化,得到拟上三角矩阵A (n-1),具体算法如下: 记A(1)=A ,并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)( +==。
对于2,,2,1-=n r 执行 1. 若()n r r i a r ir,,3,2)( ++=全为零,则令A(r+1) =A(r),转5;否则转2。
2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若 )(,12r rr r r r a c c h +-=3. 令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0 。
4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。
(3)使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值,也是A 的全部特征值,具体算法如下:1. 给定精度水平0>ε和迭代最大次数L 。
北航数值分析报告大作业二
数值分析大作业(二)学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1 算法设计方案首先将矩阵A 进行拟上三角化,把矩阵A 进行QR 分解,计算出RQ 。
要得出矩阵A 的全部特征值,首先对A 进行QR 的双步位移得出特征值。
最后,采用列主元的高斯消元法求解特征向量。
1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构,因此在进行QR 分解前对矩阵A 进行拟上三角化可以大大减少计算机的计算量,提高程序的运行效率。
具体算法如下所示,记A A =)1(,并记)(r A 的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)(ΛΛ+==。
对于2,,2,1-=n r Λ执行 若()n r r i a r ir,,3,2)(Λ++=全为零,则令)()1(r r A A =+,转5;否则转2。
计算()∑+==nri r ir r a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,12r rr r r r a c c h +-=令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0ΛΛ。
计算r r T r r h u A p /)(=r r rr r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(继续。
1.2 A 的QR 分解具体算法如下所示,记)1(1-=n A A ,并记[]nn r ij r a A ⨯=)(,令I Q =1 对于1,,2,1-=n r Λ执行 1.若()n r r i a r ir ,,3,1)(Λ++=全为零,则令r r Q Q =+1r r A A =+1,转5;否则转2。
北航矩阵大作业
(2.10)
矩阵论
������21 =
− ������2 + ������1
|
+ =0 ������2
(2.5)
− 从������11 = ������22 = 0的事实出发可知,当端口 2 端接有������0 =50Ω 时������1 =0,并有 + + − ������2 =0,在这情况下就有������1 = ������1,������2 = ������2。所以在端口 1 上施加电压������1并使用 − 二次分压,就可求出������2 = ������2,他是在端口 2 上跨接在 50Ω 负载电阻上的电
[������] = [ 0 0.707
0.707 ] 0
下面说明如何从阻抗矩阵来确定[������]矩阵,或反之。首先,我们必须嘉定所 有端口的特征阻抗������0������ 是相等的。为方便起见可设������0������ =1。可知端口 n 的总电压 和电流可写为
+ − ������ ������ = ������ ������ + ������ ������ + − + − ������������ = ������������ − ������������ =������ ������ − ������ ������
− − ������3 ,������3
图 1,任意端口微波网络 图 1 中的端口可以是某种形式的传输线或单一波导传播模式的等效传输 线。在第 n 个端口上的一个指定点处,定义一个端平面������������ ,并定义了等效的入
+ + 射波电压和电流(������ ������ ,������������ )。现在,我们给出第 n 个端平面上的总电压和电
北航 矩阵论 习题2.1参考答案
则
T1
A
0
4
1 ;对于 b(2) (4,3)T ,构造 T2 使 T2b(2) b(2) e2
0 3 2
4 / 5 3 / 5 4 1 5 2
T2
3
/
5
4
/
5
,
T2
3
2
0
1
0 1 0
0 4/5 3/5
所以, T
I 0
0 T2
T1
4 3
/ /
5 5
0 0
3/5
3 )T 3
由 a3 (2, 0, 2)T ,有
12
a3
(a3, b1) (b1, b1)
b1
(a3, b2 ) (b2 , b2 )
b2
(2, 0, 2)
14 26
(3,1, 4)T
13 24
(10 , 14 , 4 )T 13 13 13
0T
13
故
k31
7 13
,
k32
1 2
3 26
i 2
1 i
6
3
0
2i 1 6 3
R
b1
b2
1
i 2
1
2
2
i 2
0
b3
0
1 0
i 3 1
0 0
3 6
0
1
2
i
6
2
3
2 0 0
1
30
1
3
2
6
3 3
1
对
P
中对应
Q1 的列向量做单位化得
P
2
1
2
3
北航并行计算矩阵相乘作业
矩阵相乘并行实现1、算法描述:设有如下矩阵相乘:C=A×B其中A,B分别是m×k和k×n矩阵,C是m×n矩阵。
若处理器个数为p,且它们的编号依次是0,1,…,p-1,则设可将矩阵A、B、C分成p个大小为mxm的子块,其中A=(Aij)m×m ,B=(Bij)m×m,和C=(Cij)m×m,其中A¬ij,Bij和Cij是n×n矩阵。
同时假设和。
定义对角块矩阵,则其中,。
利用此关系式,将节点编号从一维映射到二维,数据,,存放在中,可得到下面的在处理机结点上的算法。
该算法数据交量算法流程如下:流程图如下所示:2、程序代码:#include <stdlib.h>#include <string.h>#include <mpi.h>#include <time.h>#include <stdio.h>#include <math.h>float **A, **B, **C;float *a, *b, *c, *tmp_a, *tmp_b;int dg=1000, dl, dl2,p, sp;int my_rank, my_row, my_col;MPI_Status status;int get_index(int row, int col, int sp){return ((row+sp)%sp)*sp + (col+sp)%sp; }void random_A_B(){int i,j;srand((unsigned int)time(NULL));for(i=0; i<dg ; i++)for(j=0; j<dg ; j++){A[i][j] = rand();B[i][j] = rand();C[i][j] = 0.0;}}void scatter_A_B(){int i,j,k,l;int p_imin,p_imax,p_jmin,p_jmax;for(k=0; k<p; k++){p_jmin = (k % sp ) * dl;p_jmax = (k % sp + 1) * dl-1;p_imin = (k - (k % sp))/sp * dl;p_imax = ((k - (k % sp))/sp +1) *dl -1;l = 0;for(i=p_imin; i<=p_imax; i++){for(j=p_jmin; j<=p_jmax; j++){tmp_a[l] = A[i][j];tmp_b[l] = B[i][j];l++;}}if(k==0){memcpy(a, tmp_a, dl2 * sizeof(float));memcpy(b, tmp_b, dl2 * sizeof(float));} else{MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD);MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD);}}}void init_alignment(){MPI_Sendrecv(a, dl2, MPI_FLOAT, get_index(my_row,my_col-my_row,sp), 1,tmp_a, dl2, MPI_FLOAT, get_index(my_row,my_col+my_row,sp), 1,MPI_COMM_WORLD, &status);memcpy(a, tmp_a, dl2 * sizeof(float) );MPI_Sendrecv(b, dl2, MPI_FLOAT, get_index(my_row-my_col,my_col,sp), 1,tmp_b, dl2, MPI_FLOAT, get_index(my_row+my_col,my_col,sp), 1,MPI_COMM_WORLD, &status);memcpy(b, tmp_b, dl2 * sizeof(float) );}void main_shift(){int i,j,k,l;for(l=0; l<sp; l++){for(i=0; i<dl; i++)for(j=0; j<dl; j++)for(k=0; k<dl; k++)c[i*dl+j] += a[i*dl+k]*b[k*dl+j];MPI_Send(a , dl2, MPI_FLOAT, get_index(my_row, my_col-1, sp), 1, MPI_COMM_WORLD);MPI_Recv(a , dl2, MPI_FLOAT, get_index(my_row, my_col+1, sp), 1, MPI_COMM_WORLD, &status);MPI_Send(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD);MPI_Recv(b , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status);}}void collect_C(){int i,j,i2,j2,k;int p_imin,p_imax,p_jmin,p_jmax;for (i=0;i<dl;i++)for(j=0;j<dl;j++)C[i][j]=c[i*dl+j];for (k=1;k<p;k++){MPI_Recv(c, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &status);p_jmin = (k % sp ) *dl;p_jmax = (k % sp + 1) *dl-1;p_imin = (k - (k % sp))/sp *dl;p_imax = ((k - (k % sp))/sp +1) *dl -1;i2=0;for(i=p_imin; i<=p_imax; i++){j2=0;for(j=p_jmin; j<=p_jmax; j++){C[i][j]=c[i2*dl+j2];j2++;}i2++;}}}void print(float **m,char *str){int i,j;printf("%s",str);for(i=0;i<dg;i++){for(j=0;j<dg;j++)printf("%15.0f ",m[i][j]);printf("\n");}printf("\n");}int main(int argc, char *argv[]){int i;MPI_Init(&argc, &argv);MPI_Comm_size(MPI_COMM_WORLD, &p);MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);sp = sqrt(p);if (sp*sp != p){if (my_rank == 0)printf("Number of processors is not a quadratic number!\n");MPI_Finalize();exit(1);}if (argc != 2){if (my_rank == 0)printf("usage: mpirun -np ProcNum cannon MatrixDimension\n");MPI_Finalize();exit(1);}dg = atoi(argv[1]);dl = dg / sp;dl2 = dl * dl;my_col = my_rank % sp ;my_row = (my_rank-my_col) / sp ;a = (float *)malloc( dl2 * sizeof(float) );b = (float *)malloc( dl2 * sizeof(float) );c = (float *)malloc( dl2 * sizeof(float) );for(i=0; i<dl2 ; i++)c[i] = 0.0;tmp_a = (float *)malloc( dl2 * sizeof(float) );tmp_b = (float *)malloc( dl2 * sizeof(float) );if (my_rank == 0){A = (float **)malloc( dg * sizeof(float*) );B = (float **)malloc( dg * sizeof(float*) );C = (float **)malloc( dg * sizeof(float*) );for(i=0; i<dg; i++){A[i] = (float *)malloc( dg * sizeof(float) );B[i] = (float *)malloc( dg * sizeof(float) );C[i] = (float *)malloc( dg * sizeof(float) );}random_A_B();scatter_A_B();} else{MPI_Recv(a, dl2, MPI_FLOAT, 0 , 1, MPI_COMM_WORLD, &status);MPI_Recv(b, dl2, MPI_FLOAT, 0 , 2, MPI_COMM_WORLD, &status);}init_alignment();main_shift();if(my_rank == 0){collect_C();print(A,"random matrix A : \n");print(B,"random matrix B : \n");print(C,"Matrix C = A * B : \n");} else{MPI_Send(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD);}MPI_Barrier(MPI_COMM_WORLD);MPI_Finalize();return 0;}。
北航矩阵理论期末试卷有解析
姓名: 学号:1. (42分)填空(1)设1234=(1,1,1,1),=(1,1,1,1),=(1,1,1,1),=(1,1,1,1)T T T T αααα------是R 4的一组基, 则(1,2,1,1)T 在上述基下的坐标是___________________. (--5111(,,,)4444T ) (2)在三次多项式空间3R[x]中,由多项式组2312()142,()19f x x x x f x x =+-+=-+ 233233432,()56,()5752x x f x x x f x x x x -+=-++=+-+张成的子空间维数是___2___.(3)设矩阵123100A=024B=52000161a ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭,, 当参数a 满足_______ ( 65≠ )时,矩阵A 与B 相似.(4)A = ⎛⎫ ⎪ ⎪ ⎪⎝⎭50203121-7, 则A 的全部盖尔圆为_______________________________,且A 是一个________(可逆或者不可逆)矩阵.(5)设⨯∈nC m n A , 则矩阵A 的正奇异值有______个,_____(是或否)存在矩阵B 使得BA=I n .(6)矩阵幂级数kk ∑∞=⎪⎪⎭⎫ ⎝⎛04.07.05.02.0=_____⎪⎪⎭⎫ ⎝⎛87561310_____________。
(7)设⎪⎪⎪⎭⎫ ⎝⎛--=201034011A ,则A 的Jordan 标准形J=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛200010011100110002或 。
(8)设10022i A i ⎛⎫ ⎪= ⎪ ⎪⎝⎭,则A +=_____10210-2i 10i ⎛⎫ ⎪-⎝⎭___________。
(9)若⨯∈442C ,且A =A ,A 的秩是2,则|A-2I|A =__4__, Sin A 的迹=__2sin1__.(10)设023302230i A i i ⎛⎫⎪= ⎪ ⎪⎝⎭,则||A||1 = _6___, ||A||F___. 2.(15分) 设 A = ⎛⎫⎪ ⎪ ⎪⎝⎭1-13-3-33, 求A 的奇异值分解.解:⎛⎫= ⎪⎝⎭H 13-3-1-33A ,则 ⎛⎫-⎛⎫⎛⎫- ⎪==⎪⎪ ⎪--⎝⎭⎝⎭ ⎪⎝⎭H1113-31919AA3-31-331919-33 ()()λλλλλλ--==--=--221919A A 1919381919T Iλλ==1238,0,对λ=138 ,求⎛⎫⎛⎫= ⎪ ⎪ ⎪⎝⎭⎝⎭12191901919x x 得η⎛⎫⎪⎝⎭11=-1对λ=20 ,求⎛⎫⎛⎫-= ⎪ ⎪ ⎪-⎝⎭⎝⎭12191901919x x 得η⎛⎫⎪⎝⎭21=1分别单位化为; ⎪⎪11-11令 ⎪11=-11V而η⎛⎫⎛⎫⎛⎫ ⎪ ⎪== ⎪ ⎪ ⎪⎝⎭ ⎪ ⎪⎝⎭⎝⎭11-121A 3-36-1-33-6,补充基为⎛⎫⎛⎫⎛⎫⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭10-63,1,1-31-1令 ⎪=0311-311U所以⎫⎪= ⎪ ⎪ ⎪⎝⎭0A 0000H U V⎛⎫⎫⎫⎪⎪⎪==⎪ ⎪ ⎪⎪ ⎪⎪ ⎪⎝⎭⎝⎭1-600311A0000110000-311HU V3.(10分)设nA R n⨯∈并且A是正交矩阵,证明A的每个特征值iλ的模等于1. 课本P51推论2证明:设A,x x xλλ=为属于的特征向量,共轭转置得H T HA,x xλ=所以H T H HA A=,x x x x x xλλ=即2||=1.λ4.(18分)已知A =⎛⎫⎪⎪⎪⎝⎭-112-221-1-2011-2,b =⎛⎫⎪⎪⎪⎝⎭154.(1)求A的满秩分解,并用满秩分解求+A.(2)判断方程组Ax = b是否有解. (3)求Ax = b的极小范数解或极小最小二乘解.解:(1)101001120000A-⎡⎤⎢⎥−−→-⎢⎥⎢⎥⎣⎦行11101021=011201A FG-⎡⎤-⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦⎢⎥⎣⎦(2)()115145111363514T T TF F F F F--+--⎡⎤⎡⎤===⎢⎥⎢⎥⎣⎦⎣⎦()11612112116511124T T TG G GG G--+⎡⎤⎢⎥-⎡⎤⎢⎥===⎢⎥⎢⎥--⎣⎦⎢⎥--⎣⎦1833181191262210154162218A G F +++--⎡⎤⎢⎥⎢⎥==⎢⎥-⎢⎥---⎣⎦(3)66=33,55AA b b Ax b +⎛⎫⎪≠= ⎪ ⎪⎝⎭故无解.(4)()b 1,T A +=-,9,10,-185.(15分)设 210420101A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,求t A e .解:2||(1)I A -=-λλλ,因为()0A A I -≠所以最小多项式为2(1)-λλ, 设2()c ,()t P a b f e =++=λλλλλ.有:(0)(0)11(0)(0)(1)(1)1t t f P a a tf P t b b t f P e a b c c e t ===⎧⎧⎧⎪⎪⎪''=⇒=⇒=⎨⎨⎨⎪⎪⎪==++=--⎩⎩⎩A 2120=4210211t t t t t e aI bA cA t t e t e t e λ-⎡⎤⎢⎥++=-+⎢⎥⎢⎥-++--⎣⎦。
北航_误差分析和数据处理_第六章作业答案
6—2解:由题意可写出误差方程为()()()⎪⎩⎪⎨⎧--=--=+-=y x y x y x 329.129.039.2321ννν 设有列向量⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=321,ˆ,9.19.09.2νννV y x X L矩阵⎢⎢⎢⎣⎡=213A ⎥⎥⎥⎦⎤--321则由误差方程可列出正规方程为()L A XA AT T=ˆ即为⎩⎨⎧-=+-=-6.41454.13514y x y x 按矩阵求解,则有⎢⎣⎡-=514C ⎥⎦⎤-145⎢⎢⎢⎢⎣⎡=-1715171141C⎥⎥⎥⎥⎦⎤171141715⎥⎦⎤⎢⎣⎡-=6.44.13L A T即⎢⎢⎢⎢⎣⎡==-171517114ˆ1L A C X T ⎥⎥⎥⎥⎦⎤171141715⎥⎦⎤⎢⎣⎡≈⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡-015.096.00152.0963.06.44.13由残余误差方程可得()()()0196.00152.03963.029.10326.00152.02963.09.00042.00152.0963.039.2321=⨯-⨯-=-=⨯--==+⨯-=ννν00146.0312=∑=i iν因为是等精度测量,故得数据21,l l 的标准差相同,为0382.023312=-=∑=i iνσ不定系数()2,1,=j i d ij是矩阵1-C 中的各元素,即⎢⎣⎡=-21111d d C⎥⎦⎤2212d d ⎢⎢⎢⎢⎣⎡=171517114⎥⎥⎥⎥⎦⎤171141715则0819.0171140819.0171142211====d d可得估计量的标准差为0109.00819.00382.00109.00819.00382.0222111======d d x x σσσσ6—3解:由题意可写出误差方程()()⎪⎪⎩⎪⎪⎨⎧+--=+-=--=-=32431322113.05.04.03.0x x x x x x νννν设有列向量⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=43214321,ˆ,3.05.04.03.0ννννV x x x x X L 矩阵⎢⎢⎢⎢⎣⎡=0101A1010⎥⎥⎥⎥⎦⎤1100由式(6—21)及式(6—22)可得()L A AA XTT1ˆ-=⎢⎢⎢⎣⎡==01A A C T10 101⎥⎥⎥⎦⎤110⎢⎢⎢⎢⎣⎡0101101⎥⎥⎥⎥⎦⎤1100=⎢⎢⎢⎣⎡102 120⎥⎥⎥⎦⎤211⎢⎢⎢⎣⎡-=-213411C231-⎥⎥⎥⎦⎤--422⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=2.07.08.0L A T⎢⎢⎢⎣⎡-=21341ˆX231- ⎥⎥⎥⎦⎤--422⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-≈⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-15.042.032.015.0425.0325.06.07.13.1412.07.08.0可解得⎪⎩⎪⎨⎧=-==15.042.032.0321x x x (修正后结果)将最佳估计值代入误差方程,可得()()FF F F F F F F F F F F 03.015.042.03.003.015.032.05.002.042.04.002.032.03.04321-=+---==+-==+-=-=-=νννν2242322214120026.0Fi i=+++=∑=ννννν因为是等精度测量,故得数据4321,,,l l l l 的标准差相同,为FF i i051.00026.034412==-=∑=νσ不定系数()32,1,,=j i d ij是矩阵1-C 中的各元素,即 ⎢⎢⎢⎣⎡=-3121111d d d C232212d d d⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤21341332313d d d231-⎥⎥⎥⎦⎤--422则00.144175.034175.0341332211=⨯==⨯==⨯=d d d可得估计量的标准差为FF d FF d F F d x x x 051.000.1051.0044.075.0051.0044.075.0051.0333222111=========σσσσσσ6-5解:由题意可写出误差方程11223125.264.9410.14()x x x x ννν=-⎧⎪=-⎨⎪=-+⎩ 设有列向量112235.26ˆ4.94,10.14x L X V x ννν⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦矩阵100111A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦可得()L A AA XTT1ˆ-=10101 2 1*01011 12 11TC A A ⎡⎤⎡⎤⎡⎤⎢⎥==⨯=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎣⎦10.6667 -0.3333-0.3333 0.6667C-⎡⎤=⎢⎥⎣⎦15.400015.0800TA L ⎡⎤=⎢⎥⎣⎦1 5.2400ˆ 4.9200T X C A L -⎡⎤==⎢⎥⎣⎦可解得12 5.24000 4.9200x x =⎧⎨=⎩将最佳估计值代入误差方程,可得123 0.02m m -0.02m m -0.02m mννν===3222212310.0044ii m mνννν==++=∑因为是等精度测量,故得数据4321,,,l l l l 的标准差相同,为3210.03532ii m mνσ==≈-∑不定系数()32,1,,=j i d ij是矩阵1-C 中的各元素,即 1112121220.6667 -0.3333-0.3333 0.6667d d Cd d -⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎣⎦则11220.66670.6667d d ==可得估计量的标准差为1112220.0290.029x x d m m d m mσσσσ====6—6解:令2211,x C x C ==则由题意可写出误差方程()⎪⎩⎪⎨⎧+-=-=-=2132211121.3008.1105.2x x x x ννν 设有列向量⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=32121,ˆ,121.3008.1015.2νννV x x X L矩阵⎢⎢⎢⎣⎡=101A⎥⎥⎥⎦⎤110则由误差方程可列出正规方程为()L A XA AT T=ˆ即为⎩⎨⎧=+=+129.42226.522121x x x x按矩阵求解,则有⎢⎣⎡==12A A C T⎥⎦⎤21⎢⎢⎢⎢⎣⎡-=-31321C⎥⎥⎥⎥⎦⎤-3231⎥⎦⎤⎢⎣⎡=129.4226.5L A T⎢⎢⎢⎢⎣⎡-==-3132ˆ1L A C X T⎥⎥⎥⎥⎦⎤-3231⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡011.1108.2129.4226.5 可解得⎩⎨⎧==Fx Fx μμ011.1108.221将为佳估计值代入误差方程,可得()()FF F x x l FF F x l F F F x l μμμνμμμνμμμν002.0011.1108.2121.3003.0011.1008.1003.0108.2105.22133222111=+-=+-=-=-=-=-=-=-=32222212310.000022ii Fννννμ==++=∑因为是等精度测量,故测得数据321,,l l l 的标准差相同,为3210.0000220.0047321ii F Fνσμμ====-∑不定系数()2,1,=j i d ij是矩阵1-C中的各元素,即⎢⎣⎡=-21111d d C⎢⎣⎡-=⎥⎦⎤12312212d d⎥⎦⎤-21则667.0231667.02312211=⨯==⨯=d d可得估计量的标准差为1112220.00470.6670.00380.00470.6670.0038x x d F F d F Fσσμμσσμμ======6—8解: 正规方程为:4562.21431.5x y x y -=-+=解得 1.42.4x y ⎛⎫⎛⎫=⎪ ⎪⎝⎭⎝⎭;即 1.4x =, 2.4y = (有效位数不同会导致结果不同)求残差代入数据后得到3210.0420.0130.010.09i i i p ν==+⨯+⨯=∑于是标准差为3210.090.31i ii p n tνσ====-∑由式子5-51及所给正规方程的系数得i a i1 a i2 p i p i a i12 p i a i22 p i a i1a i2 l i p i a i1l i p i a i2l i 1 1 -3 1 1 9 -3 -5.6 -5.6 16.8 2 4 1 2 32 2 8 8.1 64.8 16.2 32 -13 12 3 -6 0.5 3-1.5∑4514-162.2 31.5451140450141x y x y x y x y -=-+=-=-+=11220.022;0.071d d ==6-9设L=x, C=y 则由题意可写出误差方程12310.8(5)510.2(2)210.3(0.5)0.5x yx y x y ννν⎧=--⎪⎪⎪=--⎨⎪⎪=---⎪⎩以0.17,3.5为待估计量x,y 的近似值,则估计量可表示为120.173.5x y δδ=+⎧⎨=+⎩'1'2'310.8(50.17)0.00715 3.510.2(20.17)0.0029 2 3.510.3(0.50.17)0.1860.5 3.5l l l ⎧=-⨯-=⎪⨯⎪⎪=-⨯-=⎨⨯⎪⎪=--⨯-=⎪⨯⎩111101202221022033310320()5,()0.016()2,()0.041()0.5,()0.16f f a a x yf f a a x yf f a a xy∂∂====∂∂∂∂====∂∂∂∂====∂∂化成线性方程为:11220.30.0220.040.30.0710.08k k d d σσσσ======1122123120.0071(50.016)0.0029(20.041)0.186(50.16)νδδνδδνδδ=-+⎧⎪=-+⎨⎪=-+⎩ 设有列向量112230.0071ˆ0.0029,,0.186L X V νδνδν⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦ 矩阵50.01620.0410.50.016A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦按矩阵求解,则有29.2500 0.24200.2420 0.0275TC A A ⎡⎤==⎢⎥⎣⎦10.0369 -0.3240 -0.3240 39.1622C-⎡⎤=⎢⎥⎣⎦-0.46 0.0279T A L ⎡⎤=⎢⎥⎣⎦ 1 -0.0261ˆ 1.2419T X C A L -⎡⎤==⎢⎥⎣⎦0.17(0.0261)0.14 3.5 1.2419 4.74x y =+-≈⎧⎨=+≈⎩代入误差方程为12310.8(5)0.142510.2(2) 0.0255 210.3(0.5)0.05190.5x yx y x y ννν⎧=--≈⎪⎪⎪=--≈⎨⎪⎪=---≈⎪⎩3222212310.0235ii νννν==++≈∑因为是等精度测量,故测得数据,L C 的标准差相同,为3210.02350.15321ii νσ===≈-∑不定系数()2,1,=j i d ij是矩阵1-C中的各元素,即1112121220.0369 -0.3240-0.3240 39.1622d d Cd d -⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎣⎦ 则11220.80110.03660.190.8011902.1351)24.06L C d d σσσσ==≈==≈当w=3时,1130.140.353 4.74x w L w C=-=⨯-≈⨯方差22222222222212222222211()()()()()11(0.14)0.0190.19()24.060.0002670.32490.12749 4.744.7430.45250.67x w l c w x x x L w wlcC wC wσσσσσσσ∂∂∂=++=+++∂∂∂=+⨯+⨯+=++⨯⨯=≈。
北航矩阵
2. (1) 设 A
1 (21) A 1 2 1 1 2 2 2 的满秩分解为 4
k k
.
(22) 备用:如果 AC , BD 有意义,则 ( A B)(C D) ( AC ) ( BD) (23) 备用:拉直公式 ABC A C T B (24) 备用: A , B 为方阵,则 AX XB C 有唯一解 A 和 B 没有公共 二.(18 分)计算下列各题
A ; ( A , 0) = 0
(9)设 A 的各列互相正交且模长为 1,则 A AH (10) A (aij ), 则 tr ( AH A) |aij | tr ( AAH ) |aij |
2 2 i, j i, j
(11) 若 tr ( A A) 0 则 A
H
(12) (正规阵性质)若 A 是上三角形正规阵,则 A 一定是
Bnn (13) 若 0
D 为正规阵, 则 D Cnn
a 0 2 1 (14) 备用: A , B , 则 A B 的特征根为 1 b 0 3
1
(15) A 0.5
A iI2 i .
令 f ( x) cos(tx)
0 cos(it ) cos(tA) cos(it )(G1 G2 ) cos(it ) I cos(it ) 0 et e t I (用Euler公式) 2
转置矩阵实验报告
一、实验目的1. 理解矩阵转置的概念和性质。
2. 掌握矩阵转置的计算方法,包括普通矩阵和稀疏矩阵的转置。
3. 通过编程实现矩阵转置算法,并分析算法的复杂度。
4. 理解矩阵转置在数值计算中的应用。
二、实验原理矩阵转置是指将矩阵的行和列互换位置得到的新矩阵。
对于任意矩阵 \( A \) ,其转置矩阵记为 \( A^T \) 。
如果 \( A \) 是一个 \( m \times n \) 的矩阵,那么 \( A^T \) 是一个 \( n \times m \) 的矩阵。
三、实验内容1. 普通矩阵转置- 使用二维数组存储矩阵,实现普通矩阵的转置。
- 输入一个 \( m \times n \) 的矩阵,输出其转置矩阵 \( A^T \) 。
2. 稀疏矩阵转置- 使用三元组表示稀疏矩阵,实现稀疏矩阵的转置。
- 输入一个稀疏矩阵,输出其转置矩阵 \( A^T \) 。
3. 算法分析- 分析普通矩阵转置和稀疏矩阵转置算法的时间复杂度。
- 比较两种算法在处理不同类型矩阵时的效率。
四、实验步骤1. 普通矩阵转置- 定义一个二维数组 \( A \) 存储矩阵元素。
- 输入矩阵 \( A \) 的行数 \( m \) 和列数 \( n \) 。
- 输入矩阵 \( A \) 的元素。
- 遍历数组 \( A \),将元素 \( A[i][j] \) 放入新数组 \( A^T[j][i] \) 。
- 输出转置矩阵 \( A^T \) 。
2. 稀疏矩阵转置- 定义一个结构体存储三元组,包括行号、列号和元素值。
- 输入稀疏矩阵的非零元素个数 \( t \) ,行数 \( m \) 和列数 \( n \) 。
- 输入稀疏矩阵的三元组表示。
- 遍历三元组表,将每个三元组 \( (i, j, e) \) 改为 \( (j, i, e) \) 。
- 输出转置矩阵 \( A^T \) 的三元组表示。
3. 算法分析- 普通矩阵转置的时间复杂度为 \( O(mn) \) ,空间复杂度为 \( O(mn) \) 。
多核数字信号处理器并行矩阵转置算法优化
多核数字信号处理器并行矩阵转置算法优化
裴向东;王庆林;廖林玉;李荣春;梅松竹;刘杰;庞征斌
【期刊名称】《国防科技大学学报》
【年(卷),期】2023(45)1
【摘要】矩阵转置是矩阵运算的基本操作,广泛应用于信号处理、科学计算以及深度学习等各种领域。
随着国防科技大学自主研制的飞腾异构多核数字信号处理器(digital signal processor, DSP)在各种领域中的推广应用,对高性能矩阵转置实现提出了强烈需求。
针对飞腾异构多核DSP的体系结构特征与矩阵转置操作的特点,提出了一种适配不同数据位宽(8 B、4 B以及2 B)矩阵的并行矩阵转置算法ftmMT。
该算法基于DSP中向量处理单元的Load/Store部件实现了向量化,同时基于矩阵分块实现了多个DSP核的并行处理,通过隐式乒乓设计实现了片上向量化转置与片外访存的重叠以及访存性能的大幅提升。
实验结果表明,ftmMT能够显著加快矩阵转置操作,与CPU上的开源转置库HPTT相比,可获得高达8.99倍的性能加速。
【总页数】10页(P57-66)
【作者】裴向东;王庆林;廖林玉;李荣春;梅松竹;刘杰;庞征斌
【作者单位】国防科技大学计算机学院;国防科技大学并行与分布处理国防科技重点实验室
【正文语种】中文
【中图分类】TP391
【相关文献】
1.多核处理器中并行自适应索引算法优化
2.基于C6678多核数字信号处理器的声纳信号并行处理设计
3.高性能多核数字信号处理器内核验证系统设计
4.高性能多核数字信号处理器
5.多核数字信号处理器矩阵乘卷积算法性能评测
因版权原因,仅展示原文概要,查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
fclose(fdA);
}
MPI_Bcast(&size,1,MPI_INT,0,MPI_COMM_WORLD);
t=(int)sqrt(p);
if (t>size)
t=size;
if(size%t!=0)
for(;;)
{
t--;
if(size%t==0)
break;
}
p=t*t;
m=size/t;
p=group_size;
if(my_rank==0)
{
starttime=MPI_Wtime();
fdA=fopen("dataIn.txt","r");
fscanf(fdA,"%d %d", &size, &N);
if(size != N)
{
puts("The input is error!");
{
v= my_rank/t;
u= my_rank%t;
for(i=0;i<m;i++)
for(j=0;j<m;j++)
b(i,j)=a(i,j);
MPI_Recv(a,m*m,MPI_FLOAT, (u*t+v), my_rank,MPI_COMM_WORLD,&status);
MPI_Send(b,m*m,MPI_FLOAT, (u*t+v),(u*t+v),MPI_COMM_WORLD);
double time1;
double time2;
int my_rank;
int p;
MPI_Status status;
FILE *fdA;
void Environment_Finalize(float *a,float *b)
{
free(a);
free(b);
}
int main(int argc, char **argv)
MPI_Finalize();
Environment_Finalize(a,b);
return(0);
}
exit(0);
}
A=(float*)malloc(floatsize*size*size);
B=(float*)malloc(floatsize*size*size);
for(i = 0; i < size; i ++)
{
for(j = 0; j < size; j ++) fscanf(fdA, "%f", A+i*size+j);
for(j=0;j<m;j++)
a(i,j)=A(i,j);
}
if (my_rank==0)
{
for(i=1;i<p;i++)
{
v=i/t;
u=i%t;
for(j=v*m;j<(v+1)*m;j++)
for(k=u*m;k<(u+1)*m;k++)
b((j%m),(k%m))=A(j,k);
MPI_Send(b,m*m,MPI_FLOAT,i,i,MPI_COMM_WORLD);
矩阵转置并行实现
1、算法描述:
若处理器个数为p,且它们的编号依次是0,1,…,p-1,则将n阶矩阵A分成p个大小为mxm的子块,m=[n/p]。p个子块组成一个子块阵列,记其中第i行第j列的子块为Aij,它含有第(i-1)m+1至第im行中的第(j-1)m+1至第jm列的所有元素。对每一处理器按行主方式赋以二维下标,记编号为i的处理器的二维下标为(v,u),其中v=[i/ ],u=imod ,将A的子块存入下标为(v,u)表示的对应处理器中。转置分为两步进行:第一步,子块转置;第二步,处理器内部转置。为了避免对应子块交换数据是处理器发生死锁,可令下三角块先向与之对应的上三角子块发送数据,然后从上三角子块接收数据;上三角子块数据先存放在缓冲区buffer中,然后从与之对应的下三角子块接收数据,最后再将缓冲区中的数据发送给下三角子块,流程图如下所示:
2、程序代码:
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#include "math.h"
#define E 0.0001
#define a(x,y) a[x*m+y]
#define b(x,y) b[x*m+y]
#define A(x,y) A[x*size+y]
}
}
else if(my_rank<p)
MPI_Send(a,m*m,MPI_FLOAT, 0,my_rank,MPI_COMM_WORLD);
if (my_rank==0)
{
printf("Input of file \"dataIn.txt\"\n");
printf("%d\t%d\n", size, size);
{
int i,j,k,my_rank,group_size;
float *a,*b;
int u,v;
float temp;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
printf("Distribute data time = %f seconds\n",time1-starttime);
printf("Parallel compute time = %f seconds\n",time2-time1);
}
MPI_Barrier(MPI_COMM_WORLD);
for(i=0;i<size;i++)
{
for(j=0;j<size;j++) printf("%f\t",A(i,j));
printf("\n");
}
printf("\nOutput of Matrix AT\n");
for(i=0;i<size;i++)
{
for(j=0;j<size;j++) printf("%f\t",B(i,j));
}
for(i=1;i<m;i++)
for(j=0;j<i;j++)
{
temp=a(i,j);
a(i,j)=a(j,i);
a(j,i)=temp;
}
if (my_rank==0)
{
for(i=0;i<m;i++)
for(j=0;j<m;j++)
B(i,j)=a(i,j);
}
if (my_rank==0)
u= my_rank%t;
MPI_Send(a,m*m,MPI_FLOAT, (u*t+v),(u*t+v),MPI_COMM_WORLD);
MPI_Recv(a,m*m,MPI_FLOAT, (u*t+v),my_rank,MPI_COMM_WORLD,&status);
}
if ((my_rank/t)<(my_rank%t)&&my_rank<p)
{
for(i=1;i<p;i++)
{
MPI_Recv(a,m*m,MPI_FLOAT,i,i,MPI_COMM_WORLD,&status);
v= i/t;
u= i%t;
for(j=v*m;j<(v+1)*m;j++)
for(k=u*m;k<(u+1)*m;k++)
B(j,k)=a((j%m),(k%m));
a=(float *)malloc(floatsize*m*m);
b=(float *)malloc(floatsize*m*m);
if (a==NULL||b==NULL)
printf("allocate space fail!");
if (my_rank==0)
{
for(i=0;i<m;i++)
printf("\n");
}
}
time2=MPI_Wtime();
if (my_rank==0)
{
printf("\n");
printf("processor number:%d\n",p);
printf("Whole running time = %f seconds\n",time2-starttime);
}
}
else if (my_rank<p)
MPI_Recv(a,m*m,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
time1=MPI_Wtime();
if ((my_rank/t)>(my_rank%t)&&my_rank<p)