1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
subroutine split(a, v, n, l, e1, e2, na, nv)
c
c!purpose
c
c given the upper hessenberg matrix a with a 2x2 block
c starting at a(l,l), split determines if the
c corresponding eigenvalues are real or complex, if they
c are real, a rotation is determined that reduces the
c block to upper triangular form with the eigenvalue
c of largest absolute value appearing first. the
c rotation is accumulated in v. the eigenvalues (real
c or complex) are returned in e1 and e2.
c!calling sequence
c
c subroutine split(a, v, n, l, e1, e2, na, nv)
c
c double precision a,v,e1,e2
c integer n,l,na,nv
c dimension a(na,n),v(nv,n)
c
c starred parameters are altered by the subroutine
c
c *a the upper hessenberg matrix whose 2x2
c block is to be dsplit.
c *v the array in which the dsplitting trans-
c formation is to be accumulated.
c n the order of the matrix a.
c l the position of the 2x2 block.
c *e1 on return if the eigenvalues are complex
c *e2 e1 contains their common real part and
c e2 contains the positive imaginary part.
c if the eigenvalues are real. e1 contains
c the one largest in absolute value and f2
c contains the other one.
c na the first dimension of the array a.
c nv the first dimension of the array v.
c!auxiliary routines
c abs sqrt (fortran)
c!
c originator
c
double precision a,v,e1,e2
integer n,l,na,nv
dimension a(na,n),v(nv,n)
c internal variables
c
c internal variables
double precision p,q,r,t,u,w,x,y,z,zero,two
integer i,j,l1
data zero, two /0.0d+0,2.0d+0/
l1 = l + 1
c
x = a(l1,l1)
y = a(l,l)
w = a(l,l1)*a(l1,l)
p = (y-x)/two
q = p**2 + w
if (q.ge.zero) go to 10
c
c complex eigenvalue.
c
e1 = p + x
e2 = sqrt(-q)
return
10 continue
c
c two real eigenvalues. set up transformation.
c
z = sqrt(q)
if (p.lt.zero) go to 20
z = p + z
go to 30
20 continue
z = p - z
30 continue
if (z.eq.zero) go to 40
r = -w/z
go to 50
40 continue
r = zero
50 continue
if (abs(x+z).ge.abs(x+r)) z = r
y = y - x - z
x = -z
t = a(l,l1)
u = a(l1,l)
if (abs(y)+abs(u).le.abs(t)+abs(x)) go to 60
q = u
p = y
go to 70
60 continue
q = x
p = t
70 continue
r = sqrt(p**2+q**2)
if (r.gt.zero) go to 80
e1 = a(l,l)
e2 = a(l1,l1)
a(l1,l) = zero
return
80 continue
p = p/r
q = q/r
c
c premultiply.
c
do 90 j=l,n
z = a(l,j)
a(l,j) = p*z + q*a(l1,j)
a(l1,j) = p*a(l1,j) - q*z
90 continue
c
c postmultiply.
c
do 100 i=1,l1
z = a(i,l)
a(i,l) = p*z + q*a(i,l1)
a(i,l1) = p*a(i,l1) - q*z
100 continue
c
c accumulate the transformation in v.
c
do 110 i=1,n
z = v(i,l)
v(i,l) = p*z + q*v(i,l1)
v(i,l1) = p*v(i,l1) - q*z
110 continue
a(l1,l) = zero
e1 = a(l,l)
e2 = a(l1,l1)
return
end
|
