LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dbdsdc.f
Go to the documentation of this file.
1*> \brief \b DBDSDC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DBDSDC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
22* WORK, IWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER COMPQ, UPLO
26* INTEGER INFO, LDU, LDVT, N
27* ..
28* .. Array Arguments ..
29* INTEGER IQ( * ), IWORK( * )
30* DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
31* $ VT( LDVT, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> DBDSDC computes the singular value decomposition (SVD) of a real
41*> N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
42*> using a divide and conquer method, where S is a diagonal matrix
43*> with non-negative diagonal elements (the singular values of B), and
44*> U and VT are orthogonal matrices of left and right singular vectors,
45*> respectively. DBDSDC can be used to compute all singular values,
46*> and optionally, singular vectors or singular vectors in compact form.
47*>
48*> The code currently calls DLASDQ if singular values only are desired.
49*> However, it can be slightly modified to compute singular values
50*> using the divide and conquer method.
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] UPLO
57*> \verbatim
58*> UPLO is CHARACTER*1
59*> = 'U': B is upper bidiagonal.
60*> = 'L': B is lower bidiagonal.
61*> \endverbatim
62*>
63*> \param[in] COMPQ
64*> \verbatim
65*> COMPQ is CHARACTER*1
66*> Specifies whether singular vectors are to be computed
67*> as follows:
68*> = 'N': Compute singular values only;
69*> = 'P': Compute singular values and compute singular
70*> vectors in compact form;
71*> = 'I': Compute singular values and singular vectors.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*> N is INTEGER
77*> The order of the matrix B. N >= 0.
78*> \endverbatim
79*>
80*> \param[in,out] D
81*> \verbatim
82*> D is DOUBLE PRECISION array, dimension (N)
83*> On entry, the n diagonal elements of the bidiagonal matrix B.
84*> On exit, if INFO=0, the singular values of B.
85*> \endverbatim
86*>
87*> \param[in,out] E
88*> \verbatim
89*> E is DOUBLE PRECISION array, dimension (N-1)
90*> On entry, the elements of E contain the offdiagonal
91*> elements of the bidiagonal matrix whose SVD is desired.
92*> On exit, E has been destroyed.
93*> \endverbatim
94*>
95*> \param[out] U
96*> \verbatim
97*> U is DOUBLE PRECISION array, dimension (LDU,N)
98*> If COMPQ = 'I', then:
99*> On exit, if INFO = 0, U contains the left singular vectors
100*> of the bidiagonal matrix.
101*> For other values of COMPQ, U is not referenced.
102*> \endverbatim
103*>
104*> \param[in] LDU
105*> \verbatim
106*> LDU is INTEGER
107*> The leading dimension of the array U. LDU >= 1.
108*> If singular vectors are desired, then LDU >= max( 1, N ).
109*> \endverbatim
110*>
111*> \param[out] VT
112*> \verbatim
113*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
114*> If COMPQ = 'I', then:
115*> On exit, if INFO = 0, VT**T contains the right singular
116*> vectors of the bidiagonal matrix.
117*> For other values of COMPQ, VT is not referenced.
118*> \endverbatim
119*>
120*> \param[in] LDVT
121*> \verbatim
122*> LDVT is INTEGER
123*> The leading dimension of the array VT. LDVT >= 1.
124*> If singular vectors are desired, then LDVT >= max( 1, N ).
125*> \endverbatim
126*>
127*> \param[out] Q
128*> \verbatim
129*> Q is DOUBLE PRECISION array, dimension (LDQ)
130*> If COMPQ = 'P', then:
131*> On exit, if INFO = 0, Q and IQ contain the left
132*> and right singular vectors in a compact form,
133*> requiring O(N log N) space instead of 2*N**2.
134*> In particular, Q contains all the DOUBLE PRECISION data in
135*> LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
136*> words of memory, where SMLSIZ is returned by ILAENV and
137*> is equal to the maximum size of the subproblems at the
138*> bottom of the computation tree (usually about 25).
139*> For other values of COMPQ, Q is not referenced.
140*> \endverbatim
141*>
142*> \param[out] IQ
143*> \verbatim
144*> IQ is INTEGER array, dimension (LDIQ)
145*> If COMPQ = 'P', then:
146*> On exit, if INFO = 0, Q and IQ contain the left
147*> and right singular vectors in a compact form,
148*> requiring O(N log N) space instead of 2*N**2.
149*> In particular, IQ contains all INTEGER data in
150*> LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
151*> words of memory, where SMLSIZ is returned by ILAENV and
152*> is equal to the maximum size of the subproblems at the
153*> bottom of the computation tree (usually about 25).
154*> For other values of COMPQ, IQ is not referenced.
155*> \endverbatim
156*>
157*> \param[out] WORK
158*> \verbatim
159*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
160*> If COMPQ = 'N' then LWORK >= (4 * N).
161*> If COMPQ = 'P' then LWORK >= (6 * N).
162*> If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
163*> \endverbatim
164*>
165*> \param[out] IWORK
166*> \verbatim
167*> IWORK is INTEGER array, dimension (8*N)
168*> \endverbatim
169*>
170*> \param[out] INFO
171*> \verbatim
172*> INFO is INTEGER
173*> = 0: successful exit.
174*> < 0: if INFO = -i, the i-th argument had an illegal value.
175*> > 0: The algorithm failed to compute a singular value.
176*> The update process of divide and conquer failed.
177*> \endverbatim
178*
179* Authors:
180* ========
181*
182*> \author Univ. of Tennessee
183*> \author Univ. of California Berkeley
184*> \author Univ. of Colorado Denver
185*> \author NAG Ltd.
186*
187*> \ingroup bdsdc
188*
189*> \par Contributors:
190* ==================
191*>
192*> Ming Gu and Huan Ren, Computer Science Division, University of
193*> California at Berkeley, USA
194*>
195* =====================================================================
196 SUBROUTINE dbdsdc( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
197 $ WORK, IWORK, INFO )
198*
199* -- LAPACK computational routine --
200* -- LAPACK is a software package provided by Univ. of Tennessee, --
201* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
202*
203* .. Scalar Arguments ..
204 CHARACTER COMPQ, UPLO
205 INTEGER INFO, LDU, LDVT, N
206* ..
207* .. Array Arguments ..
208 INTEGER IQ( * ), IWORK( * )
209 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
210 $ vt( ldvt, * ), work( * )
211* ..
212*
213* =====================================================================
214* Changed dimension statement in comment describing E from (N) to
215* (N-1). Sven, 17 Feb 05.
216* =====================================================================
217*
218* .. Parameters ..
219 DOUBLE PRECISION ZERO, ONE, TWO
220 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
221* ..
222* .. Local Scalars ..
223 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
224 $ icompq, ierr, ii, is, iu, iuplo, ivt, j, k, kk,
225 $ mlvl, nm1, nsize, perm, poles, qstart, smlsiz,
226 $ smlszp, sqre, start, wstart, z
227 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
228* ..
229* .. External Functions ..
230 LOGICAL LSAME
231 INTEGER ILAENV
232 DOUBLE PRECISION DLAMCH, DLANST
233 EXTERNAL lsame, ilaenv, dlamch, dlanst
234* ..
235* .. External Subroutines ..
236 EXTERNAL dcopy, dlartg, dlascl, dlasd0, dlasda, dlasdq,
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC abs, dble, int, log, sign
241* ..
242* .. Executable Statements ..
243*
244* Test the input parameters.
245*
246 info = 0
247*
248 iuplo = 0
249 IF( lsame( uplo, 'U' ) )
250 $ iuplo = 1
251 IF( lsame( uplo, 'L' ) )
252 $ iuplo = 2
253 IF( lsame( compq, 'N' ) ) THEN
254 icompq = 0
255 ELSE IF( lsame( compq, 'P' ) ) THEN
256 icompq = 1
257 ELSE IF( lsame( compq, 'I' ) ) THEN
258 icompq = 2
259 ELSE
260 icompq = -1
261 END IF
262 IF( iuplo.EQ.0 ) THEN
263 info = -1
264 ELSE IF( icompq.LT.0 ) THEN
265 info = -2
266 ELSE IF( n.LT.0 ) THEN
267 info = -3
268 ELSE IF( ( ldu.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldu.LT.
269 $ n ) ) ) THEN
270 info = -7
271 ELSE IF( ( ldvt.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldvt.LT.
272 $ n ) ) ) THEN
273 info = -9
274 END IF
275 IF( info.NE.0 ) THEN
276 CALL xerbla( 'DBDSDC', -info )
277 RETURN
278 END IF
279*
280* Quick return if possible
281*
282 IF( n.EQ.0 )
283 $ RETURN
284 smlsiz = ilaenv( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
285 IF( n.EQ.1 ) THEN
286 IF( icompq.EQ.1 ) THEN
287 q( 1 ) = sign( one, d( 1 ) )
288 q( 1+smlsiz*n ) = one
289 ELSE IF( icompq.EQ.2 ) THEN
290 u( 1, 1 ) = sign( one, d( 1 ) )
291 vt( 1, 1 ) = one
292 END IF
293 d( 1 ) = abs( d( 1 ) )
294 RETURN
295 END IF
296 nm1 = n - 1
297*
298* If matrix lower bidiagonal, rotate to be upper bidiagonal
299* by applying Givens rotations on the left
300*
301 wstart = 1
302 qstart = 3
303 IF( icompq.EQ.1 ) THEN
304 CALL dcopy( n, d, 1, q( 1 ), 1 )
305 CALL dcopy( n-1, e, 1, q( n+1 ), 1 )
306 END IF
307 IF( iuplo.EQ.2 ) THEN
308 qstart = 5
309 IF( icompq .EQ. 2 ) wstart = 2*n - 1
310 DO 10 i = 1, n - 1
311 CALL dlartg( d( i ), e( i ), cs, sn, r )
312 d( i ) = r
313 e( i ) = sn*d( i+1 )
314 d( i+1 ) = cs*d( i+1 )
315 IF( icompq.EQ.1 ) THEN
316 q( i+2*n ) = cs
317 q( i+3*n ) = sn
318 ELSE IF( icompq.EQ.2 ) THEN
319 work( i ) = cs
320 work( nm1+i ) = -sn
321 END IF
322 10 CONTINUE
323 END IF
324*
325* If ICOMPQ = 0, use DLASDQ to compute the singular values.
326*
327 IF( icompq.EQ.0 ) THEN
328* Ignore WSTART, instead using WORK( 1 ), since the two vectors
329* for CS and -SN above are added only if ICOMPQ == 2,
330* and adding them exceeds documented WORK size of 4*n.
331 CALL dlasdq( 'U', 0, n, 0, 0, 0, d, e, vt, ldvt, u, ldu, u,
332 $ ldu, work( 1 ), info )
333 GO TO 40
334 END IF
335*
336* If N is smaller than the minimum divide size SMLSIZ, then solve
337* the problem with another solver.
338*
339 IF( n.LE.smlsiz ) THEN
340 IF( icompq.EQ.2 ) THEN
341 CALL dlaset( 'A', n, n, zero, one, u, ldu )
342 CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
343 CALL dlasdq( 'U', 0, n, n, n, 0, d, e, vt, ldvt, u, ldu, u,
344 $ ldu, work( wstart ), info )
345 ELSE IF( icompq.EQ.1 ) THEN
346 iu = 1
347 ivt = iu + n
348 CALL dlaset( 'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
349 $ n )
350 CALL dlaset( 'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
351 $ n )
352 CALL dlasdq( 'U', 0, n, n, n, 0, d, e,
353 $ q( ivt+( qstart-1 )*n ), n,
354 $ q( iu+( qstart-1 )*n ), n,
355 $ q( iu+( qstart-1 )*n ), n, work( wstart ),
356 $ info )
357 END IF
358 GO TO 40
359 END IF
360*
361 IF( icompq.EQ.2 ) THEN
362 CALL dlaset( 'A', n, n, zero, one, u, ldu )
363 CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
364 END IF
365*
366* Scale.
367*
368 orgnrm = dlanst( 'M', n, d, e )
369 IF( orgnrm.EQ.zero )
370 $ RETURN
371 CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
372 CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
373*
374 eps = (0.9d+0)*dlamch( 'Epsilon' )
375*
376 mlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
377 smlszp = smlsiz + 1
378*
379 IF( icompq.EQ.1 ) THEN
380 iu = 1
381 ivt = 1 + smlsiz
382 difl = ivt + smlszp
383 difr = difl + mlvl
384 z = difr + mlvl*2
385 ic = z + mlvl
386 is = ic + 1
387 poles = is + 1
388 givnum = poles + 2*mlvl
389*
390 k = 1
391 givptr = 2
392 perm = 3
393 givcol = perm + mlvl
394 END IF
395*
396 DO 20 i = 1, n
397 IF( abs( d( i ) ).LT.eps ) THEN
398 d( i ) = sign( eps, d( i ) )
399 END IF
400 20 CONTINUE
401*
402 start = 1
403 sqre = 0
404*
405 DO 30 i = 1, nm1
406 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
407*
408* Subproblem found. First determine its size and then
409* apply divide and conquer on it.
410*
411 IF( i.LT.nm1 ) THEN
412*
413* A subproblem with E(I) small for I < NM1.
414*
415 nsize = i - start + 1
416 ELSE IF( abs( e( i ) ).GE.eps ) THEN
417*
418* A subproblem with E(NM1) not too small but I = NM1.
419*
420 nsize = n - start + 1
421 ELSE
422*
423* A subproblem with E(NM1) small. This implies an
424* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
425* first.
426*
427 nsize = i - start + 1
428 IF( icompq.EQ.2 ) THEN
429 u( n, n ) = sign( one, d( n ) )
430 vt( n, n ) = one
431 ELSE IF( icompq.EQ.1 ) THEN
432 q( n+( qstart-1 )*n ) = sign( one, d( n ) )
433 q( n+( smlsiz+qstart-1 )*n ) = one
434 END IF
435 d( n ) = abs( d( n ) )
436 END IF
437 IF( icompq.EQ.2 ) THEN
438 CALL dlasd0( nsize, sqre, d( start ), e( start ),
439 $ u( start, start ), ldu, vt( start, start ),
440 $ ldvt, smlsiz, iwork, work( wstart ), info )
441 ELSE
442 CALL dlasda( icompq, smlsiz, nsize, sqre, d( start ),
443 $ e( start ), q( start+( iu+qstart-2 )*n ), n,
444 $ q( start+( ivt+qstart-2 )*n ),
445 $ iq( start+k*n ), q( start+( difl+qstart-2 )*
446 $ n ), q( start+( difr+qstart-2 )*n ),
447 $ q( start+( z+qstart-2 )*n ),
448 $ q( start+( poles+qstart-2 )*n ),
449 $ iq( start+givptr*n ), iq( start+givcol*n ),
450 $ n, iq( start+perm*n ),
451 $ q( start+( givnum+qstart-2 )*n ),
452 $ q( start+( ic+qstart-2 )*n ),
453 $ q( start+( is+qstart-2 )*n ),
454 $ work( wstart ), iwork, info )
455 END IF
456 IF( info.NE.0 ) THEN
457 RETURN
458 END IF
459 start = i + 1
460 END IF
461 30 CONTINUE
462*
463* Unscale
464*
465 CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, ierr )
466 40 CONTINUE
467*
468* Use Selection Sort to minimize swaps of singular vectors
469*
470 DO 60 ii = 2, n
471 i = ii - 1
472 kk = i
473 p = d( i )
474 DO 50 j = ii, n
475 IF( d( j ).GT.p ) THEN
476 kk = j
477 p = d( j )
478 END IF
479 50 CONTINUE
480 IF( kk.NE.i ) THEN
481 d( kk ) = d( i )
482 d( i ) = p
483 IF( icompq.EQ.1 ) THEN
484 iq( i ) = kk
485 ELSE IF( icompq.EQ.2 ) THEN
486 CALL dswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
487 CALL dswap( n, vt( i, 1 ), ldvt, vt( kk, 1 ), ldvt )
488 END IF
489 ELSE IF( icompq.EQ.1 ) THEN
490 iq( i ) = i
491 END IF
492 60 CONTINUE
493*
494* If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
495*
496 IF( icompq.EQ.1 ) THEN
497 IF( iuplo.EQ.1 ) THEN
498 iq( n ) = 1
499 ELSE
500 iq( n ) = 0
501 END IF
502 END IF
503*
504* If B is lower bidiagonal, update U by those Givens rotations
505* which rotated B to be upper bidiagonal
506*
507 IF( ( iuplo.EQ.2 ) .AND. ( icompq.EQ.2 ) )
508 $ CALL dlasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u, ldu )
509*
510 RETURN
511*
512* End of DBDSDC
513*
514 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dbdsdc(uplo, compq, n, d, e, u, ldu, vt, ldvt, q, iq, work, iwork, info)
DBDSDC
Definition dbdsdc.f:198
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 dlasd0(n, sqre, d, e, u, ldu, vt, ldvt, smlsiz, iwork, work, info)
DLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and of...
Definition dlasd0.f:152
subroutine dlasda(icompq, smlsiz, n, sqre, d, e, u, ldu, vt, k, difl, difr, z, poles, givptr, givcol, ldgcol, perm, givnum, c, s, work, iwork, info)
DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition dlasda.f:273
subroutine dlasdq(uplo, sqre, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
Definition dlasdq.f:211
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
subroutine dlasr(side, pivot, direct, m, n, c, s, a, lda)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition dlasr.f:199
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82