1*> \brief \b SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLALN2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaln2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaln2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaln2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, 22* LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) 23* 24* .. Scalar Arguments .. 25* LOGICAL LTRANS 26* INTEGER INFO, LDA, LDB, LDX, NA, NW 27* REAL CA, D1, D2, SCALE, SMIN, WI, WR, XNORM 28* .. 29* .. Array Arguments .. 30* REAL A( LDA, * ), B( LDB, * ), X( LDX, * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SLALN2 solves a system of the form (ca A - w D ) X = s B 40*> or (ca A**T - w D) X = s B with possible scaling ("s") and 41*> perturbation of A. (A**T means A-transpose.) 42*> 43*> A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA 44*> real diagonal matrix, w is a real or complex value, and X and B are 45*> NA x 1 matrices -- real if w is real, complex if w is complex. NA 46*> may be 1 or 2. 47*> 48*> If w is complex, X and B are represented as NA x 2 matrices, 49*> the first column of each being the real part and the second 50*> being the imaginary part. 51*> 52*> "s" is a scaling factor (<= 1), computed by SLALN2, which is 53*> so chosen that X can be computed without overflow. X is further 54*> scaled if necessary to assure that norm(ca A - w D)*norm(X) is less 55*> than overflow. 56*> 57*> If both singular values of (ca A - w D) are less than SMIN, 58*> SMIN*identity will be used instead of (ca A - w D). If only one 59*> singular value is less than SMIN, one element of (ca A - w D) will be 60*> perturbed enough to make the smallest singular value roughly SMIN. 61*> If both singular values are at least SMIN, (ca A - w D) will not be 62*> perturbed. In any case, the perturbation will be at most some small 63*> multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values 64*> are computed by infinity-norm approximations, and thus will only be 65*> correct to a factor of 2 or so. 66*> 67*> Note: all input quantities are assumed to be smaller than overflow 68*> by a reasonable factor. (See BIGNUM.) 69*> \endverbatim 70* 71* Arguments: 72* ========== 73* 74*> \param[in] LTRANS 75*> \verbatim 76*> LTRANS is LOGICAL 77*> =.TRUE.: A-transpose will be used. 78*> =.FALSE.: A will be used (not transposed.) 79*> \endverbatim 80*> 81*> \param[in] NA 82*> \verbatim 83*> NA is INTEGER 84*> The size of the matrix A. It may (only) be 1 or 2. 85*> \endverbatim 86*> 87*> \param[in] NW 88*> \verbatim 89*> NW is INTEGER 90*> 1 if "w" is real, 2 if "w" is complex. It may only be 1 91*> or 2. 92*> \endverbatim 93*> 94*> \param[in] SMIN 95*> \verbatim 96*> SMIN is REAL 97*> The desired lower bound on the singular values of A. This 98*> should be a safe distance away from underflow or overflow, 99*> say, between (underflow/machine precision) and (machine 100*> precision * overflow ). (See BIGNUM and ULP.) 101*> \endverbatim 102*> 103*> \param[in] CA 104*> \verbatim 105*> CA is REAL 106*> The coefficient c, which A is multiplied by. 107*> \endverbatim 108*> 109*> \param[in] A 110*> \verbatim 111*> A is REAL array, dimension (LDA,NA) 112*> The NA x NA matrix A. 113*> \endverbatim 114*> 115*> \param[in] LDA 116*> \verbatim 117*> LDA is INTEGER 118*> The leading dimension of A. It must be at least NA. 119*> \endverbatim 120*> 121*> \param[in] D1 122*> \verbatim 123*> D1 is REAL 124*> The 1,1 element in the diagonal matrix D. 125*> \endverbatim 126*> 127*> \param[in] D2 128*> \verbatim 129*> D2 is REAL 130*> The 2,2 element in the diagonal matrix D. Not used if NA=1. 131*> \endverbatim 132*> 133*> \param[in] B 134*> \verbatim 135*> B is REAL array, dimension (LDB,NW) 136*> The NA x NW matrix B (right-hand side). If NW=2 ("w" is 137*> complex), column 1 contains the real part of B and column 2 138*> contains the imaginary part. 139*> \endverbatim 140*> 141*> \param[in] LDB 142*> \verbatim 143*> LDB is INTEGER 144*> The leading dimension of B. It must be at least NA. 145*> \endverbatim 146*> 147*> \param[in] WR 148*> \verbatim 149*> WR is REAL 150*> The real part of the scalar "w". 151*> \endverbatim 152*> 153*> \param[in] WI 154*> \verbatim 155*> WI is REAL 156*> The imaginary part of the scalar "w". Not used if NW=1. 157*> \endverbatim 158*> 159*> \param[out] X 160*> \verbatim 161*> X is REAL array, dimension (LDX,NW) 162*> The NA x NW matrix X (unknowns), as computed by SLALN2. 163*> If NW=2 ("w" is complex), on exit, column 1 will contain 164*> the real part of X and column 2 will contain the imaginary 165*> part. 166*> \endverbatim 167*> 168*> \param[in] LDX 169*> \verbatim 170*> LDX is INTEGER 171*> The leading dimension of X. It must be at least NA. 172*> \endverbatim 173*> 174*> \param[out] SCALE 175*> \verbatim 176*> SCALE is REAL 177*> The scale factor that B must be multiplied by to insure 178*> that overflow does not occur when computing X. Thus, 179*> (ca A - w D) X will be SCALE*B, not B (ignoring 180*> perturbations of A.) It will be at most 1. 181*> \endverbatim 182*> 183*> \param[out] XNORM 184*> \verbatim 185*> XNORM is REAL 186*> The infinity-norm of X, when X is regarded as an NA x NW 187*> real matrix. 188*> \endverbatim 189*> 190*> \param[out] INFO 191*> \verbatim 192*> INFO is INTEGER 193*> An error flag. It will be set to zero if no error occurs, 194*> a negative number if an argument is in error, or a positive 195*> number if ca A - w D had to be perturbed. 196*> The possible values are: 197*> = 0: No error occurred, and (ca A - w D) did not have to be 198*> perturbed. 199*> = 1: (ca A - w D) had to be perturbed to make its smallest 200*> (or only) singular value greater than SMIN. 201*> NOTE: In the interests of speed, this routine does not 202*> check the inputs for errors. 203*> \endverbatim 204* 205* Authors: 206* ======== 207* 208*> \author Univ. of Tennessee 209*> \author Univ. of California Berkeley 210*> \author Univ. of Colorado Denver 211*> \author NAG Ltd. 212* 213*> \ingroup realOTHERauxiliary 214* 215* ===================================================================== 216 SUBROUTINE SLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, 217 $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) 218* 219* -- LAPACK auxiliary routine -- 220* -- LAPACK is a software package provided by Univ. of Tennessee, -- 221* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 222* 223* .. Scalar Arguments .. 224 LOGICAL LTRANS 225 INTEGER INFO, LDA, LDB, LDX, NA, NW 226 REAL CA, D1, D2, SCALE, SMIN, WI, WR, XNORM 227* .. 228* .. Array Arguments .. 229 REAL A( LDA, * ), B( LDB, * ), X( LDX, * ) 230* .. 231* 232* ===================================================================== 233* 234* .. Parameters .. 235 REAL ZERO, ONE 236 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 237 REAL TWO 238 PARAMETER ( TWO = 2.0E0 ) 239* .. 240* .. Local Scalars .. 241 INTEGER ICMAX, J 242 REAL BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21, 243 $ CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21, 244 $ LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R, 245 $ UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S, 246 $ UR22, XI1, XI2, XR1, XR2 247* .. 248* .. Local Arrays .. 249 LOGICAL CSWAP( 4 ), RSWAP( 4 ) 250 INTEGER IPIVOT( 4, 4 ) 251 REAL CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 ) 252* .. 253* .. External Functions .. 254 REAL SLAMCH 255 EXTERNAL SLAMCH 256* .. 257* .. External Subroutines .. 258 EXTERNAL SLADIV 259* .. 260* .. Intrinsic Functions .. 261 INTRINSIC ABS, MAX 262* .. 263* .. Equivalences .. 264 EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ), 265 $ ( CR( 1, 1 ), CRV( 1 ) ) 266* .. 267* .. Data statements .. 268 DATA CSWAP / .FALSE., .FALSE., .TRUE., .TRUE. / 269 DATA RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. / 270 DATA IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4, 271 $ 3, 2, 1 / 272* .. 273* .. Executable Statements .. 274* 275* Compute BIGNUM 276* 277 SMLNUM = TWO*SLAMCH( 'Safe minimum' ) 278 BIGNUM = ONE / SMLNUM 279 SMINI = MAX( SMIN, SMLNUM ) 280* 281* Don't check for input errors 282* 283 INFO = 0 284* 285* Standard Initializations 286* 287 SCALE = ONE 288* 289 IF( NA.EQ.1 ) THEN 290* 291* 1 x 1 (i.e., scalar) system C X = B 292* 293 IF( NW.EQ.1 ) THEN 294* 295* Real 1x1 system. 296* 297* C = ca A - w D 298* 299 CSR = CA*A( 1, 1 ) - WR*D1 300 CNORM = ABS( CSR ) 301* 302* If | C | < SMINI, use C = SMINI 303* 304 IF( CNORM.LT.SMINI ) THEN 305 CSR = SMINI 306 CNORM = SMINI 307 INFO = 1 308 END IF 309* 310* Check scaling for X = B / C 311* 312 BNORM = ABS( B( 1, 1 ) ) 313 IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN 314 IF( BNORM.GT.BIGNUM*CNORM ) 315 $ SCALE = ONE / BNORM 316 END IF 317* 318* Compute X 319* 320 X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR 321 XNORM = ABS( X( 1, 1 ) ) 322 ELSE 323* 324* Complex 1x1 system (w is complex) 325* 326* C = ca A - w D 327* 328 CSR = CA*A( 1, 1 ) - WR*D1 329 CSI = -WI*D1 330 CNORM = ABS( CSR ) + ABS( CSI ) 331* 332* If | C | < SMINI, use C = SMINI 333* 334 IF( CNORM.LT.SMINI ) THEN 335 CSR = SMINI 336 CSI = ZERO 337 CNORM = SMINI 338 INFO = 1 339 END IF 340* 341* Check scaling for X = B / C 342* 343 BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) ) 344 IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN 345 IF( BNORM.GT.BIGNUM*CNORM ) 346 $ SCALE = ONE / BNORM 347 END IF 348* 349* Compute X 350* 351 CALL SLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI, 352 $ X( 1, 1 ), X( 1, 2 ) ) 353 XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) 354 END IF 355* 356 ELSE 357* 358* 2x2 System 359* 360* Compute the real part of C = ca A - w D (or ca A**T - w D ) 361* 362 CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1 363 CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2 364 IF( LTRANS ) THEN 365 CR( 1, 2 ) = CA*A( 2, 1 ) 366 CR( 2, 1 ) = CA*A( 1, 2 ) 367 ELSE 368 CR( 2, 1 ) = CA*A( 2, 1 ) 369 CR( 1, 2 ) = CA*A( 1, 2 ) 370 END IF 371* 372 IF( NW.EQ.1 ) THEN 373* 374* Real 2x2 system (w is real) 375* 376* Find the largest element in C 377* 378 CMAX = ZERO 379 ICMAX = 0 380* 381 DO 10 J = 1, 4 382 IF( ABS( CRV( J ) ).GT.CMAX ) THEN 383 CMAX = ABS( CRV( J ) ) 384 ICMAX = J 385 END IF 386 10 CONTINUE 387* 388* If norm(C) < SMINI, use SMINI*identity. 389* 390 IF( CMAX.LT.SMINI ) THEN 391 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) ) 392 IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN 393 IF( BNORM.GT.BIGNUM*SMINI ) 394 $ SCALE = ONE / BNORM 395 END IF 396 TEMP = SCALE / SMINI 397 X( 1, 1 ) = TEMP*B( 1, 1 ) 398 X( 2, 1 ) = TEMP*B( 2, 1 ) 399 XNORM = TEMP*BNORM 400 INFO = 1 401 RETURN 402 END IF 403* 404* Gaussian elimination with complete pivoting. 405* 406 UR11 = CRV( ICMAX ) 407 CR21 = CRV( IPIVOT( 2, ICMAX ) ) 408 UR12 = CRV( IPIVOT( 3, ICMAX ) ) 409 CR22 = CRV( IPIVOT( 4, ICMAX ) ) 410 UR11R = ONE / UR11 411 LR21 = UR11R*CR21 412 UR22 = CR22 - UR12*LR21 413* 414* If smaller pivot < SMINI, use SMINI 415* 416 IF( ABS( UR22 ).LT.SMINI ) THEN 417 UR22 = SMINI 418 INFO = 1 419 END IF 420 IF( RSWAP( ICMAX ) ) THEN 421 BR1 = B( 2, 1 ) 422 BR2 = B( 1, 1 ) 423 ELSE 424 BR1 = B( 1, 1 ) 425 BR2 = B( 2, 1 ) 426 END IF 427 BR2 = BR2 - LR21*BR1 428 BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) ) 429 IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN 430 IF( BBND.GE.BIGNUM*ABS( UR22 ) ) 431 $ SCALE = ONE / BBND 432 END IF 433* 434 XR2 = ( BR2*SCALE ) / UR22 435 XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 ) 436 IF( CSWAP( ICMAX ) ) THEN 437 X( 1, 1 ) = XR2 438 X( 2, 1 ) = XR1 439 ELSE 440 X( 1, 1 ) = XR1 441 X( 2, 1 ) = XR2 442 END IF 443 XNORM = MAX( ABS( XR1 ), ABS( XR2 ) ) 444* 445* Further scaling if norm(A) norm(X) > overflow 446* 447 IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN 448 IF( XNORM.GT.BIGNUM / CMAX ) THEN 449 TEMP = CMAX / BIGNUM 450 X( 1, 1 ) = TEMP*X( 1, 1 ) 451 X( 2, 1 ) = TEMP*X( 2, 1 ) 452 XNORM = TEMP*XNORM 453 SCALE = TEMP*SCALE 454 END IF 455 END IF 456 ELSE 457* 458* Complex 2x2 system (w is complex) 459* 460* Find the largest element in C 461* 462 CI( 1, 1 ) = -WI*D1 463 CI( 2, 1 ) = ZERO 464 CI( 1, 2 ) = ZERO 465 CI( 2, 2 ) = -WI*D2 466 CMAX = ZERO 467 ICMAX = 0 468* 469 DO 20 J = 1, 4 470 IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN 471 CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) ) 472 ICMAX = J 473 END IF 474 20 CONTINUE 475* 476* If norm(C) < SMINI, use SMINI*identity. 477* 478 IF( CMAX.LT.SMINI ) THEN 479 BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), 480 $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) 481 IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN 482 IF( BNORM.GT.BIGNUM*SMINI ) 483 $ SCALE = ONE / BNORM 484 END IF 485 TEMP = SCALE / SMINI 486 X( 1, 1 ) = TEMP*B( 1, 1 ) 487 X( 2, 1 ) = TEMP*B( 2, 1 ) 488 X( 1, 2 ) = TEMP*B( 1, 2 ) 489 X( 2, 2 ) = TEMP*B( 2, 2 ) 490 XNORM = TEMP*BNORM 491 INFO = 1 492 RETURN 493 END IF 494* 495* Gaussian elimination with complete pivoting. 496* 497 UR11 = CRV( ICMAX ) 498 UI11 = CIV( ICMAX ) 499 CR21 = CRV( IPIVOT( 2, ICMAX ) ) 500 CI21 = CIV( IPIVOT( 2, ICMAX ) ) 501 UR12 = CRV( IPIVOT( 3, ICMAX ) ) 502 UI12 = CIV( IPIVOT( 3, ICMAX ) ) 503 CR22 = CRV( IPIVOT( 4, ICMAX ) ) 504 CI22 = CIV( IPIVOT( 4, ICMAX ) ) 505 IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN 506* 507* Code when off-diagonals of pivoted C are real 508* 509 IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN 510 TEMP = UI11 / UR11 511 UR11R = ONE / ( UR11*( ONE+TEMP**2 ) ) 512 UI11R = -TEMP*UR11R 513 ELSE 514 TEMP = UR11 / UI11 515 UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) ) 516 UR11R = -TEMP*UI11R 517 END IF 518 LR21 = CR21*UR11R 519 LI21 = CR21*UI11R 520 UR12S = UR12*UR11R 521 UI12S = UR12*UI11R 522 UR22 = CR22 - UR12*LR21 523 UI22 = CI22 - UR12*LI21 524 ELSE 525* 526* Code when diagonals of pivoted C are real 527* 528 UR11R = ONE / UR11 529 UI11R = ZERO 530 LR21 = CR21*UR11R 531 LI21 = CI21*UR11R 532 UR12S = UR12*UR11R 533 UI12S = UI12*UR11R 534 UR22 = CR22 - UR12*LR21 + UI12*LI21 535 UI22 = -UR12*LI21 - UI12*LR21 536 END IF 537 U22ABS = ABS( UR22 ) + ABS( UI22 ) 538* 539* If smaller pivot < SMINI, use SMINI 540* 541 IF( U22ABS.LT.SMINI ) THEN 542 UR22 = SMINI 543 UI22 = ZERO 544 INFO = 1 545 END IF 546 IF( RSWAP( ICMAX ) ) THEN 547 BR2 = B( 1, 1 ) 548 BR1 = B( 2, 1 ) 549 BI2 = B( 1, 2 ) 550 BI1 = B( 2, 2 ) 551 ELSE 552 BR1 = B( 1, 1 ) 553 BR2 = B( 2, 1 ) 554 BI1 = B( 1, 2 ) 555 BI2 = B( 2, 2 ) 556 END IF 557 BR2 = BR2 - LR21*BR1 + LI21*BI1 558 BI2 = BI2 - LI21*BR1 - LR21*BI1 559 BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )* 560 $ ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ), 561 $ ABS( BR2 )+ABS( BI2 ) ) 562 IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN 563 IF( BBND.GE.BIGNUM*U22ABS ) THEN 564 SCALE = ONE / BBND 565 BR1 = SCALE*BR1 566 BI1 = SCALE*BI1 567 BR2 = SCALE*BR2 568 BI2 = SCALE*BI2 569 END IF 570 END IF 571* 572 CALL SLADIV( BR2, BI2, UR22, UI22, XR2, XI2 ) 573 XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2 574 XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2 575 IF( CSWAP( ICMAX ) ) THEN 576 X( 1, 1 ) = XR2 577 X( 2, 1 ) = XR1 578 X( 1, 2 ) = XI2 579 X( 2, 2 ) = XI1 580 ELSE 581 X( 1, 1 ) = XR1 582 X( 2, 1 ) = XR2 583 X( 1, 2 ) = XI1 584 X( 2, 2 ) = XI2 585 END IF 586 XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) ) 587* 588* Further scaling if norm(A) norm(X) > overflow 589* 590 IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN 591 IF( XNORM.GT.BIGNUM / CMAX ) THEN 592 TEMP = CMAX / BIGNUM 593 X( 1, 1 ) = TEMP*X( 1, 1 ) 594 X( 2, 1 ) = TEMP*X( 2, 1 ) 595 X( 1, 2 ) = TEMP*X( 1, 2 ) 596 X( 2, 2 ) = TEMP*X( 2, 2 ) 597 XNORM = TEMP*XNORM 598 SCALE = TEMP*SCALE 599 END IF 600 END IF 601 END IF 602 END IF 603* 604 RETURN 605* 606* End of SLALN2 607* 608 END 609