首页 > 代码库 > 矩阵对角化

矩阵对角化

2024-07-09 05:23:25   221人阅读

numerical recipe 里一共讲了两种实数对称矩阵的对角化，

jacobi法

tred2生成上三角阵以后用tqli对角化

前者稳定但慢易并行，后者较快但疑似不稳定，串行。

花了一下午，一点点调试终于知道了第二种方法不稳定的原因在哪里

 1         SUBROUTINE tred2(a,d,e,novectors) 2         USE nrtype; USE nrutil, ONLY : assert_eq,outerprod 3         IMPLICIT NONE 4         REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a 5         REAL(SP), DIMENSION(:), INTENT(OUT) :: d,e 6         LOGICAL(LGT), OPTIONAL, INTENT(IN) :: novectors 7         INTEGER(I4B) :: i,j,l,n 8         REAL(SP) :: f,g,h,hh,scale 9         REAL(SP), DIMENSION(size(a,1)) :: gg10         LOGICAL(LGT) :: yesvec11         n=assert_eq(size(a,1),size(a,2),size(d),size(e),‘tred2‘)12         if (present(novectors)) then13                 yesvec=.not. novectors14         else15                 yesvec=.true.16         end if17         do i=n,2,-118                 l=i-119                 h=0.020                 if (l > 1) then21                         scale=sum(abs(a(i,1:l)))22                         if (scale == 0.0) then23                                 e(i)=a(i,l)24                         else25                                 a(i,1:l)=a(i,1:l)/scale26                                 h=sum(a(i,1:l)**2)27                                 f=a(i,l)28                                 g=-sign(sqrt(h),f)29                                 e(i)=scale*g30                                 h=h-f*g31                                 a(i,l)=f-g32                                 if (yesvec) a(1:l,i)=a(i,1:l)/h33                                 do j=1,l34                                         e(j)=(dot_product(a(j,1:j),a(i,1:j)) &35                                         +dot_product(a(j+1:l,j),a(i,j+1:l)))/h36                                 end do37                                 f=dot_product(e(1:l),a(i,1:l))38                                 hh=f/(h+h)39                                 e(1:l)=e(1:l)-hh*a(i,1:l)40                                 do j=1,l41                                         a(j,1:j)=a(j,1:j)-a(i,j)*e(1:j)-e(j)*a(i,1:j)42                                 end do43                         end if44                 else45                         e(i)=a(i,l)46                 end if47                 d(i)=h48         end do49         if (yesvec) d(1)=0.050         e(1)=0.051         do i=1,n52                 if (yesvec) then53                         l=i-154                         if (d(i) /= 0.0) then55                                 gg(1:l)=matmul(a(i,1:l),a(1:l,1:l))56                                 a(1:l,1:l)=a(1:l,1:l)-outerprod(a(1:l,i),gg(1:l))57                         end if58                         d(i)=a(i,i)59                         a(i,i)=1.060                         a(i,1:l)=0.061                         a(1:l,i)=0.062                 else63                         d(i)=a(i,i)64                 end if65         end do66         END SUBROUTINE tred2

问题出在22行,scale==0.0的判断上。稍微有点常识的人都知道，计算机里做浮点数判断不能用==，必须用abs(scale-0.0)<episilon episilon为一个很小的数。没想到numerical recipe第三版的这本书里还有这样子的代码。

鉴于scale>=0.0，所以改成scale<=episilon就行了，结果就很稳定了。

吐槽一下，怪不得我总觉得结果怎么那么随机，原来是这里出的问题，这可是真随机啊！


联系
我们