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
|
subroutine wexchn(ar, ai, vr, vi, n, l, fail, na, nv)
c!purpose
c
c given the upper triangular complex matrix a ,wexchn produce a
c unitary transformation which exchange the two consecutive blocks
c starting at a(l,l),along with their eigenvalues.
c the transformation is accumulated in v.
c!calling sequence
c
c subroutine exchng(ar, ai, vr, vi, n, l, fail, na, nv)
c
c integer l, na, nv
c double precision ar, ai, vr, vi
c dimension ar(na,n) , ai(na,n) ,vr(nv,n) ,vi(nv,n)
c logical fail
c
c starred parameters are altered by the subroutine
c
c *ar,ai the matrix whose blocks are to be
c interchanged.
c *vr,vi the array into which the transformations
c are to re accumulated.
c n the order of the matrix a.
c l the position of the blocks.
c *fail a logical variable which is false on a
c normal return. if thirty iterations were
c performed without convergence, fail is set
c to true and the element
c a(l+b2,l+b2-1) cannot be assumed zero.
c na the first dimension of the array a.
c nv the first dimension of the array v.
c
c!auxiliary routines
c max sqrt abs (fortran)
c!originator
c steer i.n.r.i.a from routine exchng
c!
c
c
c Copyright INRIA
integer l,na,nv
double precision ar,ai,vr,vi
dimension ar(na, n), ai(na, n), vr(nv, n), vi(nv, n)
logical fail
c
c internal variables.
c
double precision pr,pi,qr,qi,r,sr,si,tr,ti,zero
integer i, j, l1
data zero /0.0d+0/
l1 = l + 1
c
fail = .false.
c
c interchange 1x1 and 1x1 blocks.
c
qr = ar(l1,l1) - ar(l,l)
pr = ar(l,l1)
qi = ai(l1,l1) - ai(l,l)
pi = ai(l,l1)
r = max(abs(pr),abs(pi),abs(qr),abs(qi))
if (r.eq.zero) return
pr = pr/r
qr = qr/r
pi = pi/r
qi = qi/r
r = sqrt(pr*pr + pi*pi + qr*qr + qi*qi)
pr = pr/r
qr = qr/r
pi = pi/r
qi = qi/r
do 10 j = l,n
sr = ar(l,j)
si = ai(l,j)
tr = ar(l1,j)
ti = ai(l1,j)
ar(l,j) = pr*sr + pi*si + qr*tr + qi*ti
ai(l,j) = pr*si - pi*sr + qr*ti - qi*tr
ar(l1,j) = pr*tr - pi*ti - qr*sr + qi*si
ai(l1,j) = pr*ti + pi*tr - qr*si - qi*sr
10 continue
do 20 i = 1,l1
sr = ar(i,l)
si = ai(i,l)
tr = ar(i,l1)
ti = ai(i,l1)
ar(i,l) = pr*sr + qr*tr - pi*si - qi*ti
ai(i,l) = pi*sr + qi*tr + pr*si + qr*ti
ar(i,l1) = pr*tr + pi*ti - qr*sr - qi*si
ai(i,l1) = pr*ti - pi*tr - qr*si + qi*sr
20 continue
do 30 i = 1,n
sr = vr(i,l)
si = vi(i,l)
tr = vr(i,l1)
ti = vi(i,l1)
vr(i,l) = pr*sr + qr*tr - pi*si - qi*ti
vi(i,l) = pi*sr + qi*tr + pr*si + qr*ti
vr(i,l1) = pr*tr + pi*ti - qr*sr - qi*si
vi(i,l1) = pr*ti - pi*tr - qr*si + qi*sr
30 continue
ar(l1,l) = zero
ai(l1,l) = zero
return
end
|
