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 REAL CR, CI
59
60
61 EXTERNAL clartg, sladiv
62
63
64 INTRINSIC real,
cmplx, conjg, aimag, sqrt
65
66
67
68
69
70 cs = rone
71 sn = zero
72
73 IF( c.EQ.zero ) THEN
74 GO TO 10
75
76 ELSE IF( b.EQ.zero ) THEN
77
78
79
80 cs = rzero
81 sn = one
82 temp = d
83 d = a
84 a = temp
85 b = -c
86 c = zero
87 GO TO 10
88 ELSE IF( ( a-d ).EQ.zero ) THEN
89 temp = sqrt( b*c )
90 a = a + temp
91 d = d - temp
92 IF( ( b+c ).EQ.zero ) THEN
93 cs = sqrt( half )
94 sn =
cmplx( rzero, rone )*cs
95 ELSE
96 temp = sqrt( b+c )
97 CALL sladiv( real( sqrt( b ) ), aimag( sqrt( b ) ),
98 $ real( temp ), aimag( temp ), cr, ci )
99 temp2 =
cmplx( cr, ci )
100 cs = real( temp2 )
101 CALL sladiv( real( sqrt( c ) ), aimag( sqrt( c ) ),
102 $ real( temp ), aimag( temp ), cr, ci )
104 END IF
105 b = b - c
106 c = zero
107 GO TO 10
108 ELSE
109
110
111
112 t = d
113 u = b*c
114 x = half*( a-t )
115 y = sqrt( x*x+u )
116 IF( real( x )*real( y )+aimag( x )*aimag( y ).LT.rzero )
117 $ y = -y
118 CALL sladiv( real( u ), aimag( u ),
119 $ real( x+y ), aimag( x+y ), cr, ci )
120 t = t -
cmplx( cr, ci )
121
122
123
124
125 CALL clartg( a-t, c, cs, sn, aa )
126
127 d = d - t
128 bb = cs*b + sn*d
129 dd = -conjg( sn )*b + cs*d
130
131 a = aa*cs + bb*conjg( sn ) + t
132 b = -aa*sn + bb*cs
133 c = zero
134 d = t
135
136 END IF
137
138 10 CONTINUE
139
140
141
142 rt1 = a
143 rt2 = d
144 RETURN
145
146
147