174 SUBROUTINE slasy2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
175 $ ldtr, b, ldb, scale, x, ldx, xnorm, info )
183 LOGICAL LTRANL, LTRANR
184 INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
188 REAL B( ldb, * ), TL( ldtl, * ), TR( ldtr, * ),
196 parameter( zero = 0.0e+0, one = 1.0e+0 )
197 REAL TWO, HALF, EIGHT
198 parameter( two = 2.0e+0, half = 0.5e+0, eight = 8.0e+0 )
202 INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
203 REAL BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
204 $ temp, u11, u12, u22, xmax
207 LOGICAL BSWPIV( 4 ), XSWPIV( 4 )
208 INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
210 REAL BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
215 EXTERNAL isamax, slamch
224 DATA locu12 / 3, 4, 1, 2 / , locl21 / 2, 1, 4, 3 / ,
225 $ locu22 / 4, 3, 2, 1 /
226 DATA xswpiv / .false., .false., .true., .true. /
227 DATA bswpiv / .false., .true., .false., .true. /
237 IF( n1.EQ.0 .OR. n2.EQ.0 )
243 smlnum = slamch(
'S' ) / eps
247 GO TO ( 10, 20, 30, 50 )k
252 tau1 = tl( 1, 1 ) + sgn*tr( 1, 1 )
254 IF( bet.LE.smlnum )
THEN
261 gam = abs( b( 1, 1 ) )
262 IF( smlnum*gam.GT.bet )
265 x( 1, 1 ) = ( b( 1, 1 )*scale ) / tau1
266 xnorm = abs( x( 1, 1 ) )
275 smin = max( eps*max( abs( tl( 1, 1 ) ), abs( tr( 1, 1 ) ),
276 $ abs( tr( 1, 2 ) ), abs( tr( 2, 1 ) ), abs( tr( 2, 2 ) ) ),
278 tmp( 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
279 tmp( 4 ) = tl( 1, 1 ) + sgn*tr( 2, 2 )
281 tmp( 2 ) = sgn*tr( 2, 1 )
282 tmp( 3 ) = sgn*tr( 1, 2 )
284 tmp( 2 ) = sgn*tr( 1, 2 )
285 tmp( 3 ) = sgn*tr( 2, 1 )
287 btmp( 1 ) = b( 1, 1 )
288 btmp( 2 ) = b( 1, 2 )
296 smin = max( eps*max( abs( tr( 1, 1 ) ), abs( tl( 1, 1 ) ),
297 $ abs( tl( 1, 2 ) ), abs( tl( 2, 1 ) ), abs( tl( 2, 2 ) ) ),
299 tmp( 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
300 tmp( 4 ) = tl( 2, 2 ) + sgn*tr( 1, 1 )
302 tmp( 2 ) = tl( 1, 2 )
303 tmp( 3 ) = tl( 2, 1 )
305 tmp( 2 ) = tl( 2, 1 )
306 tmp( 3 ) = tl( 1, 2 )
308 btmp( 1 ) = b( 1, 1 )
309 btmp( 2 ) = b( 2, 1 )
315 ipiv = isamax( 4, tmp, 1 )
317 IF( abs( u11 ).LE.smin )
THEN
321 u12 = tmp( locu12( ipiv ) )
322 l21 = tmp( locl21( ipiv ) ) / u11
323 u22 = tmp( locu22( ipiv ) ) - u12*l21
324 xswap = xswpiv( ipiv )
325 bswap = bswpiv( ipiv )
326 IF( abs( u22 ).LE.smin )
THEN
332 btmp( 2 ) = btmp( 1 ) - l21*temp
335 btmp( 2 ) = btmp( 2 ) - l21*btmp( 1 )
338 IF( ( two*smlnum )*abs( btmp( 2 ) ).GT.abs( u22 ) .OR.
339 $ ( two*smlnum )*abs( btmp( 1 ) ).GT.abs( u11 ) )
THEN
340 scale = half / max( abs( btmp( 1 ) ), abs( btmp( 2 ) ) )
341 btmp( 1 ) = btmp( 1 )*scale
342 btmp( 2 ) = btmp( 2 )*scale
344 x2( 2 ) = btmp( 2 ) / u22
345 x2( 1 ) = btmp( 1 ) / u11 - ( u12 / u11 )*x2( 2 )
354 xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
357 xnorm = max( abs( x( 1, 1 ) ), abs( x( 2, 1 ) ) )
369 smin = max( abs( tr( 1, 1 ) ), abs( tr( 1, 2 ) ),
370 $ abs( tr( 2, 1 ) ), abs( tr( 2, 2 ) ) )
371 smin = max( smin, abs( tl( 1, 1 ) ), abs( tl( 1, 2 ) ),
372 $ abs( tl( 2, 1 ) ), abs( tl( 2, 2 ) ) )
373 smin = max( eps*smin, smlnum )
375 CALL scopy( 16, btmp, 0, t16, 1 )
376 t16( 1, 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
377 t16( 2, 2 ) = tl( 2, 2 ) + sgn*tr( 1, 1 )
378 t16( 3, 3 ) = tl( 1, 1 ) + sgn*tr( 2, 2 )
379 t16( 4, 4 ) = tl( 2, 2 ) + sgn*tr( 2, 2 )
381 t16( 1, 2 ) = tl( 2, 1 )
382 t16( 2, 1 ) = tl( 1, 2 )
383 t16( 3, 4 ) = tl( 2, 1 )
384 t16( 4, 3 ) = tl( 1, 2 )
386 t16( 1, 2 ) = tl( 1, 2 )
387 t16( 2, 1 ) = tl( 2, 1 )
388 t16( 3, 4 ) = tl( 1, 2 )
389 t16( 4, 3 ) = tl( 2, 1 )
392 t16( 1, 3 ) = sgn*tr( 1, 2 )
393 t16( 2, 4 ) = sgn*tr( 1, 2 )
394 t16( 3, 1 ) = sgn*tr( 2, 1 )
395 t16( 4, 2 ) = sgn*tr( 2, 1 )
397 t16( 1, 3 ) = sgn*tr( 2, 1 )
398 t16( 2, 4 ) = sgn*tr( 2, 1 )
399 t16( 3, 1 ) = sgn*tr( 1, 2 )
400 t16( 4, 2 ) = sgn*tr( 1, 2 )
402 btmp( 1 ) = b( 1, 1 )
403 btmp( 2 ) = b( 2, 1 )
404 btmp( 3 ) = b( 1, 2 )
405 btmp( 4 ) = b( 2, 2 )
413 IF( abs( t16( ip, jp ) ).GE.xmax )
THEN
414 xmax = abs( t16( ip, jp ) )
421 CALL sswap( 4, t16( ipsv, 1 ), 4, t16( i, 1 ), 4 )
423 btmp( i ) = btmp( ipsv )
427 $
CALL sswap( 4, t16( 1, jpsv ), 1, t16( 1, i ), 1 )
429 IF( abs( t16( i, i ) ).LT.smin )
THEN
434 t16( j, i ) = t16( j, i ) / t16( i, i )
435 btmp( j ) = btmp( j ) - t16( j, i )*btmp( i )
437 t16( j, k ) = t16( j, k ) - t16( j, i )*t16( i, k )
441 IF( abs( t16( 4, 4 ) ).LT.smin )
444 IF( ( eight*smlnum )*abs( btmp( 1 ) ).GT.abs( t16( 1, 1 ) ) .OR.
445 $ ( eight*smlnum )*abs( btmp( 2 ) ).GT.abs( t16( 2, 2 ) ) .OR.
446 $ ( eight*smlnum )*abs( btmp( 3 ) ).GT.abs( t16( 3, 3 ) ) .OR.
447 $ ( eight*smlnum )*abs( btmp( 4 ) ).GT.abs( t16( 4, 4 ) ) )
THEN
448 scale = ( one / eight ) / max( abs( btmp( 1 ) ),
449 $ abs( btmp( 2 ) ), abs( btmp( 3 ) ), abs( btmp( 4 ) ) )
450 btmp( 1 ) = btmp( 1 )*scale
451 btmp( 2 ) = btmp( 2 )*scale
452 btmp( 3 ) = btmp( 3 )*scale
453 btmp( 4 ) = btmp( 4 )*scale
457 temp = one / t16( k, k )
458 tmp( k ) = btmp( k )*temp
460 tmp( k ) = tmp( k ) - ( temp*t16( k, j ) )*tmp( j )
464 IF( jpiv( 4-i ).NE.4-i )
THEN
466 tmp( 4-i ) = tmp( jpiv( 4-i ) )
467 tmp( jpiv( 4-i ) ) = temp
474 xnorm = max( abs( tmp( 1 ) )+abs( tmp( 3 ) ),
475 $ abs( tmp( 2 ) )+abs( tmp( 4 ) ) )
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine slasy2(LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY