LAPACK: SRC/slasy2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
00002      $                   LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            LTRANL, LTRANR
00011       INTEGER            INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
00012       REAL               SCALE, XNORM
00013 *     ..
00014 *     .. Array Arguments ..
00015       REAL               B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
00016      $                   X( LDX, * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
00023 *
00024 *         op(TL)*X + ISGN*X*op(TR) = SCALE*B,
00025 *
00026 *  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
00027 *  -1.  op(T) = T or T', where T' denotes the transpose of T.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  LTRANL  (input) LOGICAL
00033 *          On entry, LTRANL specifies the op(TL):
00034 *             = .FALSE., op(TL) = TL,
00035 *             = .TRUE., op(TL) = TL'.
00036 *
00037 *  LTRANR  (input) LOGICAL
00038 *          On entry, LTRANR specifies the op(TR):
00039 *            = .FALSE., op(TR) = TR,
00040 *            = .TRUE., op(TR) = TR'.
00041 *
00042 *  ISGN    (input) INTEGER
00043 *          On entry, ISGN specifies the sign of the equation
00044 *          as described before. ISGN may only be 1 or -1.
00045 *
00046 *  N1      (input) INTEGER
00047 *          On entry, N1 specifies the order of matrix TL.
00048 *          N1 may only be 0, 1 or 2.
00049 *
00050 *  N2      (input) INTEGER
00051 *          On entry, N2 specifies the order of matrix TR.
00052 *          N2 may only be 0, 1 or 2.
00053 *
00054 *  TL      (input) REAL array, dimension (LDTL,2)
00055 *          On entry, TL contains an N1 by N1 matrix.
00056 *
00057 *  LDTL    (input) INTEGER
00058 *          The leading dimension of the matrix TL. LDTL >= max(1,N1).
00059 *
00060 *  TR      (input) REAL array, dimension (LDTR,2)
00061 *          On entry, TR contains an N2 by N2 matrix.
00062 *
00063 *  LDTR    (input) INTEGER
00064 *          The leading dimension of the matrix TR. LDTR >= max(1,N2).
00065 *
00066 *  B       (input) REAL array, dimension (LDB,2)
00067 *          On entry, the N1 by N2 matrix B contains the right-hand
00068 *          side of the equation.
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of the matrix B. LDB >= max(1,N1).
00072 *
00073 *  SCALE   (output) REAL
00074 *          On exit, SCALE contains the scale factor. SCALE is chosen
00075 *          less than or equal to 1 to prevent the solution overflowing.
00076 *
00077 *  X       (output) REAL array, dimension (LDX,2)
00078 *          On exit, X contains the N1 by N2 solution.
00079 *
00080 *  LDX     (input) INTEGER
00081 *          The leading dimension of the matrix X. LDX >= max(1,N1).
00082 *
00083 *  XNORM   (output) REAL
00084 *          On exit, XNORM is the infinity-norm of the solution.
00085 *
00086 *  INFO    (output) INTEGER
00087 *          On exit, INFO is set to
00088 *             0: successful exit.
00089 *             1: TL and TR have too close eigenvalues, so TL or
00090 *                TR is perturbed to get a nonsingular equation.
00091 *          NOTE: In the interests of speed, this routine does not
00092 *                check the inputs for errors.
00093 *
00094 * =====================================================================
00095 *
00096 *     .. Parameters ..
00097       REAL               ZERO, ONE
00098       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00099       REAL               TWO, HALF, EIGHT
00100       PARAMETER          ( TWO = 2.0E+0, HALF = 0.5E+0, EIGHT = 8.0E+0 )
00101 *     ..
00102 *     .. Local Scalars ..
00103       LOGICAL            BSWAP, XSWAP
00104       INTEGER            I, IP, IPIV, IPSV, J, JP, JPSV, K
00105       REAL               BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
00106      $                   TEMP, U11, U12, U22, XMAX
00107 *     ..
00108 *     .. Local Arrays ..
00109       LOGICAL            BSWPIV( 4 ), XSWPIV( 4 )
00110       INTEGER            JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
00111      $                   LOCU22( 4 )
00112       REAL               BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
00113 *     ..
00114 *     .. External Functions ..
00115       INTEGER            ISAMAX
00116       REAL               SLAMCH
00117       EXTERNAL           ISAMAX, SLAMCH
00118 *     ..
00119 *     .. External Subroutines ..
00120       EXTERNAL           SCOPY, SSWAP
00121 *     ..
00122 *     .. Intrinsic Functions ..
00123       INTRINSIC          ABS, MAX
00124 *     ..
00125 *     .. Data statements ..
00126       DATA               LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / ,
00127      $                   LOCU22 / 4, 3, 2, 1 /
00128       DATA               XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
00129       DATA               BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
00130 *     ..
00131 *     .. Executable Statements ..
00132 *
00133 *     Do not check the input parameters for errors
00134 *
00135       INFO = 0
00136 *
00137 *     Quick return if possible
00138 *
00139       IF( N1.EQ.0 .OR. N2.EQ.0 )
00140      $   RETURN
00141 *
00142 *     Set constants to control overflow
00143 *
00144       EPS = SLAMCH( 'P' )
00145       SMLNUM = SLAMCH( 'S' ) / EPS
00146       SGN = ISGN
00147 *
00148       K = N1 + N1 + N2 - 2
00149       GO TO ( 10, 20, 30, 50 )K
00150 *
00151 *     1 by 1: TL11*X + SGN*X*TR11 = B11
00152 *
00153    10 CONTINUE
00154       TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
00155       BET = ABS( TAU1 )
00156       IF( BET.LE.SMLNUM ) THEN
00157          TAU1 = SMLNUM
00158          BET = SMLNUM
00159          INFO = 1
00160       END IF
00161 *
00162       SCALE = ONE
00163       GAM = ABS( B( 1, 1 ) )
00164       IF( SMLNUM*GAM.GT.BET )
00165      $   SCALE = ONE / GAM
00166 *
00167       X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
00168       XNORM = ABS( X( 1, 1 ) )
00169       RETURN
00170 *
00171 *     1 by 2:
00172 *     TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12]  = [B11 B12]
00173 *                                       [TR21 TR22]
00174 *
00175    20 CONTINUE
00176 *
00177       SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
00178      $       ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
00179      $       SMLNUM )
00180       TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
00181       TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
00182       IF( LTRANR ) THEN
00183          TMP( 2 ) = SGN*TR( 2, 1 )
00184          TMP( 3 ) = SGN*TR( 1, 2 )
00185       ELSE
00186          TMP( 2 ) = SGN*TR( 1, 2 )
00187          TMP( 3 ) = SGN*TR( 2, 1 )
00188       END IF
00189       BTMP( 1 ) = B( 1, 1 )
00190       BTMP( 2 ) = B( 1, 2 )
00191       GO TO 40
00192 *
00193 *     2 by 1:
00194 *          op[TL11 TL12]*[X11] + ISGN* [X11]*TR11  = [B11]
00195 *            [TL21 TL22] [X21]         [X21]         [B21]
00196 *
00197    30 CONTINUE
00198       SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
00199      $       ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
00200      $       SMLNUM )
00201       TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
00202       TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
00203       IF( LTRANL ) THEN
00204          TMP( 2 ) = TL( 1, 2 )
00205          TMP( 3 ) = TL( 2, 1 )
00206       ELSE
00207          TMP( 2 ) = TL( 2, 1 )
00208          TMP( 3 ) = TL( 1, 2 )
00209       END IF
00210       BTMP( 1 ) = B( 1, 1 )
00211       BTMP( 2 ) = B( 2, 1 )
00212    40 CONTINUE
00213 *
00214 *     Solve 2 by 2 system using complete pivoting.
00215 *     Set pivots less than SMIN to SMIN.
00216 *
00217       IPIV = ISAMAX( 4, TMP, 1 )
00218       U11 = TMP( IPIV )
00219       IF( ABS( U11 ).LE.SMIN ) THEN
00220          INFO = 1
00221          U11 = SMIN
00222       END IF
00223       U12 = TMP( LOCU12( IPIV ) )
00224       L21 = TMP( LOCL21( IPIV ) ) / U11
00225       U22 = TMP( LOCU22( IPIV ) ) - U12*L21
00226       XSWAP = XSWPIV( IPIV )
00227       BSWAP = BSWPIV( IPIV )
00228       IF( ABS( U22 ).LE.SMIN ) THEN
00229          INFO = 1
00230          U22 = SMIN
00231       END IF
00232       IF( BSWAP ) THEN
00233          TEMP = BTMP( 2 )
00234          BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
00235          BTMP( 1 ) = TEMP
00236       ELSE
00237          BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
00238       END IF
00239       SCALE = ONE
00240       IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
00241      $    ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
00242          SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
00243          BTMP( 1 ) = BTMP( 1 )*SCALE
00244          BTMP( 2 ) = BTMP( 2 )*SCALE
00245       END IF
00246       X2( 2 ) = BTMP( 2 ) / U22
00247       X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
00248       IF( XSWAP ) THEN
00249          TEMP = X2( 2 )
00250          X2( 2 ) = X2( 1 )
00251          X2( 1 ) = TEMP
00252       END IF
00253       X( 1, 1 ) = X2( 1 )
00254       IF( N1.EQ.1 ) THEN
00255          X( 1, 2 ) = X2( 2 )
00256          XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
00257       ELSE
00258          X( 2, 1 ) = X2( 2 )
00259          XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
00260       END IF
00261       RETURN
00262 *
00263 *     2 by 2:
00264 *     op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
00265 *       [TL21 TL22] [X21 X22]        [X21 X22]   [TR21 TR22]   [B21 B22]
00266 *
00267 *     Solve equivalent 4 by 4 system using complete pivoting.
00268 *     Set pivots less than SMIN to SMIN.
00269 *
00270    50 CONTINUE
00271       SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
00272      $       ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
00273       SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
00274      $       ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
00275       SMIN = MAX( EPS*SMIN, SMLNUM )
00276       BTMP( 1 ) = ZERO
00277       CALL SCOPY( 16, BTMP, 0, T16, 1 )
00278       T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
00279       T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
00280       T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
00281       T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
00282       IF( LTRANL ) THEN
00283          T16( 1, 2 ) = TL( 2, 1 )
00284          T16( 2, 1 ) = TL( 1, 2 )
00285          T16( 3, 4 ) = TL( 2, 1 )
00286          T16( 4, 3 ) = TL( 1, 2 )
00287       ELSE
00288          T16( 1, 2 ) = TL( 1, 2 )
00289          T16( 2, 1 ) = TL( 2, 1 )
00290          T16( 3, 4 ) = TL( 1, 2 )
00291          T16( 4, 3 ) = TL( 2, 1 )
00292       END IF
00293       IF( LTRANR ) THEN
00294          T16( 1, 3 ) = SGN*TR( 1, 2 )
00295          T16( 2, 4 ) = SGN*TR( 1, 2 )
00296          T16( 3, 1 ) = SGN*TR( 2, 1 )
00297          T16( 4, 2 ) = SGN*TR( 2, 1 )
00298       ELSE
00299          T16( 1, 3 ) = SGN*TR( 2, 1 )
00300          T16( 2, 4 ) = SGN*TR( 2, 1 )
00301          T16( 3, 1 ) = SGN*TR( 1, 2 )
00302          T16( 4, 2 ) = SGN*TR( 1, 2 )
00303       END IF
00304       BTMP( 1 ) = B( 1, 1 )
00305       BTMP( 2 ) = B( 2, 1 )
00306       BTMP( 3 ) = B( 1, 2 )
00307       BTMP( 4 ) = B( 2, 2 )
00308 *
00309 *     Perform elimination
00310 *
00311       DO 100 I = 1, 3
00312          XMAX = ZERO
00313          DO 70 IP = I, 4
00314             DO 60 JP = I, 4
00315                IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
00316                   XMAX = ABS( T16( IP, JP ) )
00317                   IPSV = IP
00318                   JPSV = JP
00319                END IF
00320    60       CONTINUE
00321    70    CONTINUE
00322          IF( IPSV.NE.I ) THEN
00323             CALL SSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
00324             TEMP = BTMP( I )
00325             BTMP( I ) = BTMP( IPSV )
00326             BTMP( IPSV ) = TEMP
00327          END IF
00328          IF( JPSV.NE.I )
00329      $      CALL SSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
00330          JPIV( I ) = JPSV
00331          IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
00332             INFO = 1
00333             T16( I, I ) = SMIN
00334          END IF
00335          DO 90 J = I + 1, 4
00336             T16( J, I ) = T16( J, I ) / T16( I, I )
00337             BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
00338             DO 80 K = I + 1, 4
00339                T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
00340    80       CONTINUE
00341    90    CONTINUE
00342   100 CONTINUE
00343       IF( ABS( T16( 4, 4 ) ).LT.SMIN )
00344      $   T16( 4, 4 ) = SMIN
00345       SCALE = ONE
00346       IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
00347      $    ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
00348      $    ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
00349      $    ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
00350          SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
00351      $           ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
00352          BTMP( 1 ) = BTMP( 1 )*SCALE
00353          BTMP( 2 ) = BTMP( 2 )*SCALE
00354          BTMP( 3 ) = BTMP( 3 )*SCALE
00355          BTMP( 4 ) = BTMP( 4 )*SCALE
00356       END IF
00357       DO 120 I = 1, 4
00358          K = 5 - I
00359          TEMP = ONE / T16( K, K )
00360          TMP( K ) = BTMP( K )*TEMP
00361          DO 110 J = K + 1, 4
00362             TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
00363   110    CONTINUE
00364   120 CONTINUE
00365       DO 130 I = 1, 3
00366          IF( JPIV( 4-I ).NE.4-I ) THEN
00367             TEMP = TMP( 4-I )
00368             TMP( 4-I ) = TMP( JPIV( 4-I ) )
00369             TMP( JPIV( 4-I ) ) = TEMP
00370          END IF
00371   130 CONTINUE
00372       X( 1, 1 ) = TMP( 1 )
00373       X( 2, 1 ) = TMP( 2 )
00374       X( 1, 2 ) = TMP( 3 )
00375       X( 2, 2 ) = TMP( 4 )
00376       XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
00377      $        ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
00378       RETURN
00379 *
00380 *     End of SLASY2
00381 *
00382       END