1 SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
2 $ DSIGMA, WORK, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.0) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2010
8 *
9 * .. Scalar Arguments ..
10 INTEGER ICOMPQ, INFO, K, LDDIFR
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
14 $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
15 $ Z( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLASD8 finds the square roots of the roots of the secular equation,
22 * as defined by the values in DSIGMA and Z. It makes the appropriate
23 * calls to DLASD4, and stores, for each element in D, the distance
24 * to its two nearest poles (elements in DSIGMA). It also updates
25 * the arrays VF and VL, the first and last components of all the
26 * right singular vectors of the original bidiagonal matrix.
27 *
28 * DLASD8 is called from DLASD6.
29 *
30 * Arguments
31 * =========
32 *
33 * ICOMPQ (input) INTEGER
34 * Specifies whether singular vectors are to be computed in
35 * factored form in the calling routine:
36 * = 0: Compute singular values only.
37 * = 1: Compute singular vectors in factored form as well.
38 *
39 * K (input) INTEGER
40 * The number of terms in the rational function to be solved
41 * by DLASD4. K >= 1.
42 *
43 * D (output) DOUBLE PRECISION array, dimension ( K )
44 * On output, D contains the updated singular values.
45 *
46 * Z (input/output) DOUBLE PRECISION array, dimension ( K )
47 * On entry, the first K elements of this array contain the
48 * components of the deflation-adjusted updating row vector.
49 * On exit, Z is updated.
50 *
51 * VF (input/output) DOUBLE PRECISION array, dimension ( K )
52 * On entry, VF contains information passed through DBEDE8.
53 * On exit, VF contains the first K components of the first
54 * components of all right singular vectors of the bidiagonal
55 * matrix.
56 *
57 * VL (input/output) DOUBLE PRECISION array, dimension ( K )
58 * On entry, VL contains information passed through DBEDE8.
59 * On exit, VL contains the first K components of the last
60 * components of all right singular vectors of the bidiagonal
61 * matrix.
62 *
63 * DIFL (output) DOUBLE PRECISION array, dimension ( K )
64 * On exit, DIFL(I) = D(I) - DSIGMA(I).
65 *
66 * DIFR (output) DOUBLE PRECISION array,
67 * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
68 * dimension ( K ) if ICOMPQ = 0.
69 * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
70 * defined and will not be referenced.
71 *
72 * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
73 * normalizing factors for the right singular vector matrix.
74 *
75 * LDDIFR (input) INTEGER
76 * The leading dimension of DIFR, must be at least K.
77 *
78 * DSIGMA (input/output) DOUBLE PRECISION array, dimension ( K )
79 * On entry, the first K elements of this array contain the old
80 * roots of the deflated updating problem. These are the poles
81 * of the secular equation.
82 * On exit, the elements of DSIGMA may be very slightly altered
83 * in value.
84 *
85 * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K
86 *
87 * INFO (output) INTEGER
88 * = 0: successful exit.
89 * < 0: if INFO = -i, the i-th argument had an illegal value.
90 * > 0: if INFO = 1, a singular value did not converge
91 *
92 * Further Details
93 * ===============
94 *
95 * Based on contributions by
96 * Ming Gu and Huan Ren, Computer Science Division, University of
97 * California at Berkeley, USA
98 *
99 * =====================================================================
100 *
101 * .. Parameters ..
102 DOUBLE PRECISION ONE
103 PARAMETER ( ONE = 1.0D+0 )
104 * ..
105 * .. Local Scalars ..
106 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
107 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
108 * ..
109 * .. External Subroutines ..
110 EXTERNAL DCOPY, DLASCL, DLASD4, DLASET, XERBLA
111 * ..
112 * .. External Functions ..
113 DOUBLE PRECISION DDOT, DLAMC3, DNRM2
114 EXTERNAL DDOT, DLAMC3, DNRM2
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( 'DLASD8', -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 ) = DLAMC3( 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 = DNRM2( K, Z, 1 )
181 CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
182 RHO = RHO*RHO
183 *
184 * Initialize WORK(IWK3).
185 *
186 CALL DLASET( '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 DLASD4( 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 CALL XERBLA( 'DLASD4', -INFO )
199 RETURN
200 END IF
201 WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
202 DIFL( J ) = -WORK( J )
203 DIFR( J, 1 ) = -WORK( J+1 )
204 DO 20 I = 1, J - 1
205 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
206 $ WORK( IWK2I+I ) / ( DSIGMA( I )-
207 $ DSIGMA( J ) ) / ( DSIGMA( I )+
208 $ DSIGMA( J ) )
209 20 CONTINUE
210 DO 30 I = J + 1, K
211 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
212 $ WORK( IWK2I+I ) / ( DSIGMA( I )-
213 $ DSIGMA( J ) ) / ( DSIGMA( I )+
214 $ DSIGMA( J ) )
215 30 CONTINUE
216 40 CONTINUE
217 *
218 * Compute updated Z.
219 *
220 DO 50 I = 1, K
221 Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
222 50 CONTINUE
223 *
224 * Update VF and VL.
225 *
226 DO 80 J = 1, K
227 DIFLJ = DIFL( J )
228 DJ = D( J )
229 DSIGJ = -DSIGMA( J )
230 IF( J.LT.K ) THEN
231 DIFRJ = -DIFR( J, 1 )
232 DSIGJP = -DSIGMA( J+1 )
233 END IF
234 WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
235 DO 60 I = 1, J - 1
236 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
237 $ / ( DSIGMA( I )+DJ )
238 60 CONTINUE
239 DO 70 I = J + 1, K
240 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
241 $ / ( DSIGMA( I )+DJ )
242 70 CONTINUE
243 TEMP = DNRM2( K, WORK, 1 )
244 WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP
245 WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP
246 IF( ICOMPQ.EQ.1 ) THEN
247 DIFR( J, 2 ) = TEMP
248 END IF
249 80 CONTINUE
250 *
251 CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 )
252 CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 )
253 *
254 RETURN
255 *
256 * End of DLASD8
257 *
258 END
259
2 $ DSIGMA, WORK, INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.3.0) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2010
8 *
9 * .. Scalar Arguments ..
10 INTEGER ICOMPQ, INFO, K, LDDIFR
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
14 $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
15 $ Z( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLASD8 finds the square roots of the roots of the secular equation,
22 * as defined by the values in DSIGMA and Z. It makes the appropriate
23 * calls to DLASD4, and stores, for each element in D, the distance
24 * to its two nearest poles (elements in DSIGMA). It also updates
25 * the arrays VF and VL, the first and last components of all the
26 * right singular vectors of the original bidiagonal matrix.
27 *
28 * DLASD8 is called from DLASD6.
29 *
30 * Arguments
31 * =========
32 *
33 * ICOMPQ (input) INTEGER
34 * Specifies whether singular vectors are to be computed in
35 * factored form in the calling routine:
36 * = 0: Compute singular values only.
37 * = 1: Compute singular vectors in factored form as well.
38 *
39 * K (input) INTEGER
40 * The number of terms in the rational function to be solved
41 * by DLASD4. K >= 1.
42 *
43 * D (output) DOUBLE PRECISION array, dimension ( K )
44 * On output, D contains the updated singular values.
45 *
46 * Z (input/output) DOUBLE PRECISION array, dimension ( K )
47 * On entry, the first K elements of this array contain the
48 * components of the deflation-adjusted updating row vector.
49 * On exit, Z is updated.
50 *
51 * VF (input/output) DOUBLE PRECISION array, dimension ( K )
52 * On entry, VF contains information passed through DBEDE8.
53 * On exit, VF contains the first K components of the first
54 * components of all right singular vectors of the bidiagonal
55 * matrix.
56 *
57 * VL (input/output) DOUBLE PRECISION array, dimension ( K )
58 * On entry, VL contains information passed through DBEDE8.
59 * On exit, VL contains the first K components of the last
60 * components of all right singular vectors of the bidiagonal
61 * matrix.
62 *
63 * DIFL (output) DOUBLE PRECISION array, dimension ( K )
64 * On exit, DIFL(I) = D(I) - DSIGMA(I).
65 *
66 * DIFR (output) DOUBLE PRECISION array,
67 * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
68 * dimension ( K ) if ICOMPQ = 0.
69 * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
70 * defined and will not be referenced.
71 *
72 * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
73 * normalizing factors for the right singular vector matrix.
74 *
75 * LDDIFR (input) INTEGER
76 * The leading dimension of DIFR, must be at least K.
77 *
78 * DSIGMA (input/output) DOUBLE PRECISION array, dimension ( K )
79 * On entry, the first K elements of this array contain the old
80 * roots of the deflated updating problem. These are the poles
81 * of the secular equation.
82 * On exit, the elements of DSIGMA may be very slightly altered
83 * in value.
84 *
85 * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K
86 *
87 * INFO (output) INTEGER
88 * = 0: successful exit.
89 * < 0: if INFO = -i, the i-th argument had an illegal value.
90 * > 0: if INFO = 1, a singular value did not converge
91 *
92 * Further Details
93 * ===============
94 *
95 * Based on contributions by
96 * Ming Gu and Huan Ren, Computer Science Division, University of
97 * California at Berkeley, USA
98 *
99 * =====================================================================
100 *
101 * .. Parameters ..
102 DOUBLE PRECISION ONE
103 PARAMETER ( ONE = 1.0D+0 )
104 * ..
105 * .. Local Scalars ..
106 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
107 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
108 * ..
109 * .. External Subroutines ..
110 EXTERNAL DCOPY, DLASCL, DLASD4, DLASET, XERBLA
111 * ..
112 * .. External Functions ..
113 DOUBLE PRECISION DDOT, DLAMC3, DNRM2
114 EXTERNAL DDOT, DLAMC3, DNRM2
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( 'DLASD8', -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 ) = DLAMC3( 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 = DNRM2( K, Z, 1 )
181 CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
182 RHO = RHO*RHO
183 *
184 * Initialize WORK(IWK3).
185 *
186 CALL DLASET( '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 DLASD4( 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 CALL XERBLA( 'DLASD4', -INFO )
199 RETURN
200 END IF
201 WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
202 DIFL( J ) = -WORK( J )
203 DIFR( J, 1 ) = -WORK( J+1 )
204 DO 20 I = 1, J - 1
205 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
206 $ WORK( IWK2I+I ) / ( DSIGMA( I )-
207 $ DSIGMA( J ) ) / ( DSIGMA( I )+
208 $ DSIGMA( J ) )
209 20 CONTINUE
210 DO 30 I = J + 1, K
211 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
212 $ WORK( IWK2I+I ) / ( DSIGMA( I )-
213 $ DSIGMA( J ) ) / ( DSIGMA( I )+
214 $ DSIGMA( J ) )
215 30 CONTINUE
216 40 CONTINUE
217 *
218 * Compute updated Z.
219 *
220 DO 50 I = 1, K
221 Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
222 50 CONTINUE
223 *
224 * Update VF and VL.
225 *
226 DO 80 J = 1, K
227 DIFLJ = DIFL( J )
228 DJ = D( J )
229 DSIGJ = -DSIGMA( J )
230 IF( J.LT.K ) THEN
231 DIFRJ = -DIFR( J, 1 )
232 DSIGJP = -DSIGMA( J+1 )
233 END IF
234 WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
235 DO 60 I = 1, J - 1
236 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
237 $ / ( DSIGMA( I )+DJ )
238 60 CONTINUE
239 DO 70 I = J + 1, K
240 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
241 $ / ( DSIGMA( I )+DJ )
242 70 CONTINUE
243 TEMP = DNRM2( K, WORK, 1 )
244 WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP
245 WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP
246 IF( ICOMPQ.EQ.1 ) THEN
247 DIFR( J, 2 ) = TEMP
248 END IF
249 80 CONTINUE
250 *
251 CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 )
252 CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 )
253 *
254 RETURN
255 *
256 * End of DLASD8
257 *
258 END
259