LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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] 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*> \endverbatim
131*>
132*> \param[out] WORK
133*> \verbatim
134*> WORK is DOUBLE PRECISION array, dimension (3*K)
135*> \endverbatim
136*>
137*> \param[out] INFO
138*> \verbatim
139*> INFO is INTEGER
140*> = 0: successful exit.
141*> < 0: if INFO = -i, the i-th argument had an illegal value.
142*> > 0: if INFO = 1, a singular value did not converge
143*> \endverbatim
144*
145* Authors:
146* ========
147*
148*> \author Univ. of Tennessee
149*> \author Univ. of California Berkeley
150*> \author Univ. of Colorado Denver
151*> \author NAG Ltd.
152*
153*> \ingroup lasd8
154*
155*> \par Contributors:
156* ==================
157*>
158*> Ming Gu and Huan Ren, Computer Science Division, University of
159*> California at Berkeley, USA
160*>
161* =====================================================================
162 SUBROUTINE dlasd8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
163 $ DSIGMA, WORK, INFO )
164*
165* -- LAPACK auxiliary routine --
166* -- LAPACK is a software package provided by Univ. of Tennessee, --
167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168*
169* .. Scalar Arguments ..
170 INTEGER ICOMPQ, INFO, K, LDDIFR
171* ..
172* .. Array Arguments ..
173 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
174 $ dsigma( * ), vf( * ), vl( * ), work( * ),
175 $ z( * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 DOUBLE PRECISION ONE
182 parameter( one = 1.0d+0 )
183* ..
184* .. Local Scalars ..
185 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
186 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
187* ..
188* .. External Subroutines ..
189 EXTERNAL dcopy, dlascl, dlasd4, dlaset, xerbla
190* ..
191* .. External Functions ..
192 DOUBLE PRECISION DDOT, DLAMC3, DNRM2
193 EXTERNAL ddot, dlamc3, dnrm2
194* ..
195* .. Intrinsic Functions ..
196 INTRINSIC abs, sign, sqrt
197* ..
198* .. Executable Statements ..
199*
200* Test the input parameters.
201*
202 info = 0
203*
204 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
205 info = -1
206 ELSE IF( k.LT.1 ) THEN
207 info = -2
208 ELSE IF( lddifr.LT.k ) THEN
209 info = -9
210 END IF
211 IF( info.NE.0 ) THEN
212 CALL xerbla( 'DLASD8', -info )
213 RETURN
214 END IF
215*
216* Quick return if possible
217*
218 IF( k.EQ.1 ) THEN
219 d( 1 ) = abs( z( 1 ) )
220 difl( 1 ) = d( 1 )
221 IF( icompq.EQ.1 ) THEN
222 difl( 2 ) = one
223 difr( 1, 2 ) = one
224 END IF
225 RETURN
226 END IF
227*
228* Book keeping.
229*
230 iwk1 = 1
231 iwk2 = iwk1 + k
232 iwk3 = iwk2 + k
233 iwk2i = iwk2 - 1
234 iwk3i = iwk3 - 1
235*
236* Normalize Z.
237*
238 rho = dnrm2( k, z, 1 )
239 CALL dlascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
240 rho = rho*rho
241*
242* Initialize WORK(IWK3).
243*
244 CALL dlaset( 'A', k, 1, one, one, work( iwk3 ), k )
245*
246* Compute the updated singular values, the arrays DIFL, DIFR,
247* and the updated Z.
248*
249 DO 40 j = 1, k
250 CALL dlasd4( k, j, dsigma, z, work( iwk1 ), rho, d( j ),
251 $ work( iwk2 ), info )
252*
253* If the root finder fails, report the convergence failure.
254*
255 IF( info.NE.0 ) THEN
256 RETURN
257 END IF
258 work( iwk3i+j ) = work( iwk3i+j )*work( j )*work( iwk2i+j )
259 difl( j ) = -work( j )
260 difr( j, 1 ) = -work( j+1 )
261 DO 20 i = 1, j - 1
262 work( iwk3i+i ) = work( iwk3i+i )*work( i )*
263 $ work( iwk2i+i ) / ( dsigma( i )-
264 $ dsigma( j ) ) / ( dsigma( i )+
265 $ dsigma( j ) )
266 20 CONTINUE
267 DO 30 i = j + 1, k
268 work( iwk3i+i ) = work( iwk3i+i )*work( i )*
269 $ work( iwk2i+i ) / ( dsigma( i )-
270 $ dsigma( j ) ) / ( dsigma( i )+
271 $ dsigma( j ) )
272 30 CONTINUE
273 40 CONTINUE
274*
275* Compute updated Z.
276*
277 DO 50 i = 1, k
278 z( i ) = sign( sqrt( abs( work( iwk3i+i ) ) ), z( i ) )
279 50 CONTINUE
280*
281* Update VF and VL.
282*
283 DO 80 j = 1, k
284 diflj = difl( j )
285 dj = d( j )
286 dsigj = -dsigma( j )
287 IF( j.LT.k ) THEN
288 difrj = -difr( j, 1 )
289 dsigjp = -dsigma( j+1 )
290 END IF
291 work( j ) = -z( j ) / diflj / ( dsigma( j )+dj )
292*
293* Use calls to the subroutine DLAMC3 to enforce the parentheses
294* (x+y)+z. The goal is to prevent optimizing compilers
295* from doing x+(y+z).
296*
297 DO 60 i = 1, j - 1
298 work( i ) = z( i ) / ( dlamc3( dsigma( i ), dsigj )-diflj )
299 $ / ( dsigma( i )+dj )
300 60 CONTINUE
301 DO 70 i = j + 1, k
302 work( i ) = z( i ) / ( dlamc3( dsigma( i ), dsigjp )+difrj )
303 $ / ( dsigma( i )+dj )
304 70 CONTINUE
305 temp = dnrm2( k, work, 1 )
306 work( iwk2i+j ) = ddot( k, work, 1, vf, 1 ) / temp
307 work( iwk3i+j ) = ddot( k, work, 1, vl, 1 ) / temp
308 IF( icompq.EQ.1 ) THEN
309 difr( j, 2 ) = temp
310 END IF
311 80 CONTINUE
312*
313 CALL dcopy( k, work( iwk2 ), 1, vf, 1 )
314 CALL dcopy( k, work( iwk3 ), 1, vl, 1 )
315*
316 RETURN
317*
318* End of DLASD8
319*
320 END
321
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dlasd4(n, i, d, z, delta, rho, sigma, work, info)
DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...
Definition dlasd4.f:153
subroutine dlasd8(icompq, k, d, z, vf, vl, difl, difr, lddifr, dsigma, work, info)
DLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D...
Definition dlasd8.f:164
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110