Return to dlartgs.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dlartgs.f, revision 1.4

1.4     ! bertrand    1: *> \brief \b DLARTGS
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at 
        !             6: *            http://www.netlib.org/lapack/explore-html/ 
        !             7: *
        !             8: *> \htmlonly
        !             9: *> Download DLARTGS + dependencies 
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlartgs.f"> 
        !            11: *> [TGZ]</a> 
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlartgs.f"> 
        !            13: *> [ZIP]</a> 
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlartgs.f"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *       SUBROUTINE DLARTGS( X, Y, SIGMA, CS, SN )
        !            22: * 
        !            23: *       .. Scalar Arguments ..
        !            24: *       DOUBLE PRECISION        CS, SIGMA, SN, X, Y
        !            25: *       ..
        !            26: *  
        !            27: *
        !            28: *> \par Purpose:
        !            29: *  =============
        !            30: *>
        !            31: *> \verbatim
        !            32: *>
        !            33: *> DLARTGS generates a plane rotation designed to introduce a bulge in
        !            34: *> Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD
        !            35: *> problem. X and Y are the top-row entries, and SIGMA is the shift.
        !            36: *> The computed CS and SN define a plane rotation satisfying
        !            37: *>
        !            38: *>    [  CS  SN  ]  .  [ X^2 - SIGMA ]  =  [ R ],
        !            39: *>    [ -SN  CS  ]     [    X * Y    ]     [ 0 ]
        !            40: *>
        !            41: *> with R nonnegative.  If X^2 - SIGMA and X * Y are 0, then the
        !            42: *> rotation is by PI/2.
        !            43: *> \endverbatim
        !            44: *
        !            45: *  Arguments:
        !            46: *  ==========
        !            47: *
        !            48: *> \param[in] X
        !            49: *> \verbatim
        !            50: *>          X is DOUBLE PRECISION
        !            51: *>          The (1,1) entry of an upper bidiagonal matrix.
        !            52: *> \endverbatim
        !            53: *>
        !            54: *> \param[in] Y
        !            55: *> \verbatim
        !            56: *>          Y is DOUBLE PRECISION
        !            57: *>          The (1,2) entry of an upper bidiagonal matrix.
        !            58: *> \endverbatim
        !            59: *>
        !            60: *> \param[in] SIGMA
        !            61: *> \verbatim
        !            62: *>          SIGMA is DOUBLE PRECISION
        !            63: *>          The shift.
        !            64: *> \endverbatim
        !            65: *>
        !            66: *> \param[out] CS
        !            67: *> \verbatim
        !            68: *>          CS is DOUBLE PRECISION
        !            69: *>          The cosine of the rotation.
        !            70: *> \endverbatim
        !            71: *>
        !            72: *> \param[out] SN
        !            73: *> \verbatim
        !            74: *>          SN is DOUBLE PRECISION
        !            75: *>          The sine of the rotation.
        !            76: *> \endverbatim
        !            77: *
        !            78: *  Authors:
        !            79: *  ========
        !            80: *
        !            81: *> \author Univ. of Tennessee 
        !            82: *> \author Univ. of California Berkeley 
        !            83: *> \author Univ. of Colorado Denver 
        !            84: *> \author NAG Ltd. 
        !            85: *
        !            86: *> \date November 2011
1.1       bertrand   87: *
1.4     ! bertrand   88: *> \ingroup auxOTHERcomputational
1.1       bertrand   89: *
1.4     ! bertrand   90: *  =====================================================================
        !            91:       SUBROUTINE DLARTGS( X, Y, SIGMA, CS, SN )
1.1       bertrand   92: *
1.4     ! bertrand   93: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand   94: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
1.4     ! bertrand   95: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !            96: *     November 2011
1.1       bertrand   97: *
                     98: *     .. Scalar Arguments ..
                     99:       DOUBLE PRECISION        CS, SIGMA, SN, X, Y
                    100: *     ..
                    101: *
                    102: *  ===================================================================
                    103: *
                    104: *     .. Parameters ..
                    105:       DOUBLE PRECISION        NEGONE, ONE, ZERO
                    106:       PARAMETER          ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
                    107: *     ..
                    108: *     .. Local Scalars ..
                    109:       DOUBLE PRECISION        R, S, THRESH, W, Z
                    110: *     ..
                    111: *     .. External Functions ..
                    112:       DOUBLE PRECISION        DLAMCH
                    113:       EXTERNAL           DLAMCH
                    114: *     .. Executable Statements ..
                    115: *
                    116:       THRESH = DLAMCH('E')
                    117: *
1.3       bertrand  118: *     Compute the first column of B**T*B - SIGMA^2*I, up to a scale
1.1       bertrand  119: *     factor.
                    120: *
                    121:       IF( (SIGMA .EQ. ZERO .AND. ABS(X) .LT. THRESH) .OR.
                    122:      $          (ABS(X) .EQ. SIGMA .AND. Y .EQ. ZERO) ) THEN
                    123:          Z = ZERO
                    124:          W = ZERO
                    125:       ELSE IF( SIGMA .EQ. ZERO ) THEN
                    126:          IF( X .GE. ZERO ) THEN
                    127:             Z = X
                    128:             W = Y
                    129:          ELSE
                    130:             Z = -X
                    131:             W = -Y
                    132:          END IF
                    133:       ELSE IF( ABS(X) .LT. THRESH ) THEN
                    134:          Z = -SIGMA*SIGMA
                    135:          W = ZERO
                    136:       ELSE
                    137:          IF( X .GE. ZERO ) THEN
                    138:             S = ONE
                    139:          ELSE
                    140:             S = NEGONE
                    141:          END IF
                    142:          Z = S * (ABS(X)-SIGMA) * (S+SIGMA/X)
                    143:          W = S * Y
                    144:       END IF
                    145: *
                    146: *     Generate the rotation.
                    147: *     CALL DLARTGP( Z, W, CS, SN, R ) might seem more natural;
                    148: *     reordering the arguments ensures that if Z = 0 then the rotation
                    149: *     is by PI/2.
                    150: *
                    151:       CALL DLARTGP( W, Z, SN, CS, R )
                    152: *
                    153:       RETURN
                    154: *
                    155: *     End DLARTGS
                    156: *
                    157:       END
                    158: