LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlasd8.f
Go to the documentation of this file.
1 *> \brief \b DLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLASD8 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd8.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd8.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd8.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
22 * DSIGMA, WORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER ICOMPQ, INFO, K, LDDIFR
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
29 * $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
30 * $ Z( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DLASD8 finds the square roots of the roots of the secular equation,
40 *> as defined by the values in DSIGMA and Z. It makes the appropriate
41 *> calls to DLASD4, and stores, for each element in D, the distance
42 *> to its two nearest poles (elements in DSIGMA). It also updates
43 *> the arrays VF and VL, the first and last components of all the
44 *> right singular vectors of the original bidiagonal matrix.
45 *>
46 *> DLASD8 is called from DLASD6.
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \param[in] ICOMPQ
53 *> \verbatim
54 *> ICOMPQ is INTEGER
55 *> Specifies whether singular vectors are to be computed in
56 *> factored form in the calling routine:
57 *> = 0: Compute singular values only.
58 *> = 1: Compute singular vectors in factored form as well.
59 *> \endverbatim
60 *>
61 *> \param[in] K
62 *> \verbatim
63 *> K is INTEGER
64 *> The number of terms in the rational function to be solved
65 *> by DLASD4. K >= 1.
66 *> \endverbatim
67 *>
68 *> \param[out] D
69 *> \verbatim
70 *> D is DOUBLE PRECISION array, dimension ( K )
71 *> On output, D contains the updated singular values.
72 *> \endverbatim
73 *>
74 *> \param[in,out] Z
75 *> \verbatim
76 *> Z is DOUBLE PRECISION array, dimension ( K )
77 *> On entry, the first K elements of this array contain the
78 *> components of the deflation-adjusted updating row vector.
79 *> On exit, Z is updated.
80 *> \endverbatim
81 *>
82 *> \param[in,out] VF
83 *> \verbatim
84 *> VF is DOUBLE PRECISION array, dimension ( K )
85 *> On entry, VF contains information passed through DBEDE8.
86 *> On exit, VF contains the first K components of the first
87 *> components of all right singular vectors of the bidiagonal
88 *> matrix.
89 *> \endverbatim
90 *>
91 *> \param[in,out] VL
92 *> \verbatim
93 *> VL is DOUBLE PRECISION array, dimension ( K )
94 *> On entry, VL contains information passed through DBEDE8.
95 *> On exit, VL contains the first K components of the last
96 *> components of all right singular vectors of the bidiagonal
97 *> matrix.
98 *> \endverbatim
99 *>
100 *> \param[out] DIFL
101 *> \verbatim
102 *> DIFL is DOUBLE PRECISION array, dimension ( K )
103 *> On exit, DIFL(I) = D(I) - DSIGMA(I).
104 *> \endverbatim
105 *>
106 *> \param[out] DIFR
107 *> \verbatim
108 *> DIFR is DOUBLE PRECISION array,
109 *> dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
110 *> dimension ( K ) if ICOMPQ = 0.
111 *> On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
112 *> defined and will not be referenced.
113 *>
114 *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
115 *> normalizing factors for the right singular vector matrix.
116 *> \endverbatim
117 *>
118 *> \param[in] LDDIFR
119 *> \verbatim
120 *> LDDIFR is INTEGER
121 *> The leading dimension of DIFR, must be at least K.
122 *> \endverbatim
123 *>
124 *> \param[in,out] DSIGMA
125 *> \verbatim
126 *> DSIGMA is DOUBLE PRECISION array, dimension ( K )
127 *> On entry, the first K elements of this array contain the old
128 *> roots of the deflated updating problem. These are the poles
129 *> of the secular equation.
130 *> On exit, the elements of DSIGMA may be very slightly altered
131 *> in value.
132 *> \endverbatim
133 *>
134 *> \param[out] WORK
135 *> \verbatim
136 *> WORK is DOUBLE PRECISION array, dimension at least 3 * K
137 *> \endverbatim
138 *>
139 *> \param[out] INFO
140 *> \verbatim
141 *> INFO is INTEGER
142 *> = 0: successful exit.
143 *> < 0: if INFO = -i, the i-th argument had an illegal value.
144 *> > 0: if INFO = 1, a singular value did not converge
145 *> \endverbatim
146 *
147 * Authors:
148 * ========
149 *
150 *> \author Univ. of Tennessee
151 *> \author Univ. of California Berkeley
152 *> \author Univ. of Colorado Denver
153 *> \author NAG Ltd.
154 *
155 *> \date September 2012
156 *
157 *> \ingroup auxOTHERauxiliary
158 *
159 *> \par Contributors:
160 * ==================
161 *>
162 *> Ming Gu and Huan Ren, Computer Science Division, University of
163 *> California at Berkeley, USA
164 *>
165 * =====================================================================
166  SUBROUTINE dlasd8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
167  $ dsigma, work, info )
168 *
169 * -- LAPACK auxiliary routine (version 3.4.2) --
170 * -- LAPACK is a software package provided by Univ. of Tennessee, --
171 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172 * September 2012
173 *
174 * .. Scalar Arguments ..
175  INTEGER icompq, info, k, lddifr
176 * ..
177 * .. Array Arguments ..
178  DOUBLE PRECISION d( * ), difl( * ), difr( lddifr, * ),
179  $ dsigma( * ), vf( * ), vl( * ), work( * ),
180  $ z( * )
181 * ..
182 *
183 * =====================================================================
184 *
185 * .. Parameters ..
186  DOUBLE PRECISION one
187  parameter( one = 1.0d+0 )
188 * ..
189 * .. Local Scalars ..
190  INTEGER i, iwk1, iwk2, iwk2i, iwk3, iwk3i, j
191  DOUBLE PRECISION diflj, difrj, dj, dsigj, dsigjp, rho, temp
192 * ..
193 * .. External Subroutines ..
194  EXTERNAL dcopy, dlascl, dlasd4, dlaset, xerbla
195 * ..
196 * .. External Functions ..
197  DOUBLE PRECISION ddot, dlamc3, dnrm2
198  EXTERNAL ddot, dlamc3, dnrm2
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC abs, sign, sqrt
202 * ..
203 * .. Executable Statements ..
204 *
205 * Test the input parameters.
206 *
207  info = 0
208 *
209  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
210  info = -1
211  ELSE IF( k.LT.1 ) THEN
212  info = -2
213  ELSE IF( lddifr.LT.k ) THEN
214  info = -9
215  END IF
216  IF( info.NE.0 ) THEN
217  CALL xerbla( 'DLASD8', -info )
218  return
219  END IF
220 *
221 * Quick return if possible
222 *
223  IF( k.EQ.1 ) THEN
224  d( 1 ) = abs( z( 1 ) )
225  difl( 1 ) = d( 1 )
226  IF( icompq.EQ.1 ) THEN
227  difl( 2 ) = one
228  difr( 1, 2 ) = one
229  END IF
230  return
231  END IF
232 *
233 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
234 * be computed with high relative accuracy (barring over/underflow).
235 * This is a problem on machines without a guard digit in
236 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
237 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
238 * which on any of these machines zeros out the bottommost
239 * bit of DSIGMA(I) if it is 1; this makes the subsequent
240 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
241 * occurs. On binary machines with a guard digit (almost all
242 * machines) it does not change DSIGMA(I) at all. On hexadecimal
243 * and decimal machines with a guard digit, it slightly
244 * changes the bottommost bits of DSIGMA(I). It does not account
245 * for hexadecimal or decimal machines without guard digits
246 * (we know of none). We use a subroutine call to compute
247 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
248 * this code.
249 *
250  DO 10 i = 1, k
251  dsigma( i ) = dlamc3( dsigma( i ), dsigma( i ) ) - dsigma( i )
252  10 continue
253 *
254 * Book keeping.
255 *
256  iwk1 = 1
257  iwk2 = iwk1 + k
258  iwk3 = iwk2 + k
259  iwk2i = iwk2 - 1
260  iwk3i = iwk3 - 1
261 *
262 * Normalize Z.
263 *
264  rho = dnrm2( k, z, 1 )
265  CALL dlascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
266  rho = rho*rho
267 *
268 * Initialize WORK(IWK3).
269 *
270  CALL dlaset( 'A', k, 1, one, one, work( iwk3 ), k )
271 *
272 * Compute the updated singular values, the arrays DIFL, DIFR,
273 * and the updated Z.
274 *
275  DO 40 j = 1, k
276  CALL dlasd4( k, j, dsigma, z, work( iwk1 ), rho, d( j ),
277  $ work( iwk2 ), info )
278 *
279 * If the root finder fails, the computation is terminated.
280 *
281  IF( info.NE.0 ) THEN
282  CALL xerbla( 'DLASD4', -info )
283  return
284  END IF
285  work( iwk3i+j ) = work( iwk3i+j )*work( j )*work( iwk2i+j )
286  difl( j ) = -work( j )
287  difr( j, 1 ) = -work( j+1 )
288  DO 20 i = 1, j - 1
289  work( iwk3i+i ) = work( iwk3i+i )*work( i )*
290  $ work( iwk2i+i ) / ( dsigma( i )-
291  $ dsigma( j ) ) / ( dsigma( i )+
292  $ dsigma( j ) )
293  20 continue
294  DO 30 i = j + 1, k
295  work( iwk3i+i ) = work( iwk3i+i )*work( i )*
296  $ work( iwk2i+i ) / ( dsigma( i )-
297  $ dsigma( j ) ) / ( dsigma( i )+
298  $ dsigma( j ) )
299  30 continue
300  40 continue
301 *
302 * Compute updated Z.
303 *
304  DO 50 i = 1, k
305  z( i ) = sign( sqrt( abs( work( iwk3i+i ) ) ), z( i ) )
306  50 continue
307 *
308 * Update VF and VL.
309 *
310  DO 80 j = 1, k
311  diflj = difl( j )
312  dj = d( j )
313  dsigj = -dsigma( j )
314  IF( j.LT.k ) THEN
315  difrj = -difr( j, 1 )
316  dsigjp = -dsigma( j+1 )
317  END IF
318  work( j ) = -z( j ) / diflj / ( dsigma( j )+dj )
319  DO 60 i = 1, j - 1
320  work( i ) = z( i ) / ( dlamc3( dsigma( i ), dsigj )-diflj )
321  $ / ( dsigma( i )+dj )
322  60 continue
323  DO 70 i = j + 1, k
324  work( i ) = z( i ) / ( dlamc3( dsigma( i ), dsigjp )+difrj )
325  $ / ( dsigma( i )+dj )
326  70 continue
327  temp = dnrm2( k, work, 1 )
328  work( iwk2i+j ) = ddot( k, work, 1, vf, 1 ) / temp
329  work( iwk3i+j ) = ddot( k, work, 1, vl, 1 ) / temp
330  IF( icompq.EQ.1 ) THEN
331  difr( j, 2 ) = temp
332  END IF
333  80 continue
334 *
335  CALL dcopy( k, work( iwk2 ), 1, vf, 1 )
336  CALL dcopy( k, work( iwk3 ), 1, vl, 1 )
337 *
338  return
339 *
340 * End of DLASD8
341 *
342  END
343