Return to dlasv2.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
   10: *     ..
   11: *
   12: *  Purpose
   13: *  =======
   14: *
   15: *  DLASV2 computes the singular value decomposition of a 2-by-2
   16: *  triangular matrix
   17: *     [  F   G  ]
   18: *     [  0   H  ].
   19: *  On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
   20: *  smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
   21: *  right singular vectors for abs(SSMAX), giving the decomposition
   22: *
   23: *     [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
   24: *     [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
   25: *
   26: *  Arguments
   27: *  =========
   28: *
   29: *  F       (input) DOUBLE PRECISION
   30: *          The (1,1) element of the 2-by-2 matrix.
   31: *
   32: *  G       (input) DOUBLE PRECISION
   33: *          The (1,2) element of the 2-by-2 matrix.
   34: *
   35: *  H       (input) DOUBLE PRECISION
   36: *          The (2,2) element of the 2-by-2 matrix.
   37: *
   38: *  SSMIN   (output) DOUBLE PRECISION
   39: *          abs(SSMIN) is the smaller singular value.
   40: *
   41: *  SSMAX   (output) DOUBLE PRECISION
   42: *          abs(SSMAX) is the larger singular value.
   43: *
   44: *  SNL     (output) DOUBLE PRECISION
   45: *  CSL     (output) DOUBLE PRECISION
   46: *          The vector (CSL, SNL) is a unit left singular vector for the
   47: *          singular value abs(SSMAX).
   48: *
   49: *  SNR     (output) DOUBLE PRECISION
   50: *  CSR     (output) DOUBLE PRECISION
   51: *          The vector (CSR, SNR) is a unit right singular vector for the
   52: *          singular value abs(SSMAX).
   53: *
   54: *  Further Details
   55: *  ===============
   56: *
   57: *  Any input parameter may be aliased with any output parameter.
   58: *
   59: *  Barring over/underflow and assuming a guard digit in subtraction, all
   60: *  output quantities are correct to within a few units in the last
   61: *  place (ulps).
   62: *
   63: *  In IEEE arithmetic, the code works correctly if one matrix element is
   64: *  infinite.
   65: *
   66: *  Overflow will not occur unless the largest singular value itself
   67: *  overflows or is within a few ulps of overflow. (On machines with
   68: *  partial overflow, like the Cray, overflow may occur if the largest
   69: *  singular value is within a factor of 2 of overflow.)
   70: *
   71: *  Underflow is harmless if underflow is gradual. Otherwise, results
   72: *  may correspond to a matrix modified by perturbations of size near
   73: *  the underflow threshold.
   74: *
   75: * =====================================================================
   76: *
   77: *     .. Parameters ..
   78:       DOUBLE PRECISION   ZERO
   79:       PARAMETER          ( ZERO = 0.0D0 )
   80:       DOUBLE PRECISION   HALF
   81:       PARAMETER          ( HALF = 0.5D0 )
   82:       DOUBLE PRECISION   ONE
   83:       PARAMETER          ( ONE = 1.0D0 )
   84:       DOUBLE PRECISION   TWO
   85:       PARAMETER          ( TWO = 2.0D0 )
   86:       DOUBLE PRECISION   FOUR
   87:       PARAMETER          ( FOUR = 4.0D0 )
   88: *     ..
   89: *     .. Local Scalars ..
   90:       LOGICAL            GASMAL, SWAP
   91:       INTEGER            PMAX
   92:       DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
   93:      $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
   94: *     ..
   95: *     .. Intrinsic Functions ..
   96:       INTRINSIC          ABS, SIGN, SQRT
   97: *     ..
   98: *     .. External Functions ..
   99:       DOUBLE PRECISION   DLAMCH
  100:       EXTERNAL           DLAMCH
  101: *     ..
  102: *     .. Executable Statements ..
  103: *
  104:       FT = F
  105:       FA = ABS( FT )
  106:       HT = H
  107:       HA = ABS( H )
  108: *
  109: *     PMAX points to the maximum absolute element of matrix
  110: *       PMAX = 1 if F largest in absolute values
  111: *       PMAX = 2 if G largest in absolute values
  112: *       PMAX = 3 if H largest in absolute values
  113: *
  114:       PMAX = 1
  115:       SWAP = ( HA.GT.FA )
  116:       IF( SWAP ) THEN
  117:          PMAX = 3
  118:          TEMP = FT
  119:          FT = HT
  120:          HT = TEMP
  121:          TEMP = FA
  122:          FA = HA
  123:          HA = TEMP
  124: *
  125: *        Now FA .ge. HA
  126: *
  127:       END IF
  128:       GT = G
  129:       GA = ABS( GT )
  130:       IF( GA.EQ.ZERO ) THEN
  131: *
  132: *        Diagonal matrix
  133: *
  134:          SSMIN = HA
  135:          SSMAX = FA
  136:          CLT = ONE
  137:          CRT = ONE
  138:          SLT = ZERO
  139:          SRT = ZERO
  140:       ELSE
  141:          GASMAL = .TRUE.
  142:          IF( GA.GT.FA ) THEN
  143:             PMAX = 2
  144:             IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
  145: *
  146: *              Case of very large GA
  147: *
  148:                GASMAL = .FALSE.
  149:                SSMAX = GA
  150:                IF( HA.GT.ONE ) THEN
  151:                   SSMIN = FA / ( GA / HA )
  152:                ELSE
  153:                   SSMIN = ( FA / GA )*HA
  154:                END IF
  155:                CLT = ONE
  156:                SLT = HT / GT
  157:                SRT = ONE
  158:                CRT = FT / GT
  159:             END IF
  160:          END IF
  161:          IF( GASMAL ) THEN
  162: *
  163: *           Normal case
  164: *
  165:             D = FA - HA
  166:             IF( D.EQ.FA ) THEN
  167: *
  168: *              Copes with infinite F or H
  169: *
  170:                L = ONE
  171:             ELSE
  172:                L = D / FA
  173:             END IF
  174: *
  175: *           Note that 0 .le. L .le. 1
  176: *
  177:             M = GT / FT
  178: *
  179: *           Note that abs(M) .le. 1/macheps
  180: *
  181:             T = TWO - L
  182: *
  183: *           Note that T .ge. 1
  184: *
  185:             MM = M*M
  186:             TT = T*T
  187:             S = SQRT( TT+MM )
  188: *
  189: *           Note that 1 .le. S .le. 1 + 1/macheps
  190: *
  191:             IF( L.EQ.ZERO ) THEN
  192:                R = ABS( M )
  193:             ELSE
  194:                R = SQRT( L*L+MM )
  195:             END IF
  196: *
  197: *           Note that 0 .le. R .le. 1 + 1/macheps
  198: *
  199:             A = HALF*( S+R )
  200: *
  201: *           Note that 1 .le. A .le. 1 + abs(M)
  202: *
  203:             SSMIN = HA / A
  204:             SSMAX = FA*A
  205:             IF( MM.EQ.ZERO ) THEN
  206: *
  207: *              Note that M is very tiny
  208: *
  209:                IF( L.EQ.ZERO ) THEN
  210:                   T = SIGN( TWO, FT )*SIGN( ONE, GT )
  211:                ELSE
  212:                   T = GT / SIGN( D, FT ) + M / T
  213:                END IF
  214:             ELSE
  215:                T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
  216:             END IF
  217:             L = SQRT( T*T+FOUR )
  218:             CRT = TWO / L
  219:             SRT = T / L
  220:             CLT = ( CRT+SRT*M ) / A
  221:             SLT = ( HT / FT )*SRT / A
  222:          END IF
  223:       END IF
  224:       IF( SWAP ) THEN
  225:          CSL = SRT
  226:          SNL = CRT
  227:          CSR = SLT
  228:          SNR = CLT
  229:       ELSE
  230:          CSL = CLT
  231:          SNL = SLT
  232:          CSR = CRT
  233:          SNR = SRT
  234:       END IF
  235: *
  236: *     Correct signs of SSMAX and SSMIN
  237: *
  238:       IF( PMAX.EQ.1 )
  239:      $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
  240:       IF( PMAX.EQ.2 )
  241:      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
  242:       IF( PMAX.EQ.3 )
  243:      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
  244:       SSMAX = SIGN( SSMAX, TSIGN )
  245:       SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
  246:       RETURN
  247: *
  248: *     End of DLASV2
  249: *
  250:       END