LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlasd3.f
Go to the documentation of this file.
1*> \brief \b DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and Z, and then updates the singular vectors by matrix multiplication. Used by sbdsdc.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DLASD3 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd3.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd3.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd3.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
20* LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
21* INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
25* $ SQRE
26* ..
27* .. Array Arguments ..
28* INTEGER CTOT( * ), IDXC( * )
29* DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
30* $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
31* $ Z( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> DLASD3 finds all the square roots of the roots of the secular
41*> equation, as defined by the values in D and Z. It makes the
42*> appropriate calls to DLASD4 and then updates the singular
43*> vectors by matrix multiplication.
44*>
45*> DLASD3 is called from DLASD1.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] NL
52*> \verbatim
53*> NL is INTEGER
54*> The row dimension of the upper block. NL >= 1.
55*> \endverbatim
56*>
57*> \param[in] NR
58*> \verbatim
59*> NR is INTEGER
60*> The row dimension of the lower block. NR >= 1.
61*> \endverbatim
62*>
63*> \param[in] SQRE
64*> \verbatim
65*> SQRE is INTEGER
66*> = 0: the lower block is an NR-by-NR square matrix.
67*> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
68*>
69*> The bidiagonal matrix has N = NL + NR + 1 rows and
70*> M = N + SQRE >= N columns.
71*> \endverbatim
72*>
73*> \param[in] K
74*> \verbatim
75*> K is INTEGER
76*> The size of the secular equation, 1 =< K = < N.
77*> \endverbatim
78*>
79*> \param[out] D
80*> \verbatim
81*> D is DOUBLE PRECISION array, dimension(K)
82*> On exit the square roots of the roots of the secular equation,
83*> in ascending order.
84*> \endverbatim
85*>
86*> \param[out] Q
87*> \verbatim
88*> Q is DOUBLE PRECISION array, dimension (LDQ,K)
89*> \endverbatim
90*>
91*> \param[in] LDQ
92*> \verbatim
93*> LDQ is INTEGER
94*> The leading dimension of the array Q. LDQ >= K.
95*> \endverbatim
96*>
97*> \param[in] DSIGMA
98*> \verbatim
99*> DSIGMA is DOUBLE PRECISION array, dimension(K)
100*> The first K elements of this array contain the old roots
101*> of the deflated updating problem. These are the poles
102*> of the secular equation.
103*> \endverbatim
104*>
105*> \param[out] U
106*> \verbatim
107*> U is DOUBLE PRECISION array, dimension (LDU, N)
108*> The last N - K columns of this matrix contain the deflated
109*> left singular vectors.
110*> \endverbatim
111*>
112*> \param[in] LDU
113*> \verbatim
114*> LDU is INTEGER
115*> The leading dimension of the array U. LDU >= N.
116*> \endverbatim
117*>
118*> \param[in] U2
119*> \verbatim
120*> U2 is DOUBLE PRECISION array, dimension (LDU2, N)
121*> The first K columns of this matrix contain the non-deflated
122*> left singular vectors for the split problem.
123*> \endverbatim
124*>
125*> \param[in] LDU2
126*> \verbatim
127*> LDU2 is INTEGER
128*> The leading dimension of the array U2. LDU2 >= N.
129*> \endverbatim
130*>
131*> \param[out] VT
132*> \verbatim
133*> VT is DOUBLE PRECISION array, dimension (LDVT, M)
134*> The last M - K columns of VT**T contain the deflated
135*> right singular vectors.
136*> \endverbatim
137*>
138*> \param[in] LDVT
139*> \verbatim
140*> LDVT is INTEGER
141*> The leading dimension of the array VT. LDVT >= N.
142*> \endverbatim
143*>
144*> \param[in,out] VT2
145*> \verbatim
146*> VT2 is DOUBLE PRECISION array, dimension (LDVT2, N)
147*> The first K columns of VT2**T contain the non-deflated
148*> right singular vectors for the split problem.
149*> \endverbatim
150*>
151*> \param[in] LDVT2
152*> \verbatim
153*> LDVT2 is INTEGER
154*> The leading dimension of the array VT2. LDVT2 >= N.
155*> \endverbatim
156*>
157*> \param[in] IDXC
158*> \verbatim
159*> IDXC is INTEGER array, dimension ( N )
160*> The permutation used to arrange the columns of U (and rows of
161*> VT) into three groups: the first group contains non-zero
162*> entries only at and above (or before) NL +1; the second
163*> contains non-zero entries only at and below (or after) NL+2;
164*> and the third is dense. The first column of U and the row of
165*> VT are treated separately, however.
166*>
167*> The rows of the singular vectors found by DLASD4
168*> must be likewise permuted before the matrix multiplies can
169*> take place.
170*> \endverbatim
171*>
172*> \param[in] CTOT
173*> \verbatim
174*> CTOT is INTEGER array, dimension ( 4 )
175*> A count of the total number of the various types of columns
176*> in U (or rows in VT), as described in IDXC. The fourth column
177*> type is any column which has been deflated.
178*> \endverbatim
179*>
180*> \param[in,out] Z
181*> \verbatim
182*> Z is DOUBLE PRECISION array, dimension (K)
183*> The first K elements of this array contain the components
184*> of the deflation-adjusted updating row vector.
185*> \endverbatim
186*>
187*> \param[out] INFO
188*> \verbatim
189*> INFO is INTEGER
190*> = 0: successful exit.
191*> < 0: if INFO = -i, the i-th argument had an illegal value.
192*> > 0: if INFO = 1, a singular value did not converge
193*> \endverbatim
194*
195* Authors:
196* ========
197*
198*> \author Univ. of Tennessee
199*> \author Univ. of California Berkeley
200*> \author Univ. of Colorado Denver
201*> \author NAG Ltd.
202*
203*> \ingroup lasd3
204*
205*> \par Contributors:
206* ==================
207*>
208*> Ming Gu and Huan Ren, Computer Science Division, University of
209*> California at Berkeley, USA
210*>
211* =====================================================================
212 SUBROUTINE dlasd3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU,
213 $ U2,
214 $ LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
215 $ INFO )
216*
217* -- LAPACK auxiliary routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
223 $ SQRE
224* ..
225* .. Array Arguments ..
226 INTEGER CTOT( * ), IDXC( * )
227 DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
228 $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
229 $ z( * )
230* ..
231*
232* =====================================================================
233*
234* .. Parameters ..
235 DOUBLE PRECISION ONE, ZERO, NEGONE
236 PARAMETER ( ONE = 1.0d+0, zero = 0.0d+0,
237 $ negone = -1.0d+0 )
238* ..
239* .. Local Scalars ..
240 INTEGER CTEMP, I, J, JC, KTEMP, M, N, NLP1, NLP2, NRP1
241 DOUBLE PRECISION RHO, TEMP
242* ..
243* .. External Functions ..
244 DOUBLE PRECISION DNRM2
245 EXTERNAL DNRM2
246* ..
247* .. External Subroutines ..
248 EXTERNAL dcopy, dgemm, dlacpy, dlascl, dlasd4,
249 $ xerbla
250* ..
251* .. Intrinsic Functions ..
252 INTRINSIC abs, sign, sqrt
253* ..
254* .. Executable Statements ..
255*
256* Test the input parameters.
257*
258 info = 0
259*
260 IF( nl.LT.1 ) THEN
261 info = -1
262 ELSE IF( nr.LT.1 ) THEN
263 info = -2
264 ELSE IF( ( sqre.NE.1 ) .AND. ( sqre.NE.0 ) ) THEN
265 info = -3
266 END IF
267*
268 n = nl + nr + 1
269 m = n + sqre
270 nlp1 = nl + 1
271 nlp2 = nl + 2
272*
273 IF( ( k.LT.1 ) .OR. ( k.GT.n ) ) THEN
274 info = -4
275 ELSE IF( ldq.LT.k ) THEN
276 info = -7
277 ELSE IF( ldu.LT.n ) THEN
278 info = -10
279 ELSE IF( ldu2.LT.n ) THEN
280 info = -12
281 ELSE IF( ldvt.LT.m ) THEN
282 info = -14
283 ELSE IF( ldvt2.LT.m ) THEN
284 info = -16
285 END IF
286 IF( info.NE.0 ) THEN
287 CALL xerbla( 'DLASD3', -info )
288 RETURN
289 END IF
290*
291* Quick return if possible
292*
293 IF( k.EQ.1 ) THEN
294 d( 1 ) = abs( z( 1 ) )
295 CALL dcopy( m, vt2( 1, 1 ), ldvt2, vt( 1, 1 ), ldvt )
296 IF( z( 1 ).GT.zero ) THEN
297 CALL dcopy( n, u2( 1, 1 ), 1, u( 1, 1 ), 1 )
298 ELSE
299 DO 10 i = 1, n
300 u( i, 1 ) = -u2( i, 1 )
301 10 CONTINUE
302 END IF
303 RETURN
304 END IF
305*
306* Keep a copy of Z.
307*
308 CALL dcopy( k, z, 1, q, 1 )
309*
310* Normalize Z.
311*
312 rho = dnrm2( k, z, 1 )
313 CALL dlascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
314 rho = rho*rho
315*
316* Find the new singular values.
317*
318 DO 30 j = 1, k
319 CALL dlasd4( k, j, dsigma, z, u( 1, j ), rho, d( j ),
320 $ vt( 1, j ), info )
321*
322* If the zero finder fails, report the convergence failure.
323*
324 IF( info.NE.0 ) THEN
325 RETURN
326 END IF
327 30 CONTINUE
328*
329* Compute updated Z.
330*
331 DO 60 i = 1, k
332 z( i ) = u( i, k )*vt( i, k )
333 DO 40 j = 1, i - 1
334 z( i ) = z( i )*( u( i, j )*vt( i, j ) /
335 $ ( dsigma( i )-dsigma( j ) ) /
336 $ ( dsigma( i )+dsigma( j ) ) )
337 40 CONTINUE
338 DO 50 j = i, k - 1
339 z( i ) = z( i )*( u( i, j )*vt( i, j ) /
340 $ ( dsigma( i )-dsigma( j+1 ) ) /
341 $ ( dsigma( i )+dsigma( j+1 ) ) )
342 50 CONTINUE
343 z( i ) = sign( sqrt( abs( z( i ) ) ), q( i, 1 ) )
344 60 CONTINUE
345*
346* Compute left singular vectors of the modified diagonal matrix,
347* and store related information for the right singular vectors.
348*
349 DO 90 i = 1, k
350 vt( 1, i ) = z( 1 ) / u( 1, i ) / vt( 1, i )
351 u( 1, i ) = negone
352 DO 70 j = 2, k
353 vt( j, i ) = z( j ) / u( j, i ) / vt( j, i )
354 u( j, i ) = dsigma( j )*vt( j, i )
355 70 CONTINUE
356 temp = dnrm2( k, u( 1, i ), 1 )
357 q( 1, i ) = u( 1, i ) / temp
358 DO 80 j = 2, k
359 jc = idxc( j )
360 q( j, i ) = u( jc, i ) / temp
361 80 CONTINUE
362 90 CONTINUE
363*
364* Update the left singular vector matrix.
365*
366 IF( k.EQ.2 ) THEN
367 CALL dgemm( 'N', 'N', n, k, k, one, u2, ldu2, q, ldq, zero,
368 $ u,
369 $ ldu )
370 GO TO 100
371 END IF
372 IF( ctot( 1 ).GT.0 ) THEN
373 CALL dgemm( 'N', 'N', nl, k, ctot( 1 ), one, u2( 1, 2 ),
374 $ ldu2,
375 $ q( 2, 1 ), ldq, zero, u( 1, 1 ), ldu )
376 IF( ctot( 3 ).GT.0 ) THEN
377 ktemp = 2 + ctot( 1 ) + ctot( 2 )
378 CALL dgemm( 'N', 'N', nl, k, ctot( 3 ), one, u2( 1,
379 $ ktemp ),
380 $ ldu2, q( ktemp, 1 ), ldq, one, u( 1, 1 ), ldu )
381 END IF
382 ELSE IF( ctot( 3 ).GT.0 ) THEN
383 ktemp = 2 + ctot( 1 ) + ctot( 2 )
384 CALL dgemm( 'N', 'N', nl, k, ctot( 3 ), one, u2( 1, ktemp ),
385 $ ldu2, q( ktemp, 1 ), ldq, zero, u( 1, 1 ), ldu )
386 ELSE
387 CALL dlacpy( 'F', nl, k, u2, ldu2, u, ldu )
388 END IF
389 CALL dcopy( k, q( 1, 1 ), ldq, u( nlp1, 1 ), ldu )
390 ktemp = 2 + ctot( 1 )
391 ctemp = ctot( 2 ) + ctot( 3 )
392 CALL dgemm( 'N', 'N', nr, k, ctemp, one, u2( nlp2, ktemp ),
393 $ ldu2,
394 $ q( ktemp, 1 ), ldq, zero, u( nlp2, 1 ), ldu )
395*
396* Generate the right singular vectors.
397*
398 100 CONTINUE
399 DO 120 i = 1, k
400 temp = dnrm2( k, vt( 1, i ), 1 )
401 q( i, 1 ) = vt( 1, i ) / temp
402 DO 110 j = 2, k
403 jc = idxc( j )
404 q( i, j ) = vt( jc, i ) / temp
405 110 CONTINUE
406 120 CONTINUE
407*
408* Update the right singular vector matrix.
409*
410 IF( k.EQ.2 ) THEN
411 CALL dgemm( 'N', 'N', k, m, k, one, q, ldq, vt2, ldvt2,
412 $ zero,
413 $ vt, ldvt )
414 RETURN
415 END IF
416 ktemp = 1 + ctot( 1 )
417 CALL dgemm( 'N', 'N', k, nlp1, ktemp, one, q( 1, 1 ), ldq,
418 $ vt2( 1, 1 ), ldvt2, zero, vt( 1, 1 ), ldvt )
419 ktemp = 2 + ctot( 1 ) + ctot( 2 )
420 IF( ktemp.LE.ldvt2 )
421 $ CALL dgemm( 'N', 'N', k, nlp1, ctot( 3 ), one, q( 1,
422 $ ktemp ),
423 $ ldq, vt2( ktemp, 1 ), ldvt2, one, vt( 1, 1 ),
424 $ ldvt )
425*
426 ktemp = ctot( 1 ) + 1
427 nrp1 = nr + sqre
428 IF( ktemp.GT.1 ) THEN
429 DO 130 i = 1, k
430 q( i, ktemp ) = q( i, 1 )
431 130 CONTINUE
432 DO 140 i = nlp2, m
433 vt2( ktemp, i ) = vt2( 1, i )
434 140 CONTINUE
435 END IF
436 ctemp = 1 + ctot( 2 ) + ctot( 3 )
437 CALL dgemm( 'N', 'N', k, nrp1, ctemp, one, q( 1, ktemp ), ldq,
438 $ vt2( ktemp, nlp2 ), ldvt2, zero, vt( 1, nlp2 ), ldvt )
439*
440 RETURN
441*
442* End of DLASD3
443*
444 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlascl
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:142
dlasd3
subroutine dlasd3(nl, nr, sqre, k, d, q, ldq, dsigma, u, ldu, u2, ldu2, vt, ldvt, vt2, ldvt2, idxc, ctot, z, info)
DLASD3 finds all square roots of the roots of the secular equation, as defined by the values in D and...
Definition dlasd3.f:216
dlasd4
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:151