SRC/slasd5.f

   1 *> \brief \b SLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. Used by sbdsdc.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLASD5 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd5.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd5.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd5.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            I
  25 *       REAL               DSIGMA, RHO
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> This subroutine computes the square root of the I-th eigenvalue
  38 *> of a positive symmetric rank-one modification of a 2-by-2 diagonal
  39 *> matrix
  40 *>
  41 *>            diag( D ) * diag( D ) +  RHO * Z * transpose(Z) .
  42 *>
  43 *> The diagonal entries in the array D are assumed to satisfy
  44 *>
  45 *>            0 <= D(i) < D(j)  for  i < j .
  46 *>
  47 *> We also assume RHO > 0 and that the Euclidean norm of the vector
  48 *> Z is one.
  49 *> \endverbatim
  50 *
  51 *  Arguments:
  52 *  ==========
  53 *
  54 *> \param[in] I
  55 *> \verbatim
  56 *>          I is INTEGER
  57 *>         The index of the eigenvalue to be computed.  I = 1 or I = 2.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] D
  61 *> \verbatim
  62 *>          D is REAL array, dimension (2)
  63 *>         The original eigenvalues.  We assume 0 <= D(1) < D(2).
  64 *> \endverbatim
  65 *>
  66 *> \param[in] Z
  67 *> \verbatim
  68 *>          Z is REAL array, dimension (2)
  69 *>         The components of the updating vector.
  70 *> \endverbatim
  71 *>
  72 *> \param[out] DELTA
  73 *> \verbatim
  74 *>          DELTA is REAL array, dimension (2)
  75 *>         Contains (D(j) - sigma_I) in its  j-th component.
  76 *>         The vector DELTA contains the information necessary
  77 *>         to construct the eigenvectors.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] RHO
  81 *> \verbatim
  82 *>          RHO is REAL
  83 *>         The scalar in the symmetric updating formula.
  84 *> \endverbatim
  85 *>
  86 *> \param[out] DSIGMA
  87 *> \verbatim
  88 *>          DSIGMA is REAL
  89 *>         The computed sigma_I, the I-th updated eigenvalue.
  90 *> \endverbatim
  91 *>
  92 *> \param[out] WORK
  93 *> \verbatim
  94 *>          WORK is REAL array, dimension (2)
  95 *>         WORK contains (D(j) + sigma_I) in its  j-th component.
  96 *> \endverbatim
  97 *
  98 *  Authors:
  99 *  ========
 100 *
 101 *> \author Univ. of Tennessee
 102 *> \author Univ. of California Berkeley
 103 *> \author Univ. of Colorado Denver
 104 *> \author NAG Ltd.
 105 *
 106 *> \date September 2012
 107 *
 108 *> \ingroup auxOTHERauxiliary
 109 *
 110 *> \par Contributors:
 111 *  ==================
 112 *>
 113 *>     Ren-Cang Li, Computer Science Division, University of California
 114 *>     at Berkeley, USA
 115 *>
 116 *  =====================================================================
 117       SUBROUTINE SLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
 118 *
 119 *  -- LAPACK auxiliary routine (version 3.4.2) --
 120 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 121 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 122 *     September 2012
 123 *
 124 *     .. Scalar Arguments ..
 125       INTEGER            I
 126       REAL               DSIGMA, RHO
 127 *     ..
 128 *     .. Array Arguments ..
 129       REAL               D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
 130 *     ..
 131 *
 132 *  =====================================================================
 133 *
 134 *     .. Parameters ..
 135       REAL               ZERO, ONE, TWO, THREE, FOUR
 136       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
 137      $                   THREE = 3.0E+0, FOUR = 4.0E+0 )
 138 *     ..
 139 *     .. Local Scalars ..
 140       REAL               B, C, DEL, DELSQ, TAU, W
 141 *     ..
 142 *     .. Intrinsic Functions ..
 143       INTRINSIC          ABS, SQRT
 144 *     ..
 145 *     .. Executable Statements ..
 146 *
 147       DEL = D( 2 ) - D( 1 )
 148       DELSQ = DEL*( D( 2 )+D( 1 ) )
 149       IF( I.EQ.1 ) THEN
 150          W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )-
 151      $       Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL
 152          IF( W.GT.ZERO ) THEN
 153             B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
 154             C = RHO*Z( 1 )*Z( 1 )*DELSQ
 155 *
 156 *           B > ZERO, always
 157 *
 158 *           The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 )
 159 *
 160             TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) )
 161 *
 162 *           The following TAU is DSIGMA - D( 1 )
 163 *
 164             TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) )
 165             DSIGMA = D( 1 ) + TAU
 166             DELTA( 1 ) = -TAU
 167             DELTA( 2 ) = DEL - TAU
 168             WORK( 1 ) = TWO*D( 1 ) + TAU
 169             WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 )
 170 *           DELTA( 1 ) = -Z( 1 ) / TAU
 171 *           DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )
 172          ELSE
 173             B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
 174             C = RHO*Z( 2 )*Z( 2 )*DELSQ
 175 *
 176 *           The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
 177 *
 178             IF( B.GT.ZERO ) THEN
 179                TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) )
 180             ELSE
 181                TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO
 182             END IF
 183 *
 184 *           The following TAU is DSIGMA - D( 2 )
 185 *
 186             TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) )
 187             DSIGMA = D( 2 ) + TAU
 188             DELTA( 1 ) = -( DEL+TAU )
 189             DELTA( 2 ) = -TAU
 190             WORK( 1 ) = D( 1 ) + TAU + D( 2 )
 191             WORK( 2 ) = TWO*D( 2 ) + TAU
 192 *           DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
 193 *           DELTA( 2 ) = -Z( 2 ) / TAU
 194          END IF
 195 *        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
 196 *        DELTA( 1 ) = DELTA( 1 ) / TEMP
 197 *        DELTA( 2 ) = DELTA( 2 ) / TEMP
 198       ELSE
 199 *
 200 *        Now I=2
 201 *
 202          B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
 203          C = RHO*Z( 2 )*Z( 2 )*DELSQ
 204 *
 205 *        The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
 206 *
 207          IF( B.GT.ZERO ) THEN
 208             TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO
 209          ELSE
 210             TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) )
 211          END IF
 212 *
 213 *        The following TAU is DSIGMA - D( 2 )
 214 *
 215          TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) )
 216          DSIGMA = D( 2 ) + TAU
 217          DELTA( 1 ) = -( DEL+TAU )
 218          DELTA( 2 ) = -TAU
 219          WORK( 1 ) = D( 1 ) + TAU + D( 2 )
 220          WORK( 2 ) = TWO*D( 2 ) + TAU
 221 *        DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
 222 *        DELTA( 2 ) = -Z( 2 ) / TAU
 223 *        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
 224 *        DELTA( 1 ) = DELTA( 1 ) / TEMP
 225 *        DELTA( 2 ) = DELTA( 2 ) / TEMP
 226       END IF
 227       RETURN
 228 *
 229 *     End of SLASD5
 230 *
 231       END