dlaed9.f

*> \brief \b DLAED9 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAED9 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
*                          S, LDS, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
*       DOUBLE PRECISION   RHO
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
*      $                   W( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAED9 finds the roots of the secular equation, as defined by the
*> values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
*> appropriate calls to DLAED4 and then stores the new matrix of
*> eigenvectors for use in calculating the next level of Z vectors.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] K
*> \verbatim
*>          K is INTEGER
*>          The number of terms in the rational function to be solved by
*>          DLAED4.  K >= 0.
*> \endverbatim
*>
*> \param[in] KSTART
*> \verbatim
*>          KSTART is INTEGER
*> \endverbatim
*>
*> \param[in] KSTOP
*> \verbatim
*>          KSTOP is INTEGER
*>          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
*>          are to be computed.  1 <= KSTART <= KSTOP <= K.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of rows and columns in the Q matrix.
*>          N >= K (delation may result in N > K).
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          D(I) contains the updated eigenvalues
*>          for KSTART <= I <= KSTOP.
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.  LDQ >= max( 1, N ).
*> \endverbatim
*>
*> \param[in] RHO
*> \verbatim
*>          RHO is DOUBLE PRECISION
*>          The value of the parameter in the rank one update equation.
*>          RHO >= 0 required.
*> \endverbatim
*>
*> \param[in] DLAMDA
*> \verbatim
*>          DLAMDA is DOUBLE PRECISION array, dimension (K)
*>          The first K elements of this array contain the old roots
*>          of the deflated updating problem.  These are the poles
*>          of the secular equation.
*> \endverbatim
*>
*> \param[in] W
*> \verbatim
*>          W is DOUBLE PRECISION array, dimension (K)
*>          The first K elements of this array contain the components
*>          of the deflation-adjusted updating vector.
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION array, dimension (LDS, K)
*>          Will contain the eigenvectors of the repaired matrix which
*>          will be stored for subsequent Z vector calculation and
*>          multiplied by the previously accumulated eigenvectors
*>          to update the system.
*> \endverbatim
*>
*> \param[in] LDS
*> \verbatim
*>          LDS is INTEGER
*>          The leading dimension of S.  LDS >= max( 1, K ).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*>          > 0:  if INFO = 1, an eigenvalue did not converge
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup laed9
*
*> \par Contributors:
*  ==================
*>
*> Jeff Rutter, Computer Science Division, University of California
*> at Berkeley, USA
*
*  =====================================================================
      SUBROUTINE dlaed9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
     $                   S, LDS, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
      DOUBLE PRECISION   RHO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( ldq, * ), S( lds, * ),
     $                   w( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, J
      DOUBLE PRECISION   TEMP
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMC3, DNRM2
      EXTERNAL           dlamc3, dnrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlaed4, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, sign, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      IF( k.LT.0 ) THEN
         info = -1
      ELSE IF( kstart.LT.1 .OR. kstart.GT.max( 1, k ) ) THEN
         info = -2
      ELSE IF( max( 1, kstop ).LT.kstart .OR. kstop.GT.max( 1, k ) )
     $          THEN
         info = -3
      ELSE IF( n.LT.k ) THEN
         info = -4
      ELSE IF( ldq.LT.max( 1, k ) ) THEN
         info = -7
      ELSE IF( lds.LT.max( 1, k ) ) THEN
         info = -12
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DLAED9', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( k.EQ.0 )
     $   RETURN
*
*     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
*     be computed with high relative accuracy (barring over/underflow).
*     This is a problem on machines without a guard digit in
*     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
*     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
*     which on any of these machines zeros out the bottommost
*     bit of DLAMDA(I) if it is 1; this makes the subsequent
*     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
*     occurs. On binary machines with a guard digit (almost all
*     machines) it does not change DLAMDA(I) at all. On hexadecimal
*     and decimal machines with a guard digit, it slightly
*     changes the bottommost bits of DLAMDA(I). It does not account
*     for hexadecimal or decimal machines without guard digits
*     (we know of none). We use a subroutine call to compute
*     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
*     this code.
*
      DO 10 i = 1, n
         dlamda( i ) = dlamc3( dlamda( i ), dlamda( i ) ) - dlamda( i )
   10 CONTINUE
*
      DO 20 j = kstart, kstop
         CALL dlaed4( k, j, dlamda, w, q( 1, j ), rho, d( j ), info )
*
*        If the zero finder fails, the computation is terminated.
*
         IF( info.NE.0 )
     $      GO TO 120
   20 CONTINUE
*
      IF( k.EQ.1 .OR. k.EQ.2 ) THEN
         DO 40 i = 1, k
            DO 30 j = 1, k
               s( j, i ) = q( j, i )
   30       CONTINUE
   40    CONTINUE
         GO TO 120
      END IF
*
*     Compute updated W.
*
      CALL dcopy( k, w, 1, s, 1 )
*
*     Initialize W(I) = Q(I,I)
*
      CALL dcopy( k, q, ldq+1, w, 1 )
      DO 70 j = 1, k
         DO 50 i = 1, j - 1
            w( i ) = w( i )*( q( i, j ) / ( dlamda( i )-dlamda( j ) ) )
   50    CONTINUE
         DO 60 i = j + 1, k
            w( i ) = w( i )*( q( i, j ) / ( dlamda( i )-dlamda( j ) ) )
   60    CONTINUE
   70 CONTINUE
      DO 80 i = 1, k
         w( i ) = sign( sqrt( -w( i ) ), s( i, 1 ) )
   80 CONTINUE
*
*     Compute eigenvectors of the modified rank-1 modification.
*
      DO 110 j = 1, k
         DO 90 i = 1, k
            q( i, j ) = w( i ) / q( i, j )
   90    CONTINUE
         temp = dnrm2( k, q( 1, j ), 1 )
         DO 100 i = 1, k
            s( i, j ) = q( i, j ) / temp
  100    CONTINUE
  110 CONTINUE
*
  120 CONTINUE
      RETURN
*
*     End of DLAED9
*
      END