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