SRC/dlartgs.f

   1 *> \brief \b DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem.
   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 September 2012
  87 *
  88 *> \ingroup auxOTHERcomputational
  89 *
  90 *  =====================================================================
  91       SUBROUTINE DLARTGS( X, Y, SIGMA, CS, SN )
  92 *
  93 *  -- LAPACK computational routine (version 3.4.2) --
  94 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  95 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  96 *     September 2012
  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 *
 118 *     Compute the first column of B**T*B - SIGMA^2*I, up to a scale
 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