#include "ESMF_LapackBlas.inc" !> \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 !> <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. ! !> \date December 2016 ! !> \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 (version 3.7.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! December 2016 ! ! .. 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