df/de6/dlartgs_8f_source.html

00001       SUBROUTINE DLARTGS( X, Y, SIGMA, CS, SN )
00002       IMPLICIT NONE
00003 *
00004 *  -- LAPACK routine (version 3.3.0) --
00005 *
00006 *  -- Contributed by Brian Sutton of the Randolph-Macon College --
00007 *  -- November 2010
00008 *
00009 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00010 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00011 *
00012 *     .. Scalar Arguments ..
00013       DOUBLE PRECISION        CS, SIGMA, SN, X, Y
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DLARTGS generates a plane rotation designed to introduce a bulge in
00020 *  Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD
00021 *  problem. X and Y are the top-row entries, and SIGMA is the shift.
00022 *  The computed CS and SN define a plane rotation satisfying
00023 *
00024 *     [  CS  SN  ]  .  [ X^2 - SIGMA ]  =  [ R ],
00025 *     [ -SN  CS  ]     [    X * Y    ]     [ 0 ]
00026 *
00027 *  with R nonnegative.  If X^2 - SIGMA and X * Y are 0, then the
00028 *  rotation is by PI/2.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  X       (input) DOUBLE PRECISION
00034 *          The (1,1) entry of an upper bidiagonal matrix.
00035 *
00036 *  Y       (input) DOUBLE PRECISION
00037 *          The (1,2) entry of an upper bidiagonal matrix.
00038 *
00039 *  SIGMA   (input) DOUBLE PRECISION
00040 *          The shift.
00041 *
00042 *  CS      (output) DOUBLE PRECISION
00043 *          The cosine of the rotation.
00044 *
00045 *  SN      (output) DOUBLE PRECISION
00046 *          The sine of the rotation.
00047 *
00048 *  ===================================================================
00049 *
00050 *     .. Parameters ..
00051       DOUBLE PRECISION        NEGONE, ONE, ZERO
00052       PARAMETER          ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
00053 *     ..
00054 *     .. Local Scalars ..
00055       DOUBLE PRECISION        R, S, THRESH, W, Z
00056 *     ..
00057 *     .. External Functions ..
00058       DOUBLE PRECISION        DLAMCH
00059       EXTERNAL           DLAMCH
00060 *     .. Executable Statements ..
00061 *
00062       THRESH = DLAMCH('E')
00063 *
00064 *     Compute the first column of B'*B - SIGMA^2*I, up to a scale
00065 *     factor.
00066 *
00067       IF( (SIGMA .EQ. ZERO .AND. ABS(X) .LT. THRESH) .OR.
00068      $          (ABS(X) .EQ. SIGMA .AND. Y .EQ. ZERO) ) THEN
00069          Z = ZERO
00070          W = ZERO
00071       ELSE IF( SIGMA .EQ. ZERO ) THEN
00072          IF( X .GE. ZERO ) THEN
00073             Z = X
00074             W = Y
00075          ELSE
00076             Z = -X
00077             W = -Y
00078          END IF
00079       ELSE IF( ABS(X) .LT. THRESH ) THEN
00080          Z = -SIGMA*SIGMA
00081          W = ZERO
00082       ELSE
00083          IF( X .GE. ZERO ) THEN
00084             S = ONE
00085          ELSE
00086             S = NEGONE
00087          END IF
00088          Z = S * (ABS(X)-SIGMA) * (S+SIGMA/X)
00089          W = S * Y
00090       END IF
00091 *
00092 *     Generate the rotation.
00093 *     CALL DLARTGP( Z, W, CS, SN, R ) might seem more natural;
00094 *     reordering the arguments ensures that if Z = 0 then the rotation
00095 *     is by PI/2.
00096 *
00097       CALL DLARTGP( W, Z, SN, CS, R )
00098 *
00099       RETURN
00100 *
00101 *     End DLARTGS
00102 *
00103       END
00104