001:       SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
002:      $                   CTOT, W, S, 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, LDQ, N, N1
011:       REAL               RHO
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            CTOT( * ), INDX( * )
015:       REAL               D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
016:      $                   S( * ), W( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  SLAED3 finds the roots of the secular equation, as defined by the
023: *  values in D, W, and RHO, between 1 and K.  It makes the
024: *  appropriate calls to SLAED4 and then updates the eigenvectors by
025: *  multiplying the matrix of eigenvectors of the pair of eigensystems
026: *  being combined by the matrix of eigenvectors of the K-by-K system
027: *  which is solved here.
028: *
029: *  This code makes very mild assumptions about floating point
030: *  arithmetic. It will work on machines with a guard digit in
031: *  add/subtract, or on those binary machines without guard digits
032: *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
033: *  It could conceivably fail on hexadecimal or decimal machines
034: *  without guard digits, but we know of none.
035: *
036: *  Arguments
037: *  =========
038: *
039: *  K       (input) INTEGER
040: *          The number of terms in the rational function to be solved by
041: *          SLAED4.  K >= 0.
042: *
043: *  N       (input) INTEGER
044: *          The number of rows and columns in the Q matrix.
045: *          N >= K (deflation may result in N>K).
046: *
047: *  N1      (input) INTEGER
048: *          The location of the last eigenvalue in the leading submatrix.
049: *          min(1,N) <= N1 <= N/2.
050: *
051: *  D       (output) REAL array, dimension (N)
052: *          D(I) contains the updated eigenvalues for
053: *          1 <= I <= K.
054: *
055: *  Q       (output) REAL array, dimension (LDQ,N)
056: *          Initially the first K columns are used as workspace.
057: *          On output the columns 1 to K contain
058: *          the updated eigenvectors.
059: *
060: *  LDQ     (input) INTEGER
061: *          The leading dimension of the array Q.  LDQ >= max(1,N).
062: *
063: *  RHO     (input) REAL
064: *          The value of the parameter in the rank one update equation.
065: *          RHO >= 0 required.
066: *
067: *  DLAMDA  (input/output) REAL array, dimension (K)
068: *          The first K elements of this array contain the old roots
069: *          of the deflated updating problem.  These are the poles
070: *          of the secular equation. May be changed on output by
071: *          having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
072: *          Cray-2, or Cray C-90, as described above.
073: *
074: *  Q2      (input) REAL array, dimension (LDQ2, N)
075: *          The first K columns of this matrix contain the non-deflated
076: *          eigenvectors for the split problem.
077: *
078: *  INDX    (input) INTEGER array, dimension (N)
079: *          The permutation used to arrange the columns of the deflated
080: *          Q matrix into three groups (see SLAED2).
081: *          The rows of the eigenvectors found by SLAED4 must be likewise
082: *          permuted before the matrix multiply can take place.
083: *
084: *  CTOT    (input) INTEGER array, dimension (4)
085: *          A count of the total number of the various types of columns
086: *          in Q, as described in INDX.  The fourth column type is any
087: *          column which has been deflated.
088: *
089: *  W       (input/output) REAL array, dimension (K)
090: *          The first K elements of this array contain the components
091: *          of the deflation-adjusted updating vector. Destroyed on
092: *          output.
093: *
094: *  S       (workspace) REAL array, dimension (N1 + 1)*K
095: *          Will contain the eigenvectors of the repaired matrix which
096: *          will be multiplied by the previously accumulated eigenvectors
097: *          to update the system.
098: *
099: *  LDS     (input) INTEGER
100: *          The leading dimension of S.  LDS >= max(1,K).
101: *
102: *  INFO    (output) INTEGER
103: *          = 0:  successful exit.
104: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
105: *          > 0:  if INFO = 1, an eigenvalue did not converge
106: *
107: *  Further Details
108: *  ===============
109: *
110: *  Based on contributions by
111: *     Jeff Rutter, Computer Science Division, University of California
112: *     at Berkeley, USA
113: *  Modified by Francoise Tisseur, University of Tennessee.
114: *
115: *  =====================================================================
116: *
117: *     .. Parameters ..
118:       REAL               ONE, ZERO
119:       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
120: *     ..
121: *     .. Local Scalars ..
122:       INTEGER            I, II, IQ2, J, N12, N2, N23
123:       REAL               TEMP
124: *     ..
125: *     .. External Functions ..
126:       REAL               SLAMC3, SNRM2
127:       EXTERNAL           SLAMC3, SNRM2
128: *     ..
129: *     .. External Subroutines ..
130:       EXTERNAL           SCOPY, SGEMM, SLACPY, SLAED4, SLASET, XERBLA
131: *     ..
132: *     .. Intrinsic Functions ..
133:       INTRINSIC          MAX, SIGN, SQRT
134: *     ..
135: *     .. Executable Statements ..
136: *
137: *     Test the input parameters.
138: *
139:       INFO = 0
140: *
141:       IF( K.LT.0 ) THEN
142:          INFO = -1
143:       ELSE IF( N.LT.K ) THEN
144:          INFO = -2
145:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
146:          INFO = -6
147:       END IF
148:       IF( INFO.NE.0 ) THEN
149:          CALL XERBLA( 'SLAED3', -INFO )
150:          RETURN
151:       END IF
152: *
153: *     Quick return if possible
154: *
155:       IF( K.EQ.0 )
156:      $   RETURN
157: *
158: *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
159: *     be computed with high relative accuracy (barring over/underflow).
160: *     This is a problem on machines without a guard digit in
161: *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
162: *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
163: *     which on any of these machines zeros out the bottommost
164: *     bit of DLAMDA(I) if it is 1; this makes the subsequent
165: *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
166: *     occurs. On binary machines with a guard digit (almost all
167: *     machines) it does not change DLAMDA(I) at all. On hexadecimal
168: *     and decimal machines with a guard digit, it slightly
169: *     changes the bottommost bits of DLAMDA(I). It does not account
170: *     for hexadecimal or decimal machines without guard digits
171: *     (we know of none). We use a subroutine call to compute
172: *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
173: *     this code.
174: *
175:       DO 10 I = 1, K
176:          DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
177:    10 CONTINUE
178: *
179:       DO 20 J = 1, K
180:          CALL SLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
181: *
182: *        If the zero finder fails, the computation is terminated.
183: *
184:          IF( INFO.NE.0 )
185:      $      GO TO 120
186:    20 CONTINUE
187: *
188:       IF( K.EQ.1 )
189:      $   GO TO 110
190:       IF( K.EQ.2 ) THEN
191:          DO 30 J = 1, K
192:             W( 1 ) = Q( 1, J )
193:             W( 2 ) = Q( 2, J )
194:             II = INDX( 1 )
195:             Q( 1, J ) = W( II )
196:             II = INDX( 2 )
197:             Q( 2, J ) = W( II )
198:    30    CONTINUE
199:          GO TO 110
200:       END IF
201: *
202: *     Compute updated W.
203: *
204:       CALL SCOPY( K, W, 1, S, 1 )
205: *
206: *     Initialize W(I) = Q(I,I)
207: *
208:       CALL SCOPY( K, Q, LDQ+1, W, 1 )
209:       DO 60 J = 1, K
210:          DO 40 I = 1, J - 1
211:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
212:    40    CONTINUE
213:          DO 50 I = J + 1, K
214:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
215:    50    CONTINUE
216:    60 CONTINUE
217:       DO 70 I = 1, K
218:          W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
219:    70 CONTINUE
220: *
221: *     Compute eigenvectors of the modified rank-1 modification.
222: *
223:       DO 100 J = 1, K
224:          DO 80 I = 1, K
225:             S( I ) = W( I ) / Q( I, J )
226:    80    CONTINUE
227:          TEMP = SNRM2( K, S, 1 )
228:          DO 90 I = 1, K
229:             II = INDX( I )
230:             Q( I, J ) = S( II ) / TEMP
231:    90    CONTINUE
232:   100 CONTINUE
233: *
234: *     Compute the updated eigenvectors.
235: *
236:   110 CONTINUE
237: *
238:       N2 = N - N1
239:       N12 = CTOT( 1 ) + CTOT( 2 )
240:       N23 = CTOT( 2 ) + CTOT( 3 )
241: *
242:       CALL SLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
243:       IQ2 = N1*N12 + 1
244:       IF( N23.NE.0 ) THEN
245:          CALL SGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
246:      $               ZERO, Q( N1+1, 1 ), LDQ )
247:       ELSE
248:          CALL SLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
249:       END IF
250: *
251:       CALL SLACPY( 'A', N12, K, Q, LDQ, S, N12 )
252:       IF( N12.NE.0 ) THEN
253:          CALL SGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
254:      $               LDQ )
255:       ELSE
256:          CALL SLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
257:       END IF
258: *
259: *
260:   120 CONTINUE
261:       RETURN
262: *
263: *     End of SLAED3
264: *
265:       END
266: