*> \brief \b DLAED9 used by sstedc. 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 *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \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 auxOTHERcomputational * *> \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