SRC/dlasy2.f

   1 *> \brief \b DLASY2 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 DLASY2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasy2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasy2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasy2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASY2( 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 *       DOUBLE PRECISION   SCALE, XNORM
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION   B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
  31 *      $                   X( LDX, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DLASY2 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
 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 DOUBLE PRECISION 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 DOUBLE PRECISION
 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 June 2016
 170 *
 171 *> \ingroup doubleSYauxiliary
 172 *
 173 *  =====================================================================
 174       SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
 175      $                   LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
 176 *
 177 *  -- LAPACK auxiliary routine (version 3.6.1) --
 178 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 179 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 180 *     June 2016
 181 *
 182 *     .. Scalar Arguments ..
 183       LOGICAL            LTRANL, LTRANR
 184       INTEGER            INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
 185       DOUBLE PRECISION   SCALE, XNORM
 186 *     ..
 187 *     .. Array Arguments ..
 188       DOUBLE PRECISION   B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
 189      $                   X( LDX, * )
 190 *     ..
 191 *
 192 * =====================================================================
 193 *
 194 *     .. Parameters ..
 195       DOUBLE PRECISION   ZERO, ONE
 196       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 197       DOUBLE PRECISION   TWO, HALF, EIGHT
 198       PARAMETER          ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
 199 *     ..
 200 *     .. Local Scalars ..
 201       LOGICAL            BSWAP, XSWAP
 202       INTEGER            I, IP, IPIV, IPSV, J, JP, JPSV, K
 203       DOUBLE PRECISION   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       DOUBLE PRECISION   BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
 211 *     ..
 212 *     .. External Functions ..
 213       INTEGER            IDAMAX
 214       DOUBLE PRECISION   DLAMCH
 215       EXTERNAL           IDAMAX, DLAMCH
 216 *     ..
 217 *     .. External Subroutines ..
 218       EXTERNAL           DCOPY, DSWAP
 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 = DLAMCH( 'P' )
 243       SMLNUM = DLAMCH( '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 = IDAMAX( 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 DCOPY( 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 DSWAP( 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 DSWAP( 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 ) THEN
 442          INFO = 1
 443          T16( 4, 4 ) = SMIN
 444       END IF
 445       SCALE = ONE
 446       IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
 447      $    ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
 448      $    ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
 449      $    ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
 450          SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
 451      $           ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
 452          BTMP( 1 ) = BTMP( 1 )*SCALE
 453          BTMP( 2 ) = BTMP( 2 )*SCALE
 454          BTMP( 3 ) = BTMP( 3 )*SCALE
 455          BTMP( 4 ) = BTMP( 4 )*SCALE
 456       END IF
 457       DO 120 I = 1, 4
 458          K = 5 - I
 459          TEMP = ONE / T16( K, K )
 460          TMP( K ) = BTMP( K )*TEMP
 461          DO 110 J = K + 1, 4
 462             TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
 463   110    CONTINUE
 464   120 CONTINUE
 465       DO 130 I = 1, 3
 466          IF( JPIV( 4-I ).NE.4-I ) THEN
 467             TEMP = TMP( 4-I )
 468             TMP( 4-I ) = TMP( JPIV( 4-I ) )
 469             TMP( JPIV( 4-I ) ) = TEMP
 470          END IF
 471   130 CONTINUE
 472       X( 1, 1 ) = TMP( 1 )
 473       X( 2, 1 ) = TMP( 2 )
 474       X( 1, 2 ) = TMP( 3 )
 475       X( 2, 2 ) = TMP( 4 )
 476       XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
 477      $        ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
 478       RETURN
 479 *
 480 *     End of DLASY2
 481 *
 482       END