SRC/dlaed9.f

   1 *> \brief \b DLAED9 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAED9 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
  22 *                          S, LDS, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
  26 *       DOUBLE PRECISION   RHO
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
  30 *      $                   W( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DLAED9 finds the roots of the secular equation, as defined by the
  40 *> values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
  41 *> appropriate calls to DLAED4 and then stores the new matrix of
  42 *> eigenvectors for use in calculating the next level of Z vectors.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] K
  49 *> \verbatim
  50 *>          K is INTEGER
  51 *>          The number of terms in the rational function to be solved by
  52 *>          DLAED4.  K >= 0.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] KSTART
  56 *> \verbatim
  57 *>          KSTART is INTEGER
  58 *> \endverbatim
  59 *>
  60 *> \param[in] KSTOP
  61 *> \verbatim
  62 *>          KSTOP is INTEGER
  63 *>          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
  64 *>          are to be computed.  1 <= KSTART <= KSTOP <= K.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] N
  68 *> \verbatim
  69 *>          N is INTEGER
  70 *>          The number of rows and columns in the Q matrix.
  71 *>          N >= K (delation may result in N > K).
  72 *> \endverbatim
  73 *>
  74 *> \param[out] D
  75 *> \verbatim
  76 *>          D is DOUBLE PRECISION array, dimension (N)
  77 *>          D(I) contains the updated eigenvalues
  78 *>          for KSTART <= I <= KSTOP.
  79 *> \endverbatim
  80 *>
  81 *> \param[out] Q
  82 *> \verbatim
  83 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
  84 *> \endverbatim
  85 *>
  86 *> \param[in] LDQ
  87 *> \verbatim
  88 *>          LDQ is INTEGER
  89 *>          The leading dimension of the array Q.  LDQ >= max( 1, N ).
  90 *> \endverbatim
  91 *>
  92 *> \param[in] RHO
  93 *> \verbatim
  94 *>          RHO is DOUBLE PRECISION
  95 *>          The value of the parameter in the rank one update equation.
  96 *>          RHO >= 0 required.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] DLAMDA
 100 *> \verbatim
 101 *>          DLAMDA is DOUBLE PRECISION array, dimension (K)
 102 *>          The first K elements of this array contain the old roots
 103 *>          of the deflated updating problem.  These are the poles
 104 *>          of the secular equation.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] W
 108 *> \verbatim
 109 *>          W is DOUBLE PRECISION array, dimension (K)
 110 *>          The first K elements of this array contain the components
 111 *>          of the deflation-adjusted updating vector.
 112 *> \endverbatim
 113 *>
 114 *> \param[out] S
 115 *> \verbatim
 116 *>          S is DOUBLE PRECISION array, dimension (LDS, K)
 117 *>          Will contain the eigenvectors of the repaired matrix which
 118 *>          will be stored for subsequent Z vector calculation and
 119 *>          multiplied by the previously accumulated eigenvectors
 120 *>          to update the system.
 121 *> \endverbatim
 122 *>
 123 *> \param[in] LDS
 124 *> \verbatim
 125 *>          LDS is INTEGER
 126 *>          The leading dimension of S.  LDS >= max( 1, K ).
 127 *> \endverbatim
 128 *>
 129 *> \param[out] INFO
 130 *> \verbatim
 131 *>          INFO is INTEGER
 132 *>          = 0:  successful exit.
 133 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 134 *>          > 0:  if INFO = 1, an eigenvalue did not converge
 135 *> \endverbatim
 136 *
 137 *  Authors:
 138 *  ========
 139 *
 140 *> \author Univ. of Tennessee
 141 *> \author Univ. of California Berkeley
 142 *> \author Univ. of Colorado Denver
 143 *> \author NAG Ltd.
 144 *
 145 *> \date September 2012
 146 *
 147 *> \ingroup auxOTHERcomputational
 148 *
 149 *> \par Contributors:
 150 *  ==================
 151 *>
 152 *> Jeff Rutter, Computer Science Division, University of California
 153 *> at Berkeley, USA
 154 *
 155 *  =====================================================================
 156       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
 157      $                   S, LDS, INFO )
 158 *
 159 *  -- LAPACK computational routine (version 3.4.2) --
 160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 162 *     September 2012
 163 *
 164 *     .. Scalar Arguments ..
 165       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
 166       DOUBLE PRECISION   RHO
 167 *     ..
 168 *     .. Array Arguments ..
 169       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
 170      $                   W( * )
 171 *     ..
 172 *
 173 *  =====================================================================
 174 *
 175 *     .. Local Scalars ..
 176       INTEGER            I, J
 177       DOUBLE PRECISION   TEMP
 178 *     ..
 179 *     .. External Functions ..
 180       DOUBLE PRECISION   DLAMC3, DNRM2
 181       EXTERNAL           DLAMC3, DNRM2
 182 *     ..
 183 *     .. External Subroutines ..
 184       EXTERNAL           DCOPY, DLAED4, XERBLA
 185 *     ..
 186 *     .. Intrinsic Functions ..
 187       INTRINSIC          MAX, SIGN, SQRT
 188 *     ..
 189 *     .. Executable Statements ..
 190 *
 191 *     Test the input parameters.
 192 *
 193       INFO = 0
 194 *
 195       IF( K.LT.0 ) THEN
 196          INFO = -1
 197       ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
 198          INFO = -2
 199       ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
 200      $          THEN
 201          INFO = -3
 202       ELSE IF( N.LT.K ) THEN
 203          INFO = -4
 204       ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
 205          INFO = -7
 206       ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
 207          INFO = -12
 208       END IF
 209       IF( INFO.NE.0 ) THEN
 210          CALL XERBLA( 'DLAED9', -INFO )
 211          RETURN
 212       END IF
 213 *
 214 *     Quick return if possible
 215 *
 216       IF( K.EQ.0 )
 217      $   RETURN
 218 *
 219 *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
 220 *     be computed with high relative accuracy (barring over/underflow).
 221 *     This is a problem on machines without a guard digit in
 222 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
 223 *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
 224 *     which on any of these machines zeros out the bottommost
 225 *     bit of DLAMDA(I) if it is 1; this makes the subsequent
 226 *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
 227 *     occurs. On binary machines with a guard digit (almost all
 228 *     machines) it does not change DLAMDA(I) at all. On hexadecimal
 229 *     and decimal machines with a guard digit, it slightly
 230 *     changes the bottommost bits of DLAMDA(I). It does not account
 231 *     for hexadecimal or decimal machines without guard digits
 232 *     (we know of none). We use a subroutine call to compute
 233 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
 234 *     this code.
 235 *
 236       DO 10 I = 1, N
 237          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
 238    10 CONTINUE
 239 *
 240       DO 20 J = KSTART, KSTOP
 241          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
 242 *
 243 *        If the zero finder fails, the computation is terminated.
 244 *
 245          IF( INFO.NE.0 )
 246      $      GO TO 120
 247    20 CONTINUE
 248 *
 249       IF( K.EQ.1 .OR. K.EQ.2 ) THEN
 250          DO 40 I = 1, K
 251             DO 30 J = 1, K
 252                S( J, I ) = Q( J, I )
 253    30       CONTINUE
 254    40    CONTINUE
 255          GO TO 120
 256       END IF
 257 *
 258 *     Compute updated W.
 259 *
 260       CALL DCOPY( K, W, 1, S, 1 )
 261 *
 262 *     Initialize W(I) = Q(I,I)
 263 *
 264       CALL DCOPY( K, Q, LDQ+1, W, 1 )
 265       DO 70 J = 1, K
 266          DO 50 I = 1, J - 1
 267             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
 268    50    CONTINUE
 269          DO 60 I = J + 1, K
 270             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
 271    60    CONTINUE
 272    70 CONTINUE
 273       DO 80 I = 1, K
 274          W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
 275    80 CONTINUE
 276 *
 277 *     Compute eigenvectors of the modified rank-1 modification.
 278 *
 279       DO 110 J = 1, K
 280          DO 90 I = 1, K
 281             Q( I, J ) = W( I ) / Q( I, J )
 282    90    CONTINUE
 283          TEMP = DNRM2( K, Q( 1, J ), 1 )
 284          DO 100 I = 1, K
 285             S( I, J ) = Q( I, J ) / TEMP
 286   100    CONTINUE
 287   110 CONTINUE
 288 *
 289   120 CONTINUE
 290       RETURN
 291 *
 292 *     End of DLAED9
 293 *
 294       END