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