001:       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
002:      $                   S, LDS, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
011:       DOUBLE PRECISION   RHO
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
015:      $                   W( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DLAED9 finds the roots of the secular equation, as defined by the
022: *  values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
023: *  appropriate calls to DLAED4 and then stores the new matrix of
024: *  eigenvectors for use in calculating the next level of Z vectors.
025: *
026: *  Arguments
027: *  =========
028: *
029: *  K       (input) INTEGER
030: *          The number of terms in the rational function to be solved by
031: *          DLAED4.  K >= 0.
032: *
033: *  KSTART  (input) INTEGER
034: *  KSTOP   (input) INTEGER
035: *          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
036: *          are to be computed.  1 <= KSTART <= KSTOP <= K.
037: *
038: *  N       (input) INTEGER
039: *          The number of rows and columns in the Q matrix.
040: *          N >= K (delation may result in N > K).
041: *
042: *  D       (output) DOUBLE PRECISION array, dimension (N)
043: *          D(I) contains the updated eigenvalues
044: *          for KSTART <= I <= KSTOP.
045: *
046: *  Q       (workspace) DOUBLE PRECISION array, dimension (LDQ,N)
047: *
048: *  LDQ     (input) INTEGER
049: *          The leading dimension of the array Q.  LDQ >= max( 1, N ).
050: *
051: *  RHO     (input) DOUBLE PRECISION
052: *          The value of the parameter in the rank one update equation.
053: *          RHO >= 0 required.
054: *
055: *  DLAMDA  (input) DOUBLE PRECISION array, dimension (K)
056: *          The first K elements of this array contain the old roots
057: *          of the deflated updating problem.  These are the poles
058: *          of the secular equation.
059: *
060: *  W       (input) DOUBLE PRECISION array, dimension (K)
061: *          The first K elements of this array contain the components
062: *          of the deflation-adjusted updating vector.
063: *
064: *  S       (output) DOUBLE PRECISION array, dimension (LDS, K)
065: *          Will contain the eigenvectors of the repaired matrix which
066: *          will be stored for subsequent Z vector calculation and
067: *          multiplied by the previously accumulated eigenvectors
068: *          to update the system.
069: *
070: *  LDS     (input) INTEGER
071: *          The leading dimension of S.  LDS >= max( 1, K ).
072: *
073: *  INFO    (output) INTEGER
074: *          = 0:  successful exit.
075: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
076: *          > 0:  if INFO = 1, an eigenvalue did not converge
077: *
078: *  Further Details
079: *  ===============
080: *
081: *  Based on contributions by
082: *     Jeff Rutter, Computer Science Division, University of California
083: *     at Berkeley, USA
084: *
085: *  =====================================================================
086: *
087: *     .. Local Scalars ..
088:       INTEGER            I, J
089:       DOUBLE PRECISION   TEMP
090: *     ..
091: *     .. External Functions ..
092:       DOUBLE PRECISION   DLAMC3, DNRM2
093:       EXTERNAL           DLAMC3, DNRM2
094: *     ..
095: *     .. External Subroutines ..
096:       EXTERNAL           DCOPY, DLAED4, XERBLA
097: *     ..
098: *     .. Intrinsic Functions ..
099:       INTRINSIC          MAX, SIGN, SQRT
100: *     ..
101: *     .. Executable Statements ..
102: *
103: *     Test the input parameters.
104: *
105:       INFO = 0
106: *
107:       IF( K.LT.0 ) THEN
108:          INFO = -1
109:       ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
110:          INFO = -2
111:       ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
112:      $          THEN
113:          INFO = -3
114:       ELSE IF( N.LT.K ) THEN
115:          INFO = -4
116:       ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
117:          INFO = -7
118:       ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
119:          INFO = -12
120:       END IF
121:       IF( INFO.NE.0 ) THEN
122:          CALL XERBLA( 'DLAED9', -INFO )
123:          RETURN
124:       END IF
125: *
126: *     Quick return if possible
127: *
128:       IF( K.EQ.0 )
129:      $   RETURN
130: *
131: *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
132: *     be computed with high relative accuracy (barring over/underflow).
133: *     This is a problem on machines without a guard digit in
134: *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
135: *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
136: *     which on any of these machines zeros out the bottommost
137: *     bit of DLAMDA(I) if it is 1; this makes the subsequent
138: *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
139: *     occurs. On binary machines with a guard digit (almost all
140: *     machines) it does not change DLAMDA(I) at all. On hexadecimal
141: *     and decimal machines with a guard digit, it slightly
142: *     changes the bottommost bits of DLAMDA(I). It does not account
143: *     for hexadecimal or decimal machines without guard digits
144: *     (we know of none). We use a subroutine call to compute
145: *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
146: *     this code.
147: *
148:       DO 10 I = 1, N
149:          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
150:    10 CONTINUE
151: *
152:       DO 20 J = KSTART, KSTOP
153:          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
154: *
155: *        If the zero finder fails, the computation is terminated.
156: *
157:          IF( INFO.NE.0 )
158:      $      GO TO 120
159:    20 CONTINUE
160: *
161:       IF( K.EQ.1 .OR. K.EQ.2 ) THEN
162:          DO 40 I = 1, K
163:             DO 30 J = 1, K
164:                S( J, I ) = Q( J, I )
165:    30       CONTINUE
166:    40    CONTINUE
167:          GO TO 120
168:       END IF
169: *
170: *     Compute updated W.
171: *
172:       CALL DCOPY( K, W, 1, S, 1 )
173: *
174: *     Initialize W(I) = Q(I,I)
175: *
176:       CALL DCOPY( K, Q, LDQ+1, W, 1 )
177:       DO 70 J = 1, K
178:          DO 50 I = 1, J - 1
179:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
180:    50    CONTINUE
181:          DO 60 I = J + 1, K
182:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
183:    60    CONTINUE
184:    70 CONTINUE
185:       DO 80 I = 1, K
186:          W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
187:    80 CONTINUE
188: *
189: *     Compute eigenvectors of the modified rank-1 modification.
190: *
191:       DO 110 J = 1, K
192:          DO 90 I = 1, K
193:             Q( I, J ) = W( I ) / Q( I, J )
194:    90    CONTINUE
195:          TEMP = DNRM2( K, Q( 1, J ), 1 )
196:          DO 100 I = 1, K
197:             S( I, J ) = Q( I, J ) / TEMP
198:   100    CONTINUE
199:   110 CONTINUE
200: *
201:   120 CONTINUE
202:       RETURN
203: *
204: *     End of DLAED9
205: *
206:       END
207: