| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine calcz(d,e,n,dmax,dmin,z,ierr)
|
|---|
| 3 | c
|
|---|
| 4 | c Num. Rec. 2nd Ed., p.473
|
|---|
| 5 | c
|
|---|
| 6 | c Note: d(1:n) is the diagonal of the sym. tridiagonal matrix
|
|---|
| 7 | c e(1:n) is the upper diagonal of the tridiagonal matrix,
|
|---|
| 8 | c WITH e(n) ARBITRARY (a slight shift from Num.Rec.)
|
|---|
| 9 | c
|
|---|
| 10 | c z(n:n) is the packed array of eigenvectors
|
|---|
| 11 | c
|
|---|
| 12 | real d(n),e(n),z(n,n)
|
|---|
| 13 | real smalln,small
|
|---|
| 14 | c
|
|---|
| 15 | call ident(z,n)
|
|---|
| 16 | one = 1.
|
|---|
| 17 | c
|
|---|
| 18 | c Find smallest number (pff mod to Num. Rec. 2nd Ed., p.473)
|
|---|
| 19 | c
|
|---|
| 20 | small = 0.5
|
|---|
| 21 | do i = 1,100
|
|---|
| 22 | smalln = small * small
|
|---|
| 23 | if (smalln .eq.0) then
|
|---|
| 24 | do j=1,1000
|
|---|
| 25 | smalln = 0.5*small
|
|---|
| 26 | if (smalln .eq.0) goto 10
|
|---|
| 27 | small = smalln
|
|---|
| 28 | enddo
|
|---|
| 29 | endif
|
|---|
| 30 | small = smalln
|
|---|
| 31 | enddo
|
|---|
| 32 | 10 continue
|
|---|
| 33 | small = 10.*small
|
|---|
| 34 | small = max(small,1e-99)
|
|---|
| 35 |
|
|---|
| 36 | c write(6,*) 'this is small:',small
|
|---|
| 37 | c
|
|---|
| 38 | do 15 l=1,n
|
|---|
| 39 | iter = 0
|
|---|
| 40 | c
|
|---|
| 41 | 1 do 12 m=l,n-1
|
|---|
| 42 | dd = abs( d(m) ) + abs( d(m+1) )
|
|---|
| 43 | de = e(m) + dd
|
|---|
| 44 | df = abs(dd - de)
|
|---|
| 45 | c write(6,112) iter,m,'de:',dd,de,df,small
|
|---|
| 46 | if ( df .le. small ) goto 2
|
|---|
| 47 | 12 continue
|
|---|
| 48 | 112 format(i3,i9,1x,a3,1p4e16.8)
|
|---|
| 49 | m = n
|
|---|
| 50 | c
|
|---|
| 51 | 2 if ( m .ne. l ) then
|
|---|
| 52 | if ( iter .eq. 600 ) then
|
|---|
| 53 | write (6,*) 'too many iterations in calc'
|
|---|
| 54 | c n10 = min(n,10)
|
|---|
| 55 | c do i=1,n
|
|---|
| 56 | c write(6,9) d(i),(z(i,j),j=1,n10)
|
|---|
| 57 | c enddo
|
|---|
| 58 | c 9 format(1pe12.4,' e:',1p10e12.4)
|
|---|
| 59 | c call exitt
|
|---|
| 60 | ierr=1
|
|---|
| 61 | return
|
|---|
| 62 | endif
|
|---|
| 63 | c
|
|---|
| 64 | iter = iter + 1
|
|---|
| 65 | g = ( d(l+1) - d(l) ) / ( 2.0 * e(l) )
|
|---|
| 66 | r = pythag(g,one)
|
|---|
| 67 | g = d(m) - d(l) + e(l)/(g+sign(r,g))
|
|---|
| 68 | s = 1.0
|
|---|
| 69 | c = 1.0
|
|---|
| 70 | p = 0.0
|
|---|
| 71 | c
|
|---|
| 72 | do 14 i = m-1,l,-1
|
|---|
| 73 | f = s * e(i)
|
|---|
| 74 | b = c * e(i)
|
|---|
| 75 | r = pythag(f,g)
|
|---|
| 76 | e(i+1) = r
|
|---|
| 77 | if ( abs(r) .le. small ) then
|
|---|
| 78 | d(i+1) = d(i+1) - p
|
|---|
| 79 | e(m) = 0.
|
|---|
| 80 | goto 1
|
|---|
| 81 | endif
|
|---|
| 82 | s = f/r
|
|---|
| 83 | c = g/r
|
|---|
| 84 | g = d(i+1) - p
|
|---|
| 85 | r = ( d(i)-g )*s + 2.*c*b
|
|---|
| 86 | p = s*r
|
|---|
| 87 | d(i+1) = g + p
|
|---|
| 88 | g = c*r - b
|
|---|
| 89 | c ... find eigenvectors ... (new, 11/19/94, pff, p.363. Num.Rec.I.)
|
|---|
| 90 | do 13 k=1,n
|
|---|
| 91 | f = z(k,i+1)
|
|---|
| 92 | z(k,i+1) = s*z(k,i)+c*f
|
|---|
| 93 | z(k,i ) = c*z(k,i)-s*f
|
|---|
| 94 | 13 continue
|
|---|
| 95 | c ... end of eigenvector section ...
|
|---|
| 96 | 14 continue
|
|---|
| 97 | c
|
|---|
| 98 | d(l) = d(l) - p
|
|---|
| 99 | e(l) = g
|
|---|
| 100 | e(m) = 0.0
|
|---|
| 101 | goto 1
|
|---|
| 102 | endif
|
|---|
| 103 | c
|
|---|
| 104 | 15 continue
|
|---|
| 105 | c
|
|---|
| 106 | c write (8,8) (d(j),j=1,n)
|
|---|
| 107 | c 8 format('eig:',8f10.4)
|
|---|
| 108 | c
|
|---|
| 109 | dmin = d(1)
|
|---|
| 110 | dmax = d(1)
|
|---|
| 111 | do 40 i = 1 , n
|
|---|
| 112 | dmin = min( d(i) , dmin )
|
|---|
| 113 | dmax = max( d(i) , dmax )
|
|---|
| 114 | 40 continue
|
|---|
| 115 | c
|
|---|
| 116 | c Output eigenvectors
|
|---|
| 117 | c
|
|---|
| 118 | c n10 = min(n,10)
|
|---|
| 119 | c do i=1,n
|
|---|
| 120 | c write(6,9) d(i),(z(i,j),j=1,n10)
|
|---|
| 121 | c enddo
|
|---|
| 122 | c 9 format(1pe12.4,' e:',1p10e12.4)
|
|---|
| 123 | c
|
|---|
| 124 | ierr=0
|
|---|
| 125 | return
|
|---|
| 126 | end
|
|---|
| 127 | c-----------------------------------------------------------------------
|
|---|
| 128 | function pythag(a,b)
|
|---|
| 129 | c
|
|---|
| 130 | c Compute sqrt(a*a + b*b) w/o under- or over-flow.
|
|---|
| 131 | c
|
|---|
| 132 | absa=abs(a)
|
|---|
| 133 | absb=abs(b)
|
|---|
| 134 | if (absa.gt.absb) then
|
|---|
| 135 | pythag = absa*sqrt(1. + (absb/absa)**2 )
|
|---|
| 136 | else
|
|---|
| 137 | if (absb.eq.0.) then
|
|---|
| 138 | pythag = 0.
|
|---|
| 139 | else
|
|---|
| 140 | pythag = absb*sqrt(1. + (absa/absb)**2 )
|
|---|
| 141 | endif
|
|---|
| 142 | endif
|
|---|
| 143 | return
|
|---|
| 144 | end
|
|---|
| 145 | c-----------------------------------------------------------------------
|
|---|
| 146 | subroutine ident(a,n)
|
|---|
| 147 | real a(n,n)
|
|---|
| 148 | call rzero(a,n*n)
|
|---|
| 149 | do i=1,n
|
|---|
| 150 | a(i,i) = 1.0
|
|---|
| 151 | enddo
|
|---|
| 152 | return
|
|---|
| 153 | end
|
|---|
| 154 | c-----------------------------------------------------------------------
|
|---|
