001:       SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
002:      $                   DSIGMA, WORK, INFO )
003: *
004: *  -- LAPACK auxiliary 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: *     October 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            ICOMPQ, INFO, K, LDDIFR
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               D( * ), DIFL( * ), DIFR( LDDIFR, * ),
014:      $                   DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
015:      $                   Z( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLASD8 finds the square roots of the roots of the secular equation,
022: *  as defined by the values in DSIGMA and Z. It makes the appropriate
023: *  calls to SLASD4, and stores, for each  element in D, the distance
024: *  to its two nearest poles (elements in DSIGMA). It also updates
025: *  the arrays VF and VL, the first and last components of all the
026: *  right singular vectors of the original bidiagonal matrix.
027: *
028: *  SLASD8 is called from SLASD6.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  ICOMPQ  (input) INTEGER
034: *          Specifies whether singular vectors are to be computed in
035: *          factored form in the calling routine:
036: *          = 0: Compute singular values only.
037: *          = 1: Compute singular vectors in factored form as well.
038: *
039: *  K       (input) INTEGER
040: *          The number of terms in the rational function to be solved
041: *          by SLASD4.  K >= 1.
042: *
043: *  D       (output) REAL array, dimension ( K )
044: *          On output, D contains the updated singular values.
045: *
046: *  Z       (input/output) REAL array, dimension ( K )
047: *          On entry, the first K elements of this array contain the
048: *          components of the deflation-adjusted updating row vector.
049: *          On exit, Z is updated.
050: *
051: *  VF      (input/output) REAL array, dimension ( K )
052: *          On entry, VF contains  information passed through DBEDE8.
053: *          On exit, VF contains the first K components of the first
054: *          components of all right singular vectors of the bidiagonal
055: *          matrix.
056: *
057: *  VL      (input/output) REAL array, dimension ( K )
058: *          On entry, VL contains  information passed through DBEDE8.
059: *          On exit, VL contains the first K components of the last
060: *          components of all right singular vectors of the bidiagonal
061: *          matrix.
062: *
063: *  DIFL    (output) REAL array, dimension ( K )
064: *          On exit, DIFL(I) = D(I) - DSIGMA(I).
065: *
066: *  DIFR    (output) REAL array,
067: *                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
068: *                   dimension ( K ) if ICOMPQ = 0.
069: *          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
070: *          defined and will not be referenced.
071: *
072: *          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
073: *          normalizing factors for the right singular vector matrix.
074: *
075: *  LDDIFR  (input) INTEGER
076: *          The leading dimension of DIFR, must be at least K.
077: *
078: *  DSIGMA  (input/output) REAL array, dimension ( K )
079: *          On entry, the first K elements of this array contain the old
080: *          roots of the deflated updating problem.  These are the poles
081: *          of the secular equation.
082: *          On exit, the elements of DSIGMA may be very slightly altered
083: *          in value.
084: *
085: *  WORK    (workspace) REAL array, dimension at least 3 * K
086: *
087: *  INFO    (output) INTEGER
088: *          = 0:  successful exit.
089: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
090: *          > 0:  if INFO = 1, an singular value did not converge
091: *
092: *  Further Details
093: *  ===============
094: *
095: *  Based on contributions by
096: *     Ming Gu and Huan Ren, Computer Science Division, University of
097: *     California at Berkeley, USA
098: *
099: *  =====================================================================
100: *
101: *     .. Parameters ..
102:       REAL               ONE
103:       PARAMETER          ( ONE = 1.0E+0 )
104: *     ..
105: *     .. Local Scalars ..
106:       INTEGER            I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
107:       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
108: *     ..
109: *     .. External Subroutines ..
110:       EXTERNAL           SCOPY, SLASCL, SLASD4, SLASET, XERBLA
111: *     ..
112: *     .. External Functions ..
113:       REAL               SDOT, SLAMC3, SNRM2
114:       EXTERNAL           SDOT, SLAMC3, SNRM2
115: *     ..
116: *     .. Intrinsic Functions ..
117:       INTRINSIC          ABS, SIGN, SQRT
118: *     ..
119: *     .. Executable Statements ..
120: *
121: *     Test the input parameters.
122: *
123:       INFO = 0
124: *
125:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
126:          INFO = -1
127:       ELSE IF( K.LT.1 ) THEN
128:          INFO = -2
129:       ELSE IF( LDDIFR.LT.K ) THEN
130:          INFO = -9
131:       END IF
132:       IF( INFO.NE.0 ) THEN
133:          CALL XERBLA( 'SLASD8', -INFO )
134:          RETURN
135:       END IF
136: *
137: *     Quick return if possible
138: *
139:       IF( K.EQ.1 ) THEN
140:          D( 1 ) = ABS( Z( 1 ) )
141:          DIFL( 1 ) = D( 1 )
142:          IF( ICOMPQ.EQ.1 ) THEN
143:             DIFL( 2 ) = ONE
144:             DIFR( 1, 2 ) = ONE
145:          END IF
146:          RETURN
147:       END IF
148: *
149: *     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
150: *     be computed with high relative accuracy (barring over/underflow).
151: *     This is a problem on machines without a guard digit in
152: *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
153: *     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
154: *     which on any of these machines zeros out the bottommost
155: *     bit of DSIGMA(I) if it is 1; this makes the subsequent
156: *     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
157: *     occurs. On binary machines with a guard digit (almost all
158: *     machines) it does not change DSIGMA(I) at all. On hexadecimal
159: *     and decimal machines with a guard digit, it slightly
160: *     changes the bottommost bits of DSIGMA(I). It does not account
161: *     for hexadecimal or decimal machines without guard digits
162: *     (we know of none). We use a subroutine call to compute
163: *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
164: *     this code.
165: *
166:       DO 10 I = 1, K
167:          DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I )
168:    10 CONTINUE
169: *
170: *     Book keeping.
171: *
172:       IWK1 = 1
173:       IWK2 = IWK1 + K
174:       IWK3 = IWK2 + K
175:       IWK2I = IWK2 - 1
176:       IWK3I = IWK3 - 1
177: *
178: *     Normalize Z.
179: *
180:       RHO = SNRM2( K, Z, 1 )
181:       CALL SLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
182:       RHO = RHO*RHO
183: *
184: *     Initialize WORK(IWK3).
185: *
186:       CALL SLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K )
187: *
188: *     Compute the updated singular values, the arrays DIFL, DIFR,
189: *     and the updated Z.
190: *
191:       DO 40 J = 1, K
192:          CALL SLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ),
193:      $                WORK( IWK2 ), INFO )
194: *
195: *        If the root finder fails, the computation is terminated.
196: *
197:          IF( INFO.NE.0 ) THEN
198:             RETURN
199:          END IF
200:          WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
201:          DIFL( J ) = -WORK( J )
202:          DIFR( J, 1 ) = -WORK( J+1 )
203:          DO 20 I = 1, J - 1
204:             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
205:      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
206:      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
207:      $                        DSIGMA( J ) )
208:    20    CONTINUE
209:          DO 30 I = J + 1, K
210:             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
211:      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
212:      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
213:      $                        DSIGMA( J ) )
214:    30    CONTINUE
215:    40 CONTINUE
216: *
217: *     Compute updated Z.
218: *
219:       DO 50 I = 1, K
220:          Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
221:    50 CONTINUE
222: *
223: *     Update VF and VL.
224: *
225:       DO 80 J = 1, K
226:          DIFLJ = DIFL( J )
227:          DJ = D( J )
228:          DSIGJ = -DSIGMA( J )
229:          IF( J.LT.K ) THEN
230:             DIFRJ = -DIFR( J, 1 )
231:             DSIGJP = -DSIGMA( J+1 )
232:          END IF
233:          WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
234:          DO 60 I = 1, J - 1
235:             WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
236:      $                   / ( DSIGMA( I )+DJ )
237:    60    CONTINUE
238:          DO 70 I = J + 1, K
239:             WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
240:      $                   / ( DSIGMA( I )+DJ )
241:    70    CONTINUE
242:          TEMP = SNRM2( K, WORK, 1 )
243:          WORK( IWK2I+J ) = SDOT( K, WORK, 1, VF, 1 ) / TEMP
244:          WORK( IWK3I+J ) = SDOT( K, WORK, 1, VL, 1 ) / TEMP
245:          IF( ICOMPQ.EQ.1 ) THEN
246:             DIFR( J, 2 ) = TEMP
247:          END IF
248:    80 CONTINUE
249: *
250:       CALL SCOPY( K, WORK( IWK2 ), 1, VF, 1 )
251:       CALL SCOPY( K, WORK( IWK3 ), 1, VL, 1 )
252: *
253:       RETURN
254: *
255: *     End of SLASD8
256: *
257:       END
258: 
259: