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