矩陣對角化

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)) :: gg
10         LOGICAL(LGT) :: yesvec
11         n=assert_eq(size(a,1),size(a,2),size(d),size(e),'tred2')
12         if (present(novectors)) then
13                 yesvec=.not. novectors
14         else
15                 yesvec=.true.
16         end if
17         do i=n,2,-1
18                 l=i-1
19                 h=0.0
20                 if (l > 1) then
21                         scale=sum(abs(a(i,1:l)))
22                         if (scale == 0.0) then
23                                 e(i)=a(i,l)
24                         else
25                                 a(i,1:l)=a(i,1:l)/scale
26                                 h=sum(a(i,1:l)**2)
27                                 f=a(i,l)
28                                 g=-sign(sqrt(h),f)
29                                 e(i)=scale*g
30                                 h=h-f*g
31                                 a(i,l)=f-g
32                                 if (yesvec) a(1:l,i)=a(i,1:l)/h
33                                 do j=1,l
34                                         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)))/h
36                                 end do
37                                 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,l
41                                         a(j,1:j)=a(j,1:j)-a(i,j)*e(1:j)-e(j)*a(i,1:j)
42                                 end do
43                         end if
44                 else
45                         e(i)=a(i,l)
46                 end if
47                 d(i)=h
48         end do
49         if (yesvec) d(1)=0.0
50         e(1)=0.0
51         do i=1,n
52                 if (yesvec) then
53                         l=i-1
54                         if (d(i) /= 0.0) then
55                                 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 if
58                         d(i)=a(i,i)
59                         a(i,i)=1.0
60                         a(i,1:l)=0.0
61                         a(1:l,i)=0.0
62                 else
63                         d(i)=a(i,i)
64                 end if
65         end do
66         END SUBROUTINE tred2

問題出在22行,scale==0.0的判斷上。稍微有點常識的人都知道，計算機裏作浮點數判斷不能用==，必須用abs(scale-0.0)<episilon episilon爲一個很小的數。沒想到numerical recipe第三版的這本書裏還有這樣子的代碼。

鑑於scale>=0.0，因此改爲scale<=episilon就好了，結果就很穩定了。

吐槽一下，怪不得我總以爲結果怎麼那麼隨機，原來是這裏出的問題，這但是真隨機啊！