2
3
4
5
6
7
8
9 REAL CS
10 COMPLEX A, B, C, D, RT1, RT2, SN
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 REAL RZERO, HALF, RONE
50 parameter( rzero = 0.0e+0, half = 0.5e+0,
51 $ rone = 1.0e+0 )
52 COMPLEX ZERO, ONE
53 parameter( zero = ( 0.0e+0, 0.0e+0 ),
54 $ one = ( 1.0e+0, 0.0e+0 ) )
55
56
57 COMPLEX AA, BB, DD, T, TEMP, TEMP2, U, X, Y
58
59
60 COMPLEX CLADIV
62
63
65
66
67 INTRINSIC real,
cmplx, conjg, aimag, sqrt
68
69
70
71
72
73 cs = rone
74 sn = zero
75
76 IF( c.EQ.zero ) THEN
77 GO TO 10
78
79 ELSE IF( b.EQ.zero ) THEN
80
81
82
83 cs = rzero
84 sn = one
85 temp = d
86 d = a
87 a = temp
88 b = -c
89 c = zero
90 GO TO 10
91 ELSE IF( ( a-d ).EQ.zero ) THEN
92 temp = sqrt( b*c )
93 a = a + temp
94 d = d - temp
95 IF( ( b+c ).EQ.zero ) THEN
96 cs = sqrt( half )
97 sn =
cmplx( rzero, rone )*cs
98 ELSE
99 temp = sqrt( b+c )
100 temp2 =
cladiv( sqrt( b ), temp )
101 cs = real( temp2 )
102 sn =
cladiv( sqrt( c ), temp )
103 END IF
104 b = b - c
105 c = zero
106 GO TO 10
107 ELSE
108
109
110
111 t = d
112 u = b*c
113 x = half*( a-t )
114 y = sqrt( x*x+u )
115 IF( real( x )*real( y )+aimag( x )*aimag( y ).LT.rzero )
116 $ y = -y
117 t = t -
cladiv( u, ( x+y ) )
118
119
120
121
122 CALL clartg( a-t, c, cs, sn, aa )
123
124 d = d - t
125 bb = cs*b + sn*d
126 dd = -conjg( sn )*b + cs*d
127
128 a = aa*cs + bb*conjg( sn ) + t
129 b = -aa*sn + bb*cs
130 c = zero
131 d = t
132
133 END IF
134
135 10 CONTINUE
136
137
138
139 rt1 = a
140 rt2 = d
141 RETURN
142
143
144
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.