LAPACK 3.12.0
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*> \htmlonly
9*> Download DLASD3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
22* LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
23* INFO )
24*
25* .. Scalar Arguments ..
26* INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
27* $ SQRE
28* ..
29* .. Array Arguments ..
30* INTEGER CTOT( * ), IDXC( * )
31* DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
32* $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ),
33* $ Z( * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> DLASD3 finds all the square roots of the roots of the secular
43*> equation, as defined by the values in D and Z. It makes the
44*> appropriate calls to DLASD4 and then updates the singular
45*> vectors by matrix multiplication.
46*>
47*> DLASD3 is called from DLASD1.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] NL
54*> \verbatim
55*> NL is INTEGER
56*> The row dimension of the upper block. NL >= 1.
57*> \endverbatim
58*>
59*> \param[in] NR
60*> \verbatim
61*> NR is INTEGER
62*> The row dimension of the lower block. NR >= 1.
63*> \endverbatim
64*>
65*> \param[in] SQRE
66*> \verbatim
67*> SQRE is INTEGER
68*> = 0: the lower block is an NR-by-NR square matrix.
69*> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
70*>
71*> The bidiagonal matrix has N = NL + NR + 1 rows and
72*> M = N + SQRE >= N columns.
73*> \endverbatim
74*>
75*> \param[in] K
76*> \verbatim
77*> K is INTEGER
78*> The size of the secular equation, 1 =< K = < N.
79*> \endverbatim
80*>
81*> \param[out] D
82*> \verbatim
83*> D is DOUBLE PRECISION array, dimension(K)
84*> On exit the square roots of the roots of the secular equation,
85*> in ascending order.
86*> \endverbatim
87*>
88*> \param[out] Q
89*> \verbatim
90*> Q is DOUBLE PRECISION array, dimension (LDQ,K)
91*> \endverbatim
92*>
93*> \param[in] LDQ
94*> \verbatim
95*> LDQ is INTEGER
96*> The leading dimension of the array Q. LDQ >= K.
97*> \endverbatim
98*>
99*> \param[in] DSIGMA
100*> \verbatim
101*> DSIGMA is DOUBLE PRECISION array, dimension(K)
102*> The first K elements of this array contain the old roots
103*> of the deflated updating problem. These are the poles
104*> of the secular equation.
105*> \endverbatim
106*>
107*> \param[out] U
108*> \verbatim
109*> U is DOUBLE PRECISION array, dimension (LDU, N)
110*> The last N - K columns of this matrix contain the deflated
111*> left singular vectors.
112*> \endverbatim
113*>
114*> \param[in] LDU
115*> \verbatim
116*> LDU is INTEGER
117*> The leading dimension of the array U. LDU >= N.
118*> \endverbatim
119*>
120*> \param[in] U2
121*> \verbatim
122*> U2 is DOUBLE PRECISION array, dimension (LDU2, N)
123*> The first K columns of this matrix contain the non-deflated
124*> left singular vectors for the split problem.
125*> \endverbatim
126*>
127*> \param[in] LDU2
128*> \verbatim
129*> LDU2 is INTEGER
130*> The leading dimension of the array U2. LDU2 >= N.
131*> \endverbatim
132*>
133*> \param[out] VT
134*> \verbatim
135*> VT is DOUBLE PRECISION array, dimension (LDVT, M)
136*> The last M - K columns of VT**T contain the deflated
137*> right singular vectors.
138*> \endverbatim
139*>
140*> \param[in] LDVT
141*> \verbatim
142*> LDVT is INTEGER
143*> The leading dimension of the array VT. LDVT >= N.
144*> \endverbatim
145*>
146*> \param[in,out] VT2
147*> \verbatim
148*> VT2 is DOUBLE PRECISION array, dimension (LDVT2, N)
149*> The first K columns of VT2**T contain the non-deflated
150*> right singular vectors for the split problem.
151*> \endverbatim
152*>
153*> \param[in] LDVT2
154*> \verbatim
155*> LDVT2 is INTEGER
156*> The leading dimension of the array VT2. LDVT2 >= N.
157*> \endverbatim
158*>
159*> \param[in] IDXC
160*> \verbatim
161*> IDXC is INTEGER array, dimension ( N )
162*> The permutation used to arrange the columns of U (and rows of
163*> VT) into three groups: the first group contains non-zero
164*> entries only at and above (or before) NL +1; the second
165*> contains non-zero entries only at and below (or after) NL+2;
166*> and the third is dense. The first column of U and the row of
167*> VT are treated separately, however.
168*>
169*> The rows of the singular vectors found by DLASD4
170*> must be likewise permuted before the matrix multiplies can
171*> take place.
172*> \endverbatim
173*>
174*> \param[in] CTOT
175*> \verbatim
176*> CTOT is INTEGER array, dimension ( 4 )
177*> A count of the total number of the various types of columns
178*> in U (or rows in VT), as described in IDXC. The fourth column
179*> type is any column which has been deflated.
180*> \endverbatim
181*>
182*> \param[in,out] Z
183*> \verbatim
184*> Z is DOUBLE PRECISION array, dimension (K)
185*> The first K elements of this array contain the components
186*> of the deflation-adjusted updating row vector.
187*> \endverbatim
188*>
189*> \param[out] INFO
190*> \verbatim
191*> INFO is INTEGER
192*> = 0: successful exit.
193*> < 0: if INFO = -i, the i-th argument had an illegal value.
194*> > 0: if INFO = 1, a singular value did not converge
195*> \endverbatim
196*
197* Authors:
198* ========
199*
200*> \author Univ. of Tennessee
201*> \author Univ. of California Berkeley
202*> \author Univ. of Colorado Denver
203*> \author NAG Ltd.
204*
205*> \ingroup lasd3
206*
207*> \par Contributors:
208* ==================
209*>
210*> Ming Gu and Huan Ren, Computer Science Division, University of
211*> California at Berkeley, USA
212*>
213* =====================================================================
214 SUBROUTINE dlasd3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2,
215 $ LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z,
216 $ INFO )
217*
218* -- LAPACK auxiliary routine --
219* -- LAPACK is a software package provided by Univ. of Tennessee, --
220* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221*
222* .. Scalar Arguments ..
223 INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR,
224 $ SQRE
225* ..
226* .. Array Arguments ..
227 INTEGER CTOT( * ), IDXC( * )
228 DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ),
229 $ u2( ldu2, * ), vt( ldvt, * ), vt2( ldvt2, * ),
230 $ z( * )
231* ..
232*
233* =====================================================================
234*
235* .. Parameters ..
236 DOUBLE PRECISION ONE, ZERO, NEGONE
237 PARAMETER ( ONE = 1.0d+0, zero = 0.0d+0,
238 $ negone = -1.0d+0 )
239* ..
240* .. Local Scalars ..
241 INTEGER CTEMP, I, J, JC, KTEMP, M, N, NLP1, NLP2, NRP1
242 DOUBLE PRECISION RHO, TEMP
243* ..
244* .. External Functions ..
245 DOUBLE PRECISION DNRM2
246 EXTERNAL DNRM2
247* ..
248* .. External Subroutines ..
249 EXTERNAL dcopy, dgemm, dlacpy, dlascl, dlasd4, 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, u,
368 $ ldu )
369 GO TO 100
370 END IF
371 IF( ctot( 1 ).GT.0 ) THEN
372 CALL dgemm( 'N', 'N', nl, k, ctot( 1 ), one, u2( 1, 2 ), ldu2,
373 $ q( 2, 1 ), ldq, zero, u( 1, 1 ), ldu )
374 IF( ctot( 3 ).GT.0 ) THEN
375 ktemp = 2 + ctot( 1 ) + ctot( 2 )
376 CALL dgemm( 'N', 'N', nl, k, ctot( 3 ), one, u2( 1, ktemp ),
377 $ ldu2, q( ktemp, 1 ), ldq, one, u( 1, 1 ), ldu )
378 END IF
379 ELSE IF( ctot( 3 ).GT.0 ) THEN
380 ktemp = 2 + ctot( 1 ) + ctot( 2 )
381 CALL dgemm( 'N', 'N', nl, k, ctot( 3 ), one, u2( 1, ktemp ),
382 $ ldu2, q( ktemp, 1 ), ldq, zero, u( 1, 1 ), ldu )
383 ELSE
384 CALL dlacpy( 'F', nl, k, u2, ldu2, u, ldu )
385 END IF
386 CALL dcopy( k, q( 1, 1 ), ldq, u( nlp1, 1 ), ldu )
387 ktemp = 2 + ctot( 1 )
388 ctemp = ctot( 2 ) + ctot( 3 )
389 CALL dgemm( 'N', 'N', nr, k, ctemp, one, u2( nlp2, ktemp ), ldu2,
390 $ q( ktemp, 1 ), ldq, zero, u( nlp2, 1 ), ldu )
391*
392* Generate the right singular vectors.
393*
394 100 CONTINUE
395 DO 120 i = 1, k
396 temp = dnrm2( k, vt( 1, i ), 1 )
397 q( i, 1 ) = vt( 1, i ) / temp
398 DO 110 j = 2, k
399 jc = idxc( j )
400 q( i, j ) = vt( jc, i ) / temp
401 110 CONTINUE
402 120 CONTINUE
403*
404* Update the right singular vector matrix.
405*
406 IF( k.EQ.2 ) THEN
407 CALL dgemm( 'N', 'N', k, m, k, one, q, ldq, vt2, ldvt2, zero,
408 $ vt, ldvt )
409 RETURN
410 END IF
411 ktemp = 1 + ctot( 1 )
412 CALL dgemm( 'N', 'N', k, nlp1, ktemp, one, q( 1, 1 ), ldq,
413 $ vt2( 1, 1 ), ldvt2, zero, vt( 1, 1 ), ldvt )
414 ktemp = 2 + ctot( 1 ) + ctot( 2 )
415 IF( ktemp.LE.ldvt2 )
416 $ CALL dgemm( 'N', 'N', k, nlp1, ctot( 3 ), one, q( 1, ktemp ),
417 $ ldq, vt2( ktemp, 1 ), ldvt2, one, vt( 1, 1 ),
418 $ ldvt )
419*
420 ktemp = ctot( 1 ) + 1
421 nrp1 = nr + sqre
422 IF( ktemp.GT.1 ) THEN
423 DO 130 i = 1, k
424 q( i, ktemp ) = q( i, 1 )
425 130 CONTINUE
426 DO 140 i = nlp2, m
427 vt2( ktemp, i ) = vt2( 1, i )
428 140 CONTINUE
429 END IF
430 ctemp = 1 + ctot( 2 ) + ctot( 3 )
431 CALL dgemm( 'N', 'N', k, nrp1, ctemp, one, q( 1, ktemp ), ldq,
432 $ vt2( ktemp, nlp2 ), ldvt2, zero, vt( 1, nlp2 ), ldvt )
433*
434 RETURN
435*
436* End of DLASD3
437*
438 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:103
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:143
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:217
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:153