|
7072
|
1 SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, |
|
|
2 $ DSIGMA, WORK, INFO ) |
|
|
3 * |
|
|
4 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 INTEGER ICOMPQ, INFO, K, LDDIFR |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ), |
|
|
13 $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ), |
|
|
14 $ Z( * ) |
|
|
15 * .. |
|
|
16 * |
|
|
17 * Purpose |
|
|
18 * ======= |
|
|
19 * |
|
|
20 * DLASD8 finds the square roots of the roots of the secular equation, |
|
|
21 * as defined by the values in DSIGMA and Z. It makes the appropriate |
|
|
22 * calls to DLASD4, and stores, for each element in D, the distance |
|
|
23 * to its two nearest poles (elements in DSIGMA). It also updates |
|
|
24 * the arrays VF and VL, the first and last components of all the |
|
|
25 * right singular vectors of the original bidiagonal matrix. |
|
|
26 * |
|
|
27 * DLASD8 is called from DLASD6. |
|
|
28 * |
|
|
29 * Arguments |
|
|
30 * ========= |
|
|
31 * |
|
|
32 * ICOMPQ (input) INTEGER |
|
|
33 * Specifies whether singular vectors are to be computed in |
|
|
34 * factored form in the calling routine: |
|
|
35 * = 0: Compute singular values only. |
|
|
36 * = 1: Compute singular vectors in factored form as well. |
|
|
37 * |
|
|
38 * K (input) INTEGER |
|
|
39 * The number of terms in the rational function to be solved |
|
|
40 * by DLASD4. K >= 1. |
|
|
41 * |
|
|
42 * D (output) DOUBLE PRECISION array, dimension ( K ) |
|
|
43 * On output, D contains the updated singular values. |
|
|
44 * |
|
|
45 * Z (input) DOUBLE PRECISION array, dimension ( K ) |
|
|
46 * The first K elements of this array contain the components |
|
|
47 * of the deflation-adjusted updating row vector. |
|
|
48 * |
|
|
49 * VF (input/output) DOUBLE PRECISION array, dimension ( K ) |
|
|
50 * On entry, VF contains information passed through DBEDE8. |
|
|
51 * On exit, VF contains the first K components of the first |
|
|
52 * components of all right singular vectors of the bidiagonal |
|
|
53 * matrix. |
|
|
54 * |
|
|
55 * VL (input/output) DOUBLE PRECISION array, dimension ( K ) |
|
|
56 * On entry, VL contains information passed through DBEDE8. |
|
|
57 * On exit, VL contains the first K components of the last |
|
|
58 * components of all right singular vectors of the bidiagonal |
|
|
59 * matrix. |
|
|
60 * |
|
|
61 * DIFL (output) DOUBLE PRECISION array, dimension ( K ) |
|
|
62 * On exit, DIFL(I) = D(I) - DSIGMA(I). |
|
|
63 * |
|
|
64 * DIFR (output) DOUBLE PRECISION array, |
|
|
65 * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and |
|
|
66 * dimension ( K ) if ICOMPQ = 0. |
|
|
67 * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not |
|
|
68 * defined and will not be referenced. |
|
|
69 * |
|
|
70 * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the |
|
|
71 * normalizing factors for the right singular vector matrix. |
|
|
72 * |
|
|
73 * LDDIFR (input) INTEGER |
|
|
74 * The leading dimension of DIFR, must be at least K. |
|
|
75 * |
|
|
76 * DSIGMA (input) DOUBLE PRECISION array, dimension ( K ) |
|
|
77 * The first K elements of this array contain the old roots |
|
|
78 * of the deflated updating problem. These are the poles |
|
|
79 * of the secular equation. |
|
|
80 * |
|
|
81 * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K |
|
|
82 * |
|
|
83 * INFO (output) INTEGER |
|
|
84 * = 0: successful exit. |
|
|
85 * < 0: if INFO = -i, the i-th argument had an illegal value. |
|
|
86 * > 0: if INFO = 1, an singular value did not converge |
|
|
87 * |
|
|
88 * Further Details |
|
|
89 * =============== |
|
|
90 * |
|
|
91 * Based on contributions by |
|
|
92 * Ming Gu and Huan Ren, Computer Science Division, University of |
|
|
93 * California at Berkeley, USA |
|
|
94 * |
|
|
95 * ===================================================================== |
|
|
96 * |
|
|
97 * .. Parameters .. |
|
|
98 DOUBLE PRECISION ONE |
|
|
99 PARAMETER ( ONE = 1.0D+0 ) |
|
|
100 * .. |
|
|
101 * .. Local Scalars .. |
|
|
102 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J |
|
|
103 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP |
|
|
104 * .. |
|
|
105 * .. External Subroutines .. |
|
|
106 EXTERNAL DCOPY, DLASCL, DLASD4, DLASET, XERBLA |
|
|
107 * .. |
|
|
108 * .. External Functions .. |
|
|
109 DOUBLE PRECISION DDOT, DLAMC3, DNRM2 |
|
|
110 EXTERNAL DDOT, DLAMC3, DNRM2 |
|
|
111 * .. |
|
|
112 * .. Intrinsic Functions .. |
|
|
113 INTRINSIC ABS, SIGN, SQRT |
|
|
114 * .. |
|
|
115 * .. Executable Statements .. |
|
|
116 * |
|
|
117 * Test the input parameters. |
|
|
118 * |
|
|
119 INFO = 0 |
|
|
120 * |
|
|
121 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN |
|
|
122 INFO = -1 |
|
|
123 ELSE IF( K.LT.1 ) THEN |
|
|
124 INFO = -2 |
|
|
125 ELSE IF( LDDIFR.LT.K ) THEN |
|
|
126 INFO = -9 |
|
|
127 END IF |
|
|
128 IF( INFO.NE.0 ) THEN |
|
|
129 CALL XERBLA( 'DLASD8', -INFO ) |
|
|
130 RETURN |
|
|
131 END IF |
|
|
132 * |
|
|
133 * Quick return if possible |
|
|
134 * |
|
|
135 IF( K.EQ.1 ) THEN |
|
|
136 D( 1 ) = ABS( Z( 1 ) ) |
|
|
137 DIFL( 1 ) = D( 1 ) |
|
|
138 IF( ICOMPQ.EQ.1 ) THEN |
|
|
139 DIFL( 2 ) = ONE |
|
|
140 DIFR( 1, 2 ) = ONE |
|
|
141 END IF |
|
|
142 RETURN |
|
|
143 END IF |
|
|
144 * |
|
|
145 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can |
|
|
146 * be computed with high relative accuracy (barring over/underflow). |
|
|
147 * This is a problem on machines without a guard digit in |
|
|
148 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). |
|
|
149 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), |
|
|
150 * which on any of these machines zeros out the bottommost |
|
|
151 * bit of DSIGMA(I) if it is 1; this makes the subsequent |
|
|
152 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation |
|
|
153 * occurs. On binary machines with a guard digit (almost all |
|
|
154 * machines) it does not change DSIGMA(I) at all. On hexadecimal |
|
|
155 * and decimal machines with a guard digit, it slightly |
|
|
156 * changes the bottommost bits of DSIGMA(I). It does not account |
|
|
157 * for hexadecimal or decimal machines without guard digits |
|
|
158 * (we know of none). We use a subroutine call to compute |
|
|
159 * 2*DSIGMA(I) to prevent optimizing compilers from eliminating |
|
|
160 * this code. |
|
|
161 * |
|
|
162 DO 10 I = 1, K |
|
|
163 DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) |
|
|
164 10 CONTINUE |
|
|
165 * |
|
|
166 * Book keeping. |
|
|
167 * |
|
|
168 IWK1 = 1 |
|
|
169 IWK2 = IWK1 + K |
|
|
170 IWK3 = IWK2 + K |
|
|
171 IWK2I = IWK2 - 1 |
|
|
172 IWK3I = IWK3 - 1 |
|
|
173 * |
|
|
174 * Normalize Z. |
|
|
175 * |
|
|
176 RHO = DNRM2( K, Z, 1 ) |
|
|
177 CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO ) |
|
|
178 RHO = RHO*RHO |
|
|
179 * |
|
|
180 * Initialize WORK(IWK3). |
|
|
181 * |
|
|
182 CALL DLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K ) |
|
|
183 * |
|
|
184 * Compute the updated singular values, the arrays DIFL, DIFR, |
|
|
185 * and the updated Z. |
|
|
186 * |
|
|
187 DO 40 J = 1, K |
|
|
188 CALL DLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ), |
|
|
189 $ WORK( IWK2 ), INFO ) |
|
|
190 * |
|
|
191 * If the root finder fails, the computation is terminated. |
|
|
192 * |
|
|
193 IF( INFO.NE.0 ) THEN |
|
|
194 RETURN |
|
|
195 END IF |
|
|
196 WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J ) |
|
|
197 DIFL( J ) = -WORK( J ) |
|
|
198 DIFR( J, 1 ) = -WORK( J+1 ) |
|
|
199 DO 20 I = 1, J - 1 |
|
|
200 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )* |
|
|
201 $ WORK( IWK2I+I ) / ( DSIGMA( I )- |
|
|
202 $ DSIGMA( J ) ) / ( DSIGMA( I )+ |
|
|
203 $ DSIGMA( J ) ) |
|
|
204 20 CONTINUE |
|
|
205 DO 30 I = J + 1, K |
|
|
206 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )* |
|
|
207 $ WORK( IWK2I+I ) / ( DSIGMA( I )- |
|
|
208 $ DSIGMA( J ) ) / ( DSIGMA( I )+ |
|
|
209 $ DSIGMA( J ) ) |
|
|
210 30 CONTINUE |
|
|
211 40 CONTINUE |
|
|
212 * |
|
|
213 * Compute updated Z. |
|
|
214 * |
|
|
215 DO 50 I = 1, K |
|
|
216 Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) ) |
|
|
217 50 CONTINUE |
|
|
218 * |
|
|
219 * Update VF and VL. |
|
|
220 * |
|
|
221 DO 80 J = 1, K |
|
|
222 DIFLJ = DIFL( J ) |
|
|
223 DJ = D( J ) |
|
|
224 DSIGJ = -DSIGMA( J ) |
|
|
225 IF( J.LT.K ) THEN |
|
|
226 DIFRJ = -DIFR( J, 1 ) |
|
|
227 DSIGJP = -DSIGMA( J+1 ) |
|
|
228 END IF |
|
|
229 WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ ) |
|
|
230 DO 60 I = 1, J - 1 |
|
|
231 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ ) |
|
|
232 $ / ( DSIGMA( I )+DJ ) |
|
|
233 60 CONTINUE |
|
|
234 DO 70 I = J + 1, K |
|
|
235 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ ) |
|
|
236 $ / ( DSIGMA( I )+DJ ) |
|
|
237 70 CONTINUE |
|
|
238 TEMP = DNRM2( K, WORK, 1 ) |
|
|
239 WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP |
|
|
240 WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP |
|
|
241 IF( ICOMPQ.EQ.1 ) THEN |
|
|
242 DIFR( J, 2 ) = TEMP |
|
|
243 END IF |
|
|
244 80 CONTINUE |
|
|
245 * |
|
|
246 CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 ) |
|
|
247 CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 ) |
|
|
248 * |
|
|
249 RETURN |
|
|
250 * |
|
|
251 * End of DLASD8 |
|
|
252 * |
|
|
253 END |
