127
128
129
130
131
132
133 REAL A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
134
135
136
137
138
139 REAL ZERO, HALF, ONE, TWO
140 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0,
141 $ two = 2.0e+0 )
142 REAL MULTPL
143 parameter( multpl = 4.0e+0 )
144
145
146 REAL AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
147 $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z, SAFMIN,
148 $ SAFMN2, SAFMX2
149 INTEGER COUNT
150
151
152 REAL SLAMCH, SLAPY2
154
155
156 INTRINSIC abs, max, min, sign, sqrt
157
158
159
162 safmn2 =
slamch(
'B' )**int( log( safmin / eps ) /
163 $ log(
slamch(
'B' ) ) / two )
164 safmx2 = one / safmn2
165 IF( c.EQ.zero ) THEN
166 cs = one
167 sn = zero
168
169 ELSE IF( b.EQ.zero ) THEN
170
171
172
173 cs = zero
174 sn = one
175 temp = d
176 d = a
177 a = temp
178 b = -c
179 c = zero
180
181 ELSE IF( (a-d).EQ.zero .AND. sign( one, b ).NE.
182 $ sign( one, c ) ) THEN
183 cs = one
184 sn = zero
185
186 ELSE
187
188 temp = a - d
189 p = half*temp
190 bcmax = max( abs( b ), abs( c ) )
191 bcmis = min( abs( b ), abs( c ) )*sign( one, b )*sign( one, c )
192 scale = max( abs( p ), bcmax )
193 z = ( p / scale )*p + ( bcmax / scale )*bcmis
194
195
196
197
198 IF( z.GE.multpl*eps ) THEN
199
200
201
202 z = p + sign( sqrt( scale )*sqrt( z ), p )
203 a = d + z
204 d = d - ( bcmax / z )*bcmis
205
206
207
209 cs = z / tau
210 sn = c / tau
211 b = b - c
212 c = zero
213
214 ELSE
215
216
217
218
219 count = 0
220 sigma = b + c
221 10 CONTINUE
222 count = count + 1
223 scale = max( abs(temp), abs(sigma) )
224 IF( scale.GE.safmx2 ) THEN
225 sigma = sigma * safmn2
226 temp = temp * safmn2
227 IF (count .LE. 20)
228 $ GOTO 10
229 END IF
230 IF( scale.LE.safmn2 ) THEN
231 sigma = sigma * safmx2
232 temp = temp * safmx2
233 IF (count .LE. 20)
234 $ GOTO 10
235 END IF
236 p = half*temp
237 tau =
slapy2( sigma, temp )
238 cs = sqrt( half*( one+abs( sigma ) / tau ) )
239 sn = -( p / ( tau*cs ) )*sign( one, sigma )
240
241
242
243
244 aa = a*cs + b*sn
245 bb = -a*sn + b*cs
246 cc = c*cs + d*sn
247 dd = -c*sn + d*cs
248
249
250
251
252 a = aa*cs + cc*sn
253 b = bb*cs + dd*sn
254 c = -aa*sn + cc*cs
255 d = -bb*sn + dd*cs
256
257 temp = half*( a+d )
258 a = temp
259 d = temp
260
261 IF( c.NE.zero ) THEN
262 IF( b.NE.zero ) THEN
263 IF( sign( one, b ).EQ.sign( one, c ) ) THEN
264
265
266
267 sab = sqrt( abs( b ) )
268 sac = sqrt( abs( c ) )
269 p = sign( sab*sac, c )
270 tau = one / sqrt( abs( b+c ) )
271 a = temp + p
272 d = temp - p
273 b = b - c
274 c = zero
275 cs1 = sab*tau
276 sn1 = sac*tau
277 temp = cs*cs1 - sn*sn1
278 sn = cs*sn1 + sn*cs1
279 cs = temp
280 END IF
281 ELSE
282 b = -c
283 c = zero
284 temp = cs
285 cs = -sn
286 sn = temp
287 END IF
288 END IF
289 END IF
290
291 END IF
292
293
294
295 rt1r = a
296 rt2r = d
297 IF( c.EQ.zero ) THEN
298 rt1i = zero
299 rt2i = zero
300 ELSE
301 rt1i = sqrt( abs( b ) )*sqrt( abs( c ) )
302 rt2i = -rt1i
303 END IF
304 RETURN
305
306
307
real function slamch(cmach)
SLAMCH
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).