LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slasd8.f
Go to the documentation of this file.
1*> \brief \b SLASD8 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*> Download SLASD8 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd8.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd8.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd8.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
20* DSIGMA, WORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER ICOMPQ, INFO, K, LDDIFR
24* ..
25* .. Array Arguments ..
26* REAL D( * ), DIFL( * ), DIFR( LDDIFR, * ),
27* $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
28* $ Z( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SLASD8 finds the square roots of the roots of the secular equation,
38*> as defined by the values in DSIGMA and Z. It makes the appropriate
39*> calls to SLASD4, and stores, for each element in D, the distance
40*> to its two nearest poles (elements in DSIGMA). It also updates
41*> the arrays VF and VL, the first and last components of all the
42*> right singular vectors of the original bidiagonal matrix.
43*>
44*> SLASD8 is called from SLASD6.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] ICOMPQ
51*> \verbatim
52*> ICOMPQ is INTEGER
53*> Specifies whether singular vectors are to be computed in
54*> factored form in the calling routine:
55*> = 0: Compute singular values only.
56*> = 1: Compute singular vectors in factored form as well.
57*> \endverbatim
58*>
59*> \param[in] K
60*> \verbatim
61*> K is INTEGER
62*> The number of terms in the rational function to be solved
63*> by SLASD4. K >= 1.
64*> \endverbatim
65*>
66*> \param[out] D
67*> \verbatim
68*> D is REAL array, dimension ( K )
69*> On output, D contains the updated singular values.
70*> \endverbatim
71*>
72*> \param[in,out] Z
73*> \verbatim
74*> Z is REAL array, dimension ( K )
75*> On entry, the first K elements of this array contain the
76*> components of the deflation-adjusted updating row vector.
77*> On exit, Z is updated.
78*> \endverbatim
79*>
80*> \param[in,out] VF
81*> \verbatim
82*> VF is REAL array, dimension ( K )
83*> On entry, VF contains information passed through DBEDE8.
84*> On exit, VF contains the first K components of the first
85*> components of all right singular vectors of the bidiagonal
86*> matrix.
87*> \endverbatim
88*>
89*> \param[in,out] VL
90*> \verbatim
91*> VL is REAL array, dimension ( K )
92*> On entry, VL contains information passed through DBEDE8.
93*> On exit, VL contains the first K components of the last
94*> components of all right singular vectors of the bidiagonal
95*> matrix.
96*> \endverbatim
97*>
98*> \param[out] DIFL
99*> \verbatim
100*> DIFL is REAL array, dimension ( K )
101*> On exit, DIFL(I) = D(I) - DSIGMA(I).
102*> \endverbatim
103*>
104*> \param[out] DIFR
105*> \verbatim
106*> DIFR is REAL array,
107*> dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
108*> dimension ( K ) if ICOMPQ = 0.
109*> On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
110*> defined and will not be referenced.
111*>
112*> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
113*> normalizing factors for the right singular vector matrix.
114*> \endverbatim
115*>
116*> \param[in] LDDIFR
117*> \verbatim
118*> LDDIFR is INTEGER
119*> The leading dimension of DIFR, must be at least K.
120*> \endverbatim
121*>
122*> \param[in] DSIGMA
123*> \verbatim
124*> DSIGMA is REAL array, dimension ( K )
125*> On entry, the first K elements of this array contain the old
126*> roots of the deflated updating problem. These are the poles
127*> of the secular equation.
128*> \endverbatim
129*>
130*> \param[out] WORK
131*> \verbatim
132*> WORK is REAL array, dimension (3*K)
133*> \endverbatim
134*>
135*> \param[out] INFO
136*> \verbatim
137*> INFO is INTEGER
138*> = 0: successful exit.
139*> < 0: if INFO = -i, the i-th argument had an illegal value.
140*> > 0: if INFO = 1, a singular value did not converge
141*> \endverbatim
142*
143* Authors:
144* ========
145*
146*> \author Univ. of Tennessee
147*> \author Univ. of California Berkeley
148*> \author Univ. of Colorado Denver
149*> \author NAG Ltd.
150*
151*> \ingroup lasd8
152*
153*> \par Contributors:
154* ==================
155*>
156*> Ming Gu and Huan Ren, Computer Science Division, University of
157*> California at Berkeley, USA
158*>
159* =====================================================================
160 SUBROUTINE slasd8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
161 $ DSIGMA, WORK, INFO )
162*
163* -- LAPACK auxiliary routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments ..
168 INTEGER ICOMPQ, INFO, K, LDDIFR
169* ..
170* .. Array Arguments ..
171 REAL D( * ), DIFL( * ), DIFR( LDDIFR, * ),
172 $ dsigma( * ), vf( * ), vl( * ), work( * ),
173 $ z( * )
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 REAL ONE
180 parameter( one = 1.0e+0 )
181* ..
182* .. Local Scalars ..
183 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
184 REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
185* ..
186* .. External Subroutines ..
187 EXTERNAL scopy, slascl, slasd4, slaset,
188 $ xerbla
189* ..
190* .. External Functions ..
191 REAL SDOT, SLAMC3, SNRM2
192 EXTERNAL sdot, slamc3, snrm2
193* ..
194* .. Intrinsic Functions ..
195 INTRINSIC abs, sign, sqrt
196* ..
197* .. Executable Statements ..
198*
199* Test the input parameters.
200*
201 info = 0
202*
203 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
204 info = -1
205 ELSE IF( k.LT.1 ) THEN
206 info = -2
207 ELSE IF( lddifr.LT.k ) THEN
208 info = -9
209 END IF
210 IF( info.NE.0 ) THEN
211 CALL xerbla( 'SLASD8', -info )
212 RETURN
213 END IF
214*
215* Quick return if possible
216*
217 IF( k.EQ.1 ) THEN
218 d( 1 ) = abs( z( 1 ) )
219 difl( 1 ) = d( 1 )
220 IF( icompq.EQ.1 ) THEN
221 difl( 2 ) = one
222 difr( 1, 2 ) = one
223 END IF
224 RETURN
225 END IF
226*
227* Book keeping.
228*
229 iwk1 = 1
230 iwk2 = iwk1 + k
231 iwk3 = iwk2 + k
232 iwk2i = iwk2 - 1
233 iwk3i = iwk3 - 1
234*
235* Normalize Z.
236*
237 rho = snrm2( k, z, 1 )
238 CALL slascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
239 rho = rho*rho
240*
241* Initialize WORK(IWK3).
242*
243 CALL slaset( 'A', k, 1, one, one, work( iwk3 ), k )
244*
245* Compute the updated singular values, the arrays DIFL, DIFR,
246* and the updated Z.
247*
248 DO 40 j = 1, k
249 CALL slasd4( k, j, dsigma, z, work( iwk1 ), rho, d( j ),
250 $ work( iwk2 ), info )
251*
252* If the root finder fails, report the convergence failure.
253*
254 IF( info.NE.0 ) THEN
255 RETURN
256 END IF
257 work( iwk3i+j ) = work( iwk3i+j )*work( j )*work( iwk2i+j )
258 difl( j ) = -work( j )
259 difr( j, 1 ) = -work( j+1 )
260 DO 20 i = 1, j - 1
261 work( iwk3i+i ) = work( iwk3i+i )*work( i )*
262 $ work( iwk2i+i ) / ( dsigma( i )-
263 $ dsigma( j ) ) / ( dsigma( i )+
264 $ dsigma( j ) )
265 20 CONTINUE
266 DO 30 i = j + 1, k
267 work( iwk3i+i ) = work( iwk3i+i )*work( i )*
268 $ work( iwk2i+i ) / ( dsigma( i )-
269 $ dsigma( j ) ) / ( dsigma( i )+
270 $ dsigma( j ) )
271 30 CONTINUE
272 40 CONTINUE
273*
274* Compute updated Z.
275*
276 DO 50 i = 1, k
277 z( i ) = sign( sqrt( abs( work( iwk3i+i ) ) ), z( i ) )
278 50 CONTINUE
279*
280* Update VF and VL.
281*
282 DO 80 j = 1, k
283 diflj = difl( j )
284 dj = d( j )
285 dsigj = -dsigma( j )
286 IF( j.LT.k ) THEN
287 difrj = -difr( j, 1 )
288 dsigjp = -dsigma( j+1 )
289 END IF
290 work( j ) = -z( j ) / diflj / ( dsigma( j )+dj )
291*
292* Use calls to the subroutine SLAMC3 to enforce the parentheses
293* (x+y)+z. The goal is to prevent optimizing compilers
294* from doing x+(y+z).
295*
296 DO 60 i = 1, j - 1
297 work( i ) = z( i ) / ( slamc3( dsigma( i ),
298 $ dsigj )-diflj )
299 $ / ( dsigma( i )+dj )
300 60 CONTINUE
301 DO 70 i = j + 1, k
302 work( i ) = z( i ) / ( slamc3( dsigma( i ),
303 $ dsigjp )+difrj )
304 $ / ( dsigma( i )+dj )
305 70 CONTINUE
306 temp = snrm2( k, work, 1 )
307 work( iwk2i+j ) = sdot( k, work, 1, vf, 1 ) / temp
308 work( iwk3i+j ) = sdot( k, work, 1, vl, 1 ) / temp
309 IF( icompq.EQ.1 ) THEN
310 difr( j, 2 ) = temp
311 END IF
312 80 CONTINUE
313*
314 CALL scopy( k, work( iwk2 ), 1, vf, 1 )
315 CALL scopy( k, work( iwk3 ), 1, vl, 1 )
316*
317 RETURN
318*
319* End of SLASD8
320*
321 END
322
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
subroutine slasd4(n, i, d, z, delta, rho, sigma, work, info)
SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...
Definition slasd4.f:151
subroutine slasd8(icompq, k, d, z, vf, vl, difl, difr, lddifr, dsigma, work, info)
SLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D...
Definition slasd8.f:162
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108