SRC/dlaed3.f

   1 *> \brief \b DLAED3 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAED3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
  22 *                          CTOT, W, S, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, K, LDQ, N, N1
  26 *       DOUBLE PRECISION   RHO
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            CTOT( * ), INDX( * )
  30 *       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
  31 *      $                   S( * ), W( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DLAED3 finds the roots of the secular equation, as defined by the
  41 *> values in D, W, and RHO, between 1 and K.  It makes the
  42 *> appropriate calls to DLAED4 and then updates the eigenvectors by
  43 *> multiplying the matrix of eigenvectors of the pair of eigensystems
  44 *> being combined by the matrix of eigenvectors of the K-by-K system
  45 *> which is solved here.
  46 *>
  47 *> This code makes very mild assumptions about floating point
  48 *> arithmetic. It will work on machines with a guard digit in
  49 *> add/subtract, or on those binary machines without guard digits
  50 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
  51 *> It could conceivably fail on hexadecimal or decimal machines
  52 *> without guard digits, but we know of none.
  53 *> \endverbatim
  54 *
  55 *  Arguments:
  56 *  ==========
  57 *
  58 *> \param[in] K
  59 *> \verbatim
  60 *>          K is INTEGER
  61 *>          The number of terms in the rational function to be solved by
  62 *>          DLAED4.  K >= 0.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] N
  66 *> \verbatim
  67 *>          N is INTEGER
  68 *>          The number of rows and columns in the Q matrix.
  69 *>          N >= K (deflation may result in N>K).
  70 *> \endverbatim
  71 *>
  72 *> \param[in] N1
  73 *> \verbatim
  74 *>          N1 is INTEGER
  75 *>          The location of the last eigenvalue in the leading submatrix.
  76 *>          min(1,N) <= N1 <= N/2.
  77 *> \endverbatim
  78 *>
  79 *> \param[out] D
  80 *> \verbatim
  81 *>          D is DOUBLE PRECISION array, dimension (N)
  82 *>          D(I) contains the updated eigenvalues for
  83 *>          1 <= I <= K.
  84 *> \endverbatim
  85 *>
  86 *> \param[out] Q
  87 *> \verbatim
  88 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
  89 *>          Initially the first K columns are used as workspace.
  90 *>          On output the columns 1 to K contain
  91 *>          the updated eigenvectors.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] LDQ
  95 *> \verbatim
  96 *>          LDQ is INTEGER
  97 *>          The leading dimension of the array Q.  LDQ >= max(1,N).
  98 *> \endverbatim
  99 *>
 100 *> \param[in] RHO
 101 *> \verbatim
 102 *>          RHO is DOUBLE PRECISION
 103 *>          The value of the parameter in the rank one update equation.
 104 *>          RHO >= 0 required.
 105 *> \endverbatim
 106 *>
 107 *> \param[in,out] DLAMDA
 108 *> \verbatim
 109 *>          DLAMDA is DOUBLE PRECISION array, dimension (K)
 110 *>          The first K elements of this array contain the old roots
 111 *>          of the deflated updating problem.  These are the poles
 112 *>          of the secular equation. May be changed on output by
 113 *>          having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
 114 *>          Cray-2, or Cray C-90, as described above.
 115 *> \endverbatim
 116 *>
 117 *> \param[in] Q2
 118 *> \verbatim
 119 *>          Q2 is DOUBLE PRECISION array, dimension (LDQ2, N)
 120 *>          The first K columns of this matrix contain the non-deflated
 121 *>          eigenvectors for the split problem.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] INDX
 125 *> \verbatim
 126 *>          INDX is INTEGER array, dimension (N)
 127 *>          The permutation used to arrange the columns of the deflated
 128 *>          Q matrix into three groups (see DLAED2).
 129 *>          The rows of the eigenvectors found by DLAED4 must be likewise
 130 *>          permuted before the matrix multiply can take place.
 131 *> \endverbatim
 132 *>
 133 *> \param[in] CTOT
 134 *> \verbatim
 135 *>          CTOT is INTEGER array, dimension (4)
 136 *>          A count of the total number of the various types of columns
 137 *>          in Q, as described in INDX.  The fourth column type is any
 138 *>          column which has been deflated.
 139 *> \endverbatim
 140 *>
 141 *> \param[in,out] W
 142 *> \verbatim
 143 *>          W is DOUBLE PRECISION array, dimension (K)
 144 *>          The first K elements of this array contain the components
 145 *>          of the deflation-adjusted updating vector. Destroyed on
 146 *>          output.
 147 *> \endverbatim
 148 *>
 149 *> \param[out] S
 150 *> \verbatim
 151 *>          S is DOUBLE PRECISION array, dimension (N1 + 1)*K
 152 *>          Will contain the eigenvectors of the repaired matrix which
 153 *>          will be multiplied by the previously accumulated eigenvectors
 154 *>          to update the system.
 155 *> \endverbatim
 156 *>
 157 *> \param[out] INFO
 158 *> \verbatim
 159 *>          INFO is INTEGER
 160 *>          = 0:  successful exit.
 161 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 162 *>          > 0:  if INFO = 1, an eigenvalue did not converge
 163 *> \endverbatim
 164 *
 165 *  Authors:
 166 *  ========
 167 *
 168 *> \author Univ. of Tennessee
 169 *> \author Univ. of California Berkeley
 170 *> \author Univ. of Colorado Denver
 171 *> \author NAG Ltd.
 172 *
 173 *> \date September 2012
 174 *
 175 *> \ingroup auxOTHERcomputational
 176 *
 177 *> \par Contributors:
 178 *  ==================
 179 *>
 180 *> Jeff Rutter, Computer Science Division, University of California
 181 *> at Berkeley, USA \n
 182 *>  Modified by Francoise Tisseur, University of Tennessee
 183 *>
 184 *  =====================================================================
 185       SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
 186      $                   CTOT, W, S, INFO )
 187 *
 188 *  -- LAPACK computational routine (version 3.4.2) --
 189 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 190 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 191 *     September 2012
 192 *
 193 *     .. Scalar Arguments ..
 194       INTEGER            INFO, K, LDQ, N, N1
 195       DOUBLE PRECISION   RHO
 196 *     ..
 197 *     .. Array Arguments ..
 198       INTEGER            CTOT( * ), INDX( * )
 199       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
 200      $                   S( * ), W( * )
 201 *     ..
 202 *
 203 *  =====================================================================
 204 *
 205 *     .. Parameters ..
 206       DOUBLE PRECISION   ONE, ZERO
 207       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
 208 *     ..
 209 *     .. Local Scalars ..
 210       INTEGER            I, II, IQ2, J, N12, N2, N23
 211       DOUBLE PRECISION   TEMP
 212 *     ..
 213 *     .. External Functions ..
 214       DOUBLE PRECISION   DLAMC3, DNRM2
 215       EXTERNAL           DLAMC3, DNRM2
 216 *     ..
 217 *     .. External Subroutines ..
 218       EXTERNAL           DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA
 219 *     ..
 220 *     .. Intrinsic Functions ..
 221       INTRINSIC          MAX, SIGN, SQRT
 222 *     ..
 223 *     .. Executable Statements ..
 224 *
 225 *     Test the input parameters.
 226 *
 227       INFO = 0
 228 *
 229       IF( K.LT.0 ) THEN
 230          INFO = -1
 231       ELSE IF( N.LT.K ) THEN
 232          INFO = -2
 233       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
 234          INFO = -6
 235       END IF
 236       IF( INFO.NE.0 ) THEN
 237          CALL XERBLA( 'DLAED3', -INFO )
 238          RETURN
 239       END IF
 240 *
 241 *     Quick return if possible
 242 *
 243       IF( K.EQ.0 )
 244      $   RETURN
 245 *
 246 *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
 247 *     be computed with high relative accuracy (barring over/underflow).
 248 *     This is a problem on machines without a guard digit in
 249 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
 250 *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
 251 *     which on any of these machines zeros out the bottommost
 252 *     bit of DLAMDA(I) if it is 1; this makes the subsequent
 253 *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
 254 *     occurs. On binary machines with a guard digit (almost all
 255 *     machines) it does not change DLAMDA(I) at all. On hexadecimal
 256 *     and decimal machines with a guard digit, it slightly
 257 *     changes the bottommost bits of DLAMDA(I). It does not account
 258 *     for hexadecimal or decimal machines without guard digits
 259 *     (we know of none). We use a subroutine call to compute
 260 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
 261 *     this code.
 262 *
 263       DO 10 I = 1, K
 264          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
 265    10 CONTINUE
 266 *
 267       DO 20 J = 1, K
 268          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
 269 *
 270 *        If the zero finder fails, the computation is terminated.
 271 *
 272          IF( INFO.NE.0 )
 273      $      GO TO 120
 274    20 CONTINUE
 275 *
 276       IF( K.EQ.1 )
 277      $   GO TO 110
 278       IF( K.EQ.2 ) THEN
 279          DO 30 J = 1, K
 280             W( 1 ) = Q( 1, J )
 281             W( 2 ) = Q( 2, J )
 282             II = INDX( 1 )
 283             Q( 1, J ) = W( II )
 284             II = INDX( 2 )
 285             Q( 2, J ) = W( II )
 286    30    CONTINUE
 287          GO TO 110
 288       END IF
 289 *
 290 *     Compute updated W.
 291 *
 292       CALL DCOPY( K, W, 1, S, 1 )
 293 *
 294 *     Initialize W(I) = Q(I,I)
 295 *
 296       CALL DCOPY( K, Q, LDQ+1, W, 1 )
 297       DO 60 J = 1, K
 298          DO 40 I = 1, J - 1
 299             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
 300    40    CONTINUE
 301          DO 50 I = J + 1, K
 302             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
 303    50    CONTINUE
 304    60 CONTINUE
 305       DO 70 I = 1, K
 306          W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
 307    70 CONTINUE
 308 *
 309 *     Compute eigenvectors of the modified rank-1 modification.
 310 *
 311       DO 100 J = 1, K
 312          DO 80 I = 1, K
 313             S( I ) = W( I ) / Q( I, J )
 314    80    CONTINUE
 315          TEMP = DNRM2( K, S, 1 )
 316          DO 90 I = 1, K
 317             II = INDX( I )
 318             Q( I, J ) = S( II ) / TEMP
 319    90    CONTINUE
 320   100 CONTINUE
 321 *
 322 *     Compute the updated eigenvectors.
 323 *
 324   110 CONTINUE
 325 *
 326       N2 = N - N1
 327       N12 = CTOT( 1 ) + CTOT( 2 )
 328       N23 = CTOT( 2 ) + CTOT( 3 )
 329 *
 330       CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
 331       IQ2 = N1*N12 + 1
 332       IF( N23.NE.0 ) THEN
 333          CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
 334      $               ZERO, Q( N1+1, 1 ), LDQ )
 335       ELSE
 336          CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
 337       END IF
 338 *
 339       CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 )
 340       IF( N12.NE.0 ) THEN
 341          CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
 342      $               LDQ )
 343       ELSE
 344          CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
 345       END IF
 346 *
 347 *
 348   120 CONTINUE
 349       RETURN
 350 *
 351 *     End of DLAED3
 352 *
 353       END