1*> \brief \b SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLASY2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasy2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasy2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasy2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, 22* LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO ) 23* 24* .. Scalar Arguments .. 25* LOGICAL LTRANL, LTRANR 26* INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2 27* REAL SCALE, XNORM 28* .. 29* .. Array Arguments .. 30* REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ), 31* $ X( LDX, * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in 41*> 42*> op(TL)*X + ISGN*X*op(TR) = SCALE*B, 43*> 44*> where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or 45*> -1. op(T) = T or T**T, where T**T denotes the transpose of T. 46*> \endverbatim 47* 48* Arguments: 49* ========== 50* 51*> \param[in] LTRANL 52*> \verbatim 53*> LTRANL is LOGICAL 54*> On entry, LTRANL specifies the op(TL): 55*> = .FALSE., op(TL) = TL, 56*> = .TRUE., op(TL) = TL**T. 57*> \endverbatim 58*> 59*> \param[in] LTRANR 60*> \verbatim 61*> LTRANR is LOGICAL 62*> On entry, LTRANR specifies the op(TR): 63*> = .FALSE., op(TR) = TR, 64*> = .TRUE., op(TR) = TR**T. 65*> \endverbatim 66*> 67*> \param[in] ISGN 68*> \verbatim 69*> ISGN is INTEGER 70*> On entry, ISGN specifies the sign of the equation 71*> as described before. ISGN may only be 1 or -1. 72*> \endverbatim 73*> 74*> \param[in] N1 75*> \verbatim 76*> N1 is INTEGER 77*> On entry, N1 specifies the order of matrix TL. 78*> N1 may only be 0, 1 or 2. 79*> \endverbatim 80*> 81*> \param[in] N2 82*> \verbatim 83*> N2 is INTEGER 84*> On entry, N2 specifies the order of matrix TR. 85*> N2 may only be 0, 1 or 2. 86*> \endverbatim 87*> 88*> \param[in] TL 89*> \verbatim 90*> TL is REAL array, dimension (LDTL,2) 91*> On entry, TL contains an N1 by N1 matrix. 92*> \endverbatim 93*> 94*> \param[in] LDTL 95*> \verbatim 96*> LDTL is INTEGER 97*> The leading dimension of the matrix TL. LDTL >= max(1,N1). 98*> \endverbatim 99*> 100*> \param[in] TR 101*> \verbatim 102*> TR is REAL array, dimension (LDTR,2) 103*> On entry, TR contains an N2 by N2 matrix. 104*> \endverbatim 105*> 106*> \param[in] LDTR 107*> \verbatim 108*> LDTR is INTEGER 109*> The leading dimension of the matrix TR. LDTR >= max(1,N2). 110*> \endverbatim 111*> 112*> \param[in] B 113*> \verbatim 114*> B is REAL array, dimension (LDB,2) 115*> On entry, the N1 by N2 matrix B contains the right-hand 116*> side of the equation. 117*> \endverbatim 118*> 119*> \param[in] LDB 120*> \verbatim 121*> LDB is INTEGER 122*> The leading dimension of the matrix B. LDB >= max(1,N1). 123*> \endverbatim 124*> 125*> \param[out] SCALE 126*> \verbatim 127*> SCALE is REAL 128*> On exit, SCALE contains the scale factor. SCALE is chosen 129*> less than or equal to 1 to prevent the solution overflowing. 130*> \endverbatim 131*> 132*> \param[out] X 133*> \verbatim 134*> X is REAL array, dimension (LDX,2) 135*> On exit, X contains the N1 by N2 solution. 136*> \endverbatim 137*> 138*> \param[in] LDX 139*> \verbatim 140*> LDX is INTEGER 141*> The leading dimension of the matrix X. LDX >= max(1,N1). 142*> \endverbatim 143*> 144*> \param[out] XNORM 145*> \verbatim 146*> XNORM is REAL 147*> On exit, XNORM is the infinity-norm of the solution. 148*> \endverbatim 149*> 150*> \param[out] INFO 151*> \verbatim 152*> INFO is INTEGER 153*> On exit, INFO is set to 154*> 0: successful exit. 155*> 1: TL and TR have too close eigenvalues, so TL or 156*> TR is perturbed to get a nonsingular equation. 157*> NOTE: In the interests of speed, this routine does not 158*> check the inputs for errors. 159*> \endverbatim 160* 161* Authors: 162* ======== 163* 164*> \author Univ. of Tennessee 165*> \author Univ. of California Berkeley 166*> \author Univ. of Colorado Denver 167*> \author NAG Ltd. 168* 169*> \date September 2012 170* 171*> \ingroup realSYauxiliary 172* 173* ===================================================================== 174 SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, 175 $ LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO ) 176* 177* -- LAPACK auxiliary routine (version 3.4.2) -- 178* -- LAPACK is a software package provided by Univ. of Tennessee, -- 179* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 180* September 2012 181* 182* .. Scalar Arguments .. 183 LOGICAL LTRANL, LTRANR 184 INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2 185 REAL SCALE, XNORM 186* .. 187* .. Array Arguments .. 188 REAL B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ), 189 $ X( LDX, * ) 190* .. 191* 192* ===================================================================== 193* 194* .. Parameters .. 195 REAL ZERO, ONE 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 ) 199* .. 200* .. Local Scalars .. 201 LOGICAL BSWAP, XSWAP 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 205* .. 206* .. Local Arrays .. 207 LOGICAL BSWPIV( 4 ), XSWPIV( 4 ) 208 INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ), 209 $ LOCU22( 4 ) 210 REAL BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 ) 211* .. 212* .. External Functions .. 213 INTEGER ISAMAX 214 REAL SLAMCH 215 EXTERNAL ISAMAX, SLAMCH 216* .. 217* .. External Subroutines .. 218 EXTERNAL SCOPY, SSWAP 219* .. 220* .. Intrinsic Functions .. 221 INTRINSIC ABS, MAX 222* .. 223* .. Data statements .. 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. / 228* .. 229* .. Executable Statements .. 230* 231* Do not check the input parameters for errors 232* 233 INFO = 0 234* 235* Quick return if possible 236* 237 IF( N1.EQ.0 .OR. N2.EQ.0 ) 238 $ RETURN 239* 240* Set constants to control overflow 241* 242 EPS = SLAMCH( 'P' ) 243 SMLNUM = SLAMCH( 'S' ) / EPS 244 SGN = ISGN 245* 246 K = N1 + N1 + N2 - 2 247 GO TO ( 10, 20, 30, 50 )K 248* 249* 1 by 1: TL11*X + SGN*X*TR11 = B11 250* 251 10 CONTINUE 252 TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 ) 253 BET = ABS( TAU1 ) 254 IF( BET.LE.SMLNUM ) THEN 255 TAU1 = SMLNUM 256 BET = SMLNUM 257 INFO = 1 258 END IF 259* 260 SCALE = ONE 261 GAM = ABS( B( 1, 1 ) ) 262 IF( SMLNUM*GAM.GT.BET ) 263 $ SCALE = ONE / GAM 264* 265 X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1 266 XNORM = ABS( X( 1, 1 ) ) 267 RETURN 268* 269* 1 by 2: 270* TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12] 271* [TR21 TR22] 272* 273 20 CONTINUE 274* 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 ) ) ), 277 $ SMLNUM ) 278 TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 ) 279 TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 ) 280 IF( LTRANR ) THEN 281 TMP( 2 ) = SGN*TR( 2, 1 ) 282 TMP( 3 ) = SGN*TR( 1, 2 ) 283 ELSE 284 TMP( 2 ) = SGN*TR( 1, 2 ) 285 TMP( 3 ) = SGN*TR( 2, 1 ) 286 END IF 287 BTMP( 1 ) = B( 1, 1 ) 288 BTMP( 2 ) = B( 1, 2 ) 289 GO TO 40 290* 291* 2 by 1: 292* op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11] 293* [TL21 TL22] [X21] [X21] [B21] 294* 295 30 CONTINUE 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 ) ) ), 298 $ SMLNUM ) 299 TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 ) 300 TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 ) 301 IF( LTRANL ) THEN 302 TMP( 2 ) = TL( 1, 2 ) 303 TMP( 3 ) = TL( 2, 1 ) 304 ELSE 305 TMP( 2 ) = TL( 2, 1 ) 306 TMP( 3 ) = TL( 1, 2 ) 307 END IF 308 BTMP( 1 ) = B( 1, 1 ) 309 BTMP( 2 ) = B( 2, 1 ) 310 40 CONTINUE 311* 312* Solve 2 by 2 system using complete pivoting. 313* Set pivots less than SMIN to SMIN. 314* 315 IPIV = ISAMAX( 4, TMP, 1 ) 316 U11 = TMP( IPIV ) 317 IF( ABS( U11 ).LE.SMIN ) THEN 318 INFO = 1 319 U11 = SMIN 320 END IF 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 327 INFO = 1 328 U22 = SMIN 329 END IF 330 IF( BSWAP ) THEN 331 TEMP = BTMP( 2 ) 332 BTMP( 2 ) = BTMP( 1 ) - L21*TEMP 333 BTMP( 1 ) = TEMP 334 ELSE 335 BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 ) 336 END IF 337 SCALE = ONE 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 343 END IF 344 X2( 2 ) = BTMP( 2 ) / U22 345 X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 ) 346 IF( XSWAP ) THEN 347 TEMP = X2( 2 ) 348 X2( 2 ) = X2( 1 ) 349 X2( 1 ) = TEMP 350 END IF 351 X( 1, 1 ) = X2( 1 ) 352 IF( N1.EQ.1 ) THEN 353 X( 1, 2 ) = X2( 2 ) 354 XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) 355 ELSE 356 X( 2, 1 ) = X2( 2 ) 357 XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) ) 358 END IF 359 RETURN 360* 361* 2 by 2: 362* op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12] 363* [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22] 364* 365* Solve equivalent 4 by 4 system using complete pivoting. 366* Set pivots less than SMIN to SMIN. 367* 368 50 CONTINUE 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 ) 374 BTMP( 1 ) = ZERO 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 ) 380 IF( LTRANL ) THEN 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 ) 385 ELSE 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 ) 390 END IF 391 IF( LTRANR ) THEN 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 ) 396 ELSE 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 ) 401 END IF 402 BTMP( 1 ) = B( 1, 1 ) 403 BTMP( 2 ) = B( 2, 1 ) 404 BTMP( 3 ) = B( 1, 2 ) 405 BTMP( 4 ) = B( 2, 2 ) 406* 407* Perform elimination 408* 409 DO 100 I = 1, 3 410 XMAX = ZERO 411 DO 70 IP = I, 4 412 DO 60 JP = I, 4 413 IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN 414 XMAX = ABS( T16( IP, JP ) ) 415 IPSV = IP 416 JPSV = JP 417 END IF 418 60 CONTINUE 419 70 CONTINUE 420 IF( IPSV.NE.I ) THEN 421 CALL SSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 ) 422 TEMP = BTMP( I ) 423 BTMP( I ) = BTMP( IPSV ) 424 BTMP( IPSV ) = TEMP 425 END IF 426 IF( JPSV.NE.I ) 427 $ CALL SSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 ) 428 JPIV( I ) = JPSV 429 IF( ABS( T16( I, I ) ).LT.SMIN ) THEN 430 INFO = 1 431 T16( I, I ) = SMIN 432 END IF 433 DO 90 J = I + 1, 4 434 T16( J, I ) = T16( J, I ) / T16( I, I ) 435 BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I ) 436 DO 80 K = I + 1, 4 437 T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K ) 438 80 CONTINUE 439 90 CONTINUE 440 100 CONTINUE 441 IF( ABS( T16( 4, 4 ) ).LT.SMIN ) 442 $ T16( 4, 4 ) = SMIN 443 SCALE = ONE 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 454 END IF 455 DO 120 I = 1, 4 456 K = 5 - I 457 TEMP = ONE / T16( K, K ) 458 TMP( K ) = BTMP( K )*TEMP 459 DO 110 J = K + 1, 4 460 TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J ) 461 110 CONTINUE 462 120 CONTINUE 463 DO 130 I = 1, 3 464 IF( JPIV( 4-I ).NE.4-I ) THEN 465 TEMP = TMP( 4-I ) 466 TMP( 4-I ) = TMP( JPIV( 4-I ) ) 467 TMP( JPIV( 4-I ) ) = TEMP 468 END IF 469 130 CONTINUE 470 X( 1, 1 ) = TMP( 1 ) 471 X( 2, 1 ) = TMP( 2 ) 472 X( 1, 2 ) = TMP( 3 ) 473 X( 2, 2 ) = TMP( 4 ) 474 XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ), 475 $ ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) ) 476 RETURN 477* 478* End of SLASY2 479* 480 END 481c $Id$ 482