157
158
159
160
161
162
163 INTEGER LDA, LDB
164 REAL CSL, CSR, SNL, SNR
165
166
167 REAL A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ),
168 $ B( LDB, * ), BETA( 2 )
169
170
171
172
173
174 REAL ZERO, ONE
175 parameter( zero = 0.0e+0, one = 1.0e+0 )
176
177
178 REAL ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ,
179 $ R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1,
180 $ WR2
181
182
184
185
186 REAL SLAMCH, SLAPY2
188
189
190 INTRINSIC abs, max
191
192
193
196
197
198
199 anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
200 $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
201 ascale = one / anorm
202 a( 1, 1 ) = ascale*a( 1, 1 )
203 a( 1, 2 ) = ascale*a( 1, 2 )
204 a( 2, 1 ) = ascale*a( 2, 1 )
205 a( 2, 2 ) = ascale*a( 2, 2 )
206
207
208
209 bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
210 $ safmin )
211 bscale = one / bnorm
212 b( 1, 1 ) = bscale*b( 1, 1 )
213 b( 1, 2 ) = bscale*b( 1, 2 )
214 b( 2, 2 ) = bscale*b( 2, 2 )
215
216
217
218 IF( abs( a( 2, 1 ) ).LE.ulp ) THEN
219 csl = one
220 snl = zero
221 csr = one
222 snr = zero
223 a( 2, 1 ) = zero
224 b( 2, 1 ) = zero
225 wi = zero
226
227
228
229 ELSE IF( abs( b( 1, 1 ) ).LE.ulp ) THEN
230 CALL slartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
231 csr = one
232 snr = zero
233 CALL srot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
234 CALL srot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
235 a( 2, 1 ) = zero
236 b( 1, 1 ) = zero
237 b( 2, 1 ) = zero
238 wi = zero
239
240 ELSE IF( abs( b( 2, 2 ) ).LE.ulp ) THEN
241 CALL slartg( a( 2, 2 ), a( 2, 1 ), csr, snr, t )
242 snr = -snr
243 CALL srot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
244 CALL srot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
245 csl = one
246 snl = zero
247 a( 2, 1 ) = zero
248 b( 2, 1 ) = zero
249 b( 2, 2 ) = zero
250 wi = zero
251
252 ELSE
253
254
255
256 CALL slag2( a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2,
257 $ wi )
258
259 IF( wi.EQ.zero ) THEN
260
261
262
263 h1 = scale1*a( 1, 1 ) - wr1*b( 1, 1 )
264 h2 = scale1*a( 1, 2 ) - wr1*b( 1, 2 )
265 h3 = scale1*a( 2, 2 ) - wr1*b( 2, 2 )
266
268 qq =
slapy2( scale1*a( 2, 1 ), h3 )
269
270 IF( rr.GT.qq ) THEN
271
272
273
274
275 CALL slartg( h2, h1, csr, snr, t )
276
277 ELSE
278
279
280
281
282 CALL slartg( h3, scale1*a( 2, 1 ), csr, snr, t )
283
284 END IF
285
286 snr = -snr
287 CALL srot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
288 CALL srot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
289
290
291
292 h1 = max( abs( a( 1, 1 ) )+abs( a( 1, 2 ) ),
293 $ abs( a( 2, 1 ) )+abs( a( 2, 2 ) ) )
294 h2 = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
295 $ abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
296
297 IF( ( scale1*h1 ).GE.abs( wr1 )*h2 ) THEN
298
299
300
301 CALL slartg( b( 1, 1 ), b( 2, 1 ), csl, snl, r )
302
303 ELSE
304
305
306
307 CALL slartg( a( 1, 1 ), a( 2, 1 ), csl, snl, r )
308
309 END IF
310
311 CALL srot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
312 CALL srot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
313
314 a( 2, 1 ) = zero
315 b( 2, 1 ) = zero
316
317 ELSE
318
319
320
321
322 CALL slasv2( b( 1, 1 ), b( 1, 2 ), b( 2, 2 ), r, t, snr,
323 $ csr, snl, csl )
324
325
326
327
328 CALL srot( 2, a( 1, 1 ), lda, a( 2, 1 ), lda, csl, snl )
329 CALL srot( 2, b( 1, 1 ), ldb, b( 2, 1 ), ldb, csl, snl )
330 CALL srot( 2, a( 1, 1 ), 1, a( 1, 2 ), 1, csr, snr )
331 CALL srot( 2, b( 1, 1 ), 1, b( 1, 2 ), 1, csr, snr )
332
333 b( 2, 1 ) = zero
334 b( 1, 2 ) = zero
335
336 END IF
337
338 END IF
339
340
341
342 a( 1, 1 ) = anorm*a( 1, 1 )
343 a( 2, 1 ) = anorm*a( 2, 1 )
344 a( 1, 2 ) = anorm*a( 1, 2 )
345 a( 2, 2 ) = anorm*a( 2, 2 )
346 b( 1, 1 ) = bnorm*b( 1, 1 )
347 b( 2, 1 ) = bnorm*b( 2, 1 )
348 b( 1, 2 ) = bnorm*b( 1, 2 )
349 b( 2, 2 ) = bnorm*b( 2, 2 )
350
351 IF( wi.EQ.zero ) THEN
352 alphar( 1 ) = a( 1, 1 )
353 alphar( 2 ) = a( 2, 2 )
354 alphai( 1 ) = zero
355 alphai( 2 ) = zero
356 beta( 1 ) = b( 1, 1 )
357 beta( 2 ) = b( 2, 2 )
358 ELSE
359 alphar( 1 ) = anorm*wr1 / scale1 / bnorm
360 alphai( 1 ) = anorm*wi / scale1 / bnorm
361 alphar( 2 ) = alphar( 1 )
362 alphai( 2 ) = -alphai( 1 )
363 beta( 1 ) = one
364 beta( 2 ) = one
365 END IF
366
367 RETURN
368
369
370
subroutine slag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
real function slamch(cmach)
SLAMCH
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine slasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT