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