数值分析计算实习题（二）SY1004114 全昌彪一：算法设计方案概述：本题采用fortran90语言编写程序，依据题目要求，采用带双步位移QR分解法求出所给矩阵的所有特征值，并求出相应于其实特征值的特征向量，以及相关需要给出的中间结果。

1、矩阵的A的初始化（赋值）：利用子函数initial（a,n）来实现，返回n×n 维二维数组a。

2、A矩阵的拟上三角化：利用子函数hessenberg(a,n)，在对矩阵进行QR分解前进行拟上三角化，这样可以提高计算效率，减少计算量，返回A矩阵的相似矩阵Hessenberg阵A(n-1)。

3、对A(n-1)进行带双步位移QR分解得出Ｃｍ及A矩阵的所有特征值，这一步利用了两个子函数eigenvalue(a,n,lamda,lamdai)和qrresolve(b,c,m)带双步位移QR分解可以加速收敛。

每次QR分解前先进行判断，若可以直接得到矩阵的特征值，则对矩阵直接降阶处理；若不可以，则进行QR分解，这样就进一步减少了计算量，提高了计算效率。

考虑到矩阵A可能有复特征值，采用两个一维数组lamda(n)及lamdai(n)分别存储其实部和虚部。

在双步位移处理及降阶过程中，被分解的矩阵Ａk(m ×m)及中间矩阵M k(m×m)的维数随m不断减少而降阶，于是引入了动态矩阵C(m×m)和B(m×m)分别存储，在使用前，先声明分配内存，使用结束后立即释放内存。

返回A(n-1)经双步位移QR分解后的矩阵及Ａ矩阵的所有特征值。

4、特征向量的求解：采用子函数eigenvector(a,lamda)实现求解A矩阵的属于实特征值的特征向量。

核心算法为高斯列主元消去法，(A-λI)x=b,b=0,回代过程令x(10)=1，即可求出对应于每一实特征值的特征向量的各个元素。

5、相关输出结果：所有数据均采用e型输出，数据保留到小数点后12位。

对于Ａ矩阵的拟上三角化结果集双步位移QR分解结果比较庞大，为了更好的显示，还采用了ｆ型输出，保留５位小数。

所有计算精度水平Ｅ＝10-12。

二：算法fortran源程序!此函数用于给A赋值，即初始化A矩阵subroutine initial(a,n)integer :: i,jdimension a(n,n)double precision ado i=1,ndo j=1,nif (i==j) thena(i,j)=1.5*cos(i+1.2*j)elsea(i,j)=sin(0.5*i+0.2*j)end ifend doend doend subroutine initial!此函数用于对Ａ矩阵进行拟上三角化subroutine hessenberg(a,n)implicit noneinteger::i,j,r,ndimension a(n,n),u(n),p(n),q(n),w(n) double precision a,u,p,q,w,d,t,c,h,z,e e=1.0d-12 !精度控制outer1:do r=1,n-2inner1:do i=1,np(i)=0.0q(i)=0.0w(i)=0.0end do inner1d=0.0t=0.0z=0inner2:do i=r+2,nif (a(i,r)/=0) thenz=1end ifend do inner2if (z==1) theninner3:do i=r+1,nd=sqrt(d**2+a(i,r)**2)end do inner3if (a(r+1,r)>0) thenc=-delsec=dend ifh=c**2-c*a(r+1,r)inner4:do i=1,nif (i<=r) thenu(i)=0else if (i==r+1) thenu(i)=a(r+1,r)-celseu(i)=a(i,r)end ifend do inner4inner5:do i=1,ndo j=1,np(i)=p(i)+u(j)*a(j,i)/hq(i)=q(i)+a(i,j)*u(j)/hend doend do inner5inner6:do i=1,nt=t+p(i)*u(i)end do inner6t=t/hinner7:do i=1,nw(i)=q(i)-t*u(i)end do inner7inner8:do i=1,ndo j=1,na(i,j)=a(i,j)-w(i)*u(j)-u(i)*p(j)end doend do inner8end ifend do outer1do i=1,ndo j=1,nif (abs(a(i,j))<=e) then!a(i,j)输出精度控制a(i,j)=0end ifend doend dowrite(*,*) ('矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)（位有效数字e型输出）:') write(*,100) ((a(i,j),j=1,n),i=1,n)write(*,*) ('矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)(保留五位小数f型输出):') write(*,200) ((a(i,j),j=1,n),i=1,n)!输出格式控制100format (10e20.12)200format (10f9.5)end subroutine hessenberg!QR分解子函数，求出Ａ（ｎ－１）矩阵QR分解的结果!主要是用于eigenvalue(a,n,lamda,lamdai)的调用subroutine qrresolve(b,c,m)implicit noneinteger::i,j,r,mdimension b(m,m),c(m,m),u(m),p(m),q(m),w(m),v(m)double precision b,c,u,p,q,w,v,d,c1,h,t,zouter1:do r=1,m-1inner1:do i=1,m!中间量初始化p(i)=0.0q(i)=0.0w(i)=0.0v(i)=0.0end do inner1d=0.0t=0.0z=0inner2:do i=r+1,mif (b(i,r)/=0) thenz=1end ifend do inner2if (z==1) theninner3:do i=r,md=sqrt(d**2+b(i,r)**2) end do inner3if (b(r,r)>0) thenc1=-delsec1=dend ifh=c1**2-c1*b(r,r)inner4:do i=1,mif (i<r) thenu(i)=0else if (i==r) thenu(i)=b(r,r)-c1elseu(i)=b(i,r)end ifend do inner4inner5:do i=1,mdo j=1,mv(i)=v(i)+u(j)*b(j,i)/h end doend do inner5inner6:do i=1,mdo j=1,mb(i,j)=b(i,j)-u(i)*v(j) end doend do inner6inner7:do i=1,mdo j=1,mp(i)=p(i)+u(j)*c(j,i)/hq(i)=q(i)+c(i,j)*u(j)/hend doend do inner7inner8:do i=1,mt=t+p(i)*u(i)end do inner8t=t/hinner9:do i=1,mw(i)=q(i)-t*u(i)end do inner9inner10:do i=1,mdo j=1,mc(i,j)=c(i,j)-w(i)*u(j)-u(i)*p(j)end doend do inner10end ifend do outer1end subroutine qrresolve!此函数用于求A矩阵的特征值（包含实特征值即复特征值）!核心算法为双步位移QR方法（中间调用了QR分解法子函数）!此程序最大的特点是只用了一个goto语句!（按课本上给出个goto语句将更易编写，笔者也编写了运行结果一样）!为了清晰显示，每一个循环语句，if语句基本上都有命名subroutine eigenvalue(a,n,lamda,lamdai)implicit nonedimension a(n,n),lamda(n),lamdai(n)double precision,allocatable,dimension(:,:) :: b,c !动态数组，用于保存Mk及C矩阵元素，维数是变化的integer :: i,j,k,m,n,Ldouble precision a,s,t,temp1,temp2,temp3,deta,lamda,lamdai,EE=1.0D-12 !精度m=nL=150 !循环迭代最大次数限定lamda=0lamdai=0do k=1,Louter1:if (abs(a(m,m-1))<=E) thenlamda(m)=a(m,m)lamdai(m)=0m=m-1 !m自减，降一维10 inner1:if (m==1) thenlamda(m)=a(1,1)lamdai(m)=0exitelse if (m==0) thenexitelsecycleend if inner1else!Dk子块求根temp1=(a(m-1,m-1)+a(m,m))/2.0temp2=a(m-1,m-1)*a(m,m)-a(m,m-1)*a(m-1,m)deta=temp1**2-temp2inner2:if (deta>0) thenlamda(m)=temp1+sqrt(deta)lamdai(m)=0lamda(m-1)=temp1-sqrt(deta)lamdai(m-1)=0else if (deta<0) thenlamda(m)=temp1lamda(m-1)=temp1lamdai(m)=sqrt(-deta)lamdai(m-1)=-sqrt(-deta)elselamda(m)=temp1lamda(m-1)=temp1lamdai(m)=0lamdai(m-1)=0end if inner2inner3: if (m==2) thenexitelseinner4: if (abs(a(m-1,m-2))<=E) thenm=m-2 !m自减，降两维goto 10elseinner5: if (k==L) thenexitelseallocate(b(m,m)) !动态数组分配内存声明allocate(c(m,m))inner6: do i=1,mdo j=1,mc(i,j)=a(i,j)end doend do inner6s=a(m-1,m-1)+a(m,m)t=a(m-1,m-1)*a(m,m)-a(m,m-1)*a(m-1,m) b=matmul(c,c) !fotran90可直接用矩阵乘法inner7: do i=1,mdo j=1,mb(i,j)=b(i,j)-s*c(i,j)if (i==j) thenb(i,j)=b(i,j)+tend ifend doend do inner7call qrresolve(b,c,m) !调用qr分解子函数inner8: do i=1,mdo j=1,ma(i,j)=c(i,j)end doend do inner8 !结束，释放可分配内存deallocate(b)deallocate(c)end if inner5end if inner4end if inner3end if outer1end dodo i=1,ndo j=1,nif (abs(a(i,j))<=e) then!a(i,j)输出精度控制a(i,j)=0end ifend doend dowrite(*,*) ('矩阵A(n-1)双步位移QR分解的结果（位有效数字e型输出）:') write(*,100) ((a(i,j),j=1,n),i=1,n)write(*,*) ('矩阵A(n-1)双步位移QR分解的结果（保留五位小数f型输出）:') write(*,200) ((a(i,j),j=1,n),i=1,n)write(*,*) ('矩阵A的全部特征值:')do i=1,nwrite(*,300) lamda(i),'+',lamdai(i),'i'end do!输出格式控制语句100format (10e20.12)200format (10f9.5)300format (1x,1e20.12,a2,1x,1e20.12,a1)end subroutine eigenvalue!此函数用于求矩阵A属于实特征值的特征向量!核心算法为高斯列主元消去法，(A-λI)x=b,b=0,x(10)=1回代subroutine eigenvector(a,lamda)parameter n=10dimension a(n,n),m(n),b(n),x(n),e(n,n),a_max(n),i_k(n) double precision a,lamda,b,x,t,a_max,sum,e,i_kinteger m,i,j,k,lcall initial(a,n) !调用子函数，对A进行初始化!单位矩阵do i=1,ne(i,i)=1end do!各参数初始化b=0m=0a_max=0t=0a=a-lamda*e !A-λI!消元过程do k=1,n-1do i=k,nif (abs(a(i,k))>=a_max(k)) thena_max(k)=a(i,k) !选列主元m(k)=iend ifend dodo j=k,nt=a(k,j)a(k,j)=a(m(k),j)a(m(k),j)=tt=b(k)b(k)=b(m(k))b(m(k))=tend dodo i=k+1,ni_k(i)=a(i,k)/a(k,k)do j=k+1,na(i,j)=a(i,j)-i_k(i)*a(k,j)end dob(i)=b(i)-i_k(i)*b(k)end doend do!回代过程x(n)=1do k=n-1,1,-1sum=0do j=k+1,nsum=sum+a(k,j)*x(j)end dox(k)=(b(k)-sum)/a(k,k)end dowrite (*,100) 'A的属于lamda=',lamda,'的特征向量:' write (*,200) x(:)100format (3x,a14,1e19.12,a11)200format (1x,1e20.12)end subroutine eigenvector!主函数，用于调用各个子函数program main_functionimplicit nonedimension a(10,10),lamda(10),lamdai(10),x(10)integer :: i,j,ndouble precision a,lamda,lamdai,xn=10call initial(a,n) !A矩阵初始化call hessenberg(a,n) !A(n-1)矩阵拟上三角化call eigenvalue(a,n,lamda,lamdai) !求A矩阵的所有特征值!求A矩阵属于实特征值的特征向量call eigenvector(a,lamda(1))call eigenvector(a,lamda(4))call eigenvector(a,lamda(5))call eigenvector(a,lamda(8))call eigenvector(a,lamda(9))call eigenvector(a,lamda(10))end program main_function三：程序输出结果汇总：矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)（12位有效数字e型输出）:-0.882751733065E+00 -0.993313316411E-01 -0.110334955595E+01 -0.760044676830E+00 0.154909623840E+00-0.194659168629E+01 -0.878226438751E-01 -0.925593370941E+00 0.603254375819E+00 0.151882370093E+00-0.234787833029E+01 0.237236973917E+01 0.181929100012E+01 0.323780721082E+00 0.220580321938E+000.210269272044E+01 0.181611739616E+00 0.127884456394E+01 -0.638050158135E+00 -0.415402766732E+000.000000000000E+00 0.172827457015E+01 -0.117146801290E+01 -0.124383914708E+01 -0.639976054127E+00-0.200283235265E+01 0.292496128495E+00 -0.641286862951E+00 0.978337580632E-01 0.255770647387E+000.000000000000E+00 0.000000000000E+00 -0.129166965220E+01 -0.111160388379E+01 0.117134648592E+01-0.130735637674E+01 0.180370952121E+00 -0.424641232157E+00 0.798856288534E-01 0.160878411848E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.156012616609E+01 0.812504527179E+000.442174836453E+00 -0.358869742717E-01 0.469176590186E+00 -0.273656195621E+00 -0.735907113693E-010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 -0.770777431331E+00-0.158305209469E+01 -0.304282234567E+00 0.252871461915E+00 -0.670993649410E+00 0.254459252711E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00-0.746342926455E+00 -0.270872578791E-01 -0.948652518297E+00 0.119586353771E+00 0.192915386004E-010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 -0.770179463681E+00 -0.469763787838E+00 0.498822338663E+00 0.113768224333E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.701315606266E+00 0.158220658011E+00 0.386262135075E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.484380788945E+00 0.399280281752E+00矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)(保留五位小数f型输出):-0.88275 -0.09933 -1.10335 -0.76004 0.15491 -1.94659 -0.08782 -0.92559 0.60325 0.15188-2.34788 2.37237 1.81929 0.32378 0.22058 2.10269 0.18161 1.27884 -0.63805 -0.415400.00000 1.72827 -1.17147 -1.24384 -0.63998 -2.00283 0.29250 -0.64129 0.09783 0.255770.00000 0.00000 -1.29167 -1.11160 1.17135 -1.30736 0.18037 -0.42464 0.07989 0.160880.00000 0.00000 0.00000 1.56013 0.81250 0.44217 -0.03589 0.46918 -0.27366 -0.073590.00000 0.00000 0.00000 0.00000 -0.77078 -1.58305 -0.30428 0.25287 -0.67099 0.254460.00000 0.00000 0.00000 0.00000 0.00000 -0.74634 -0.02709 -0.94865 0.11959 0.019290.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.77018 -0.46976 0.49882 0.113770.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.70132 0.15822 0.386260.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.48438 0.39928矩阵A(n-1)双步位移QR分解的结果（12位有效数字e型输出）:0.338303924721E+01 0.894875525475E+00 -0.895678244508E+00 0.846512632866E-01 0.261267793959E+000.161043313340E+01 -0.102255900702E+01 0.937223229328E-01 -0.100257738993E+01 -0.408625114186E+000.000000000000E+00 -0.211848218734E+01 -0.236152888518E+01 -0.345573222905E-01 -0.473663572352E-010.181640642861E+01 -0.231836979306E+00 -0.143548828808E+00 -0.653703238456E+00 0.322737505939E-010.000000000000E+00 0.355512159597E+00 -0.252851058727E+01 -0.637526428692E+00 0.202388052178E-010.183862769119E+01 0.186937662765E+00 -0.293254867876E+00 0.198706730246E+01 0.100463175221E+010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.157754818551E+01 -0.139621667385E-01-0.697199695782E+00 0.155545542881E+00 0.840521541031E-02 -0.815438441765E-01 -0.108612039515E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 -0.148403977622E+01-0.100530588729E+00 0.424954671617E-01 0.262346756601E-01 0.104004221001E+00 -0.118074774235E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00-0.694860436802E+00 -0.528901458582E+00 0.267916098808E+00 -0.596236908619E+00 -0.491158600319E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+000.178846387945E+00 -0.126620169756E+01 0.471455617468E-01 0.290478168004E+00 -0.357041410381E-010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.935588710088E+00 0.187740659272E+00 0.136024981245E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.636050766399E+00 0.273703606786E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.245568085231E-04 0.565162119254E-01矩阵A(n-1)双步位移QR分解的结果（保留五位小数f型输出）:3.38304 0.89488 -0.89568 0.08465 0.26127 1.61043 -1.02256 0.09372 -1.00258 -0.408630.00000 -2.11848 -2.36153 -0.03456 -0.04737 1.81641 -0.23184 -0.14355 -0.65370 0.032270.00000 0.35551 -2.52851 -0.63753 0.02024 1.83863 0.18694 -0.29325 1.987071.004630.00000 0.00000 0.00000 1.57755 -0.01396 -0.69720 0.15555 0.00841 -0.08154 -0.108610.00000 0.00000 0.00000 0.00000 -1.48404 -0.10053 0.04250 0.02623 0.10400 -0.118070.00000 0.00000 0.00000 0.00000 0.00000 -0.69486 -0.52890 0.26792 -0.59624 -0.491160.00000 0.00000 0.00000 0.00000 0.00000 0.17885 -1.26620 0.04715 0.29048 -0.035700.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.93559 0.18774 0.136020.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.63605 0.273700.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00002 0.05652矩阵A的全部特征值:0.338303924721E+01 + 0.000000000000E+00i-0.232349638731E+01 + -0.893040543164E+00i-0.232349638731E+01 + 0.893040543164E+00i0.157754818551E+01 + 0.000000000000E+00i-0.148403977622E+01 + 0.000000000000E+00i-0.980531067180E+00 + -0.113949139468E+00i-0.980531067180E+00 + 0.113949139468E+00i0.935588710088E+00 + 0.000000000000E+00i0.565046144245E-01 + 0.000000000000E+00i0.636062363900E+00 + 0.000000000000E+00iA的属于lamda= 0.338303924721E+01的特征向量: -0.436739834602E+00-0.906357312937E+00-0.196334018248E+01-0.108285729157E+01-0.126929785310E+01-0.107244858764E+010.362510824633E+000.168290589916E+010.211274357722E+010.100000000000E+01A的属于lamda= 0.157754818551E+01的特征向量: 0.621734685009E+00-0.111511734070E+00-0.248377525794E+01-0.130686124886E+01-0.381560716897E+010.811730907979E+01-0.123917001440E+01-0.680029414194E+000.269189633346E+010.100000000000E+01A的属于lamda=-0.148403977622E+01的特征向量: -0.173077810726E+020.240818399765E+020.413384989796E+00-0.857202156181E+010.928742697832E-01-0.783266447227E-01-0.637425984951E+00-0.340318975687E+00-0.378486114544E+000.100000000000E+01A的属于lamda= 0.935588710088E+00的特征向量: 0.279242914534E+010.159824395408E+01-0.520736954810E+00-0.166788993965E+01-0.122571220310E+020.724123969831E+01-0.539823439712E+010.284100862311E+02-0.121651397115E+020.100000000000E+01A的属于lamda= 0.565046144245E-01的特征向量:-0.510500561103E+01-0.488632284214E+010.950516391275E+01-0.678833013260E+00-0.960433257896E+01-0.304574395607E+010.157487434454E+02-0.739503192243E+01-0.710966608732E+010.100000000000E+01A的属于lamda= 0.636062363900E+00的特征向量:0.474503250282E+010.315786893686E+010.172995116830E+02-0.198004944728E+01-0.318752699418E+020.779400967078E+01-0.100425554520E+020.167075308672E+020.131053073601E+020.100000000000E+01请按任意键继续. . .四：编写程序过程中发现的一个问题本程序运行正确，但是在编写过程中曾经产生一个很小的错误，将矩阵M k=A k2-s A k+t I 表达式中错写成了M k=A k2-s A k-t I，依旧可以得到Ａ矩阵的所有特征值及属于实特征值的特征向量，而且和本程序的答案完全吻合，只是Ａ(n-1)进行带双步位移QR分解结束后得到的矩阵C m相同（因为和同学对比了计算结果，发现不对），于是我调试了很久（也重新写过带双步位移QR分解子程序）方才发现这个小变动，而这个小变动只影响带双步位移QR分解，并不影响特征值的求解。