LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sbdsdc.f
Go to the documentation of this file.
1 *> \brief \b SBDSDC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SBDSDC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sbdsdc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sbdsdc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sbdsdc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SBDSDC( 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 * REAL D( * ), E( * ), Q( * ), U( LDU, * ),
31 * $ VT( LDVT, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> SBDSDC 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. SBDSDC 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 SLASD3 for details.
54 *>
55 *> The code currently calls SLASDQ 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 *> \date November 2011
195 *
196 *> \ingroup auxOTHERcomputational
197 *
198 *> \par Contributors:
199 * ==================
200 *>
201 *> Ming Gu and Huan Ren, Computer Science Division, University of
202 *> California at Berkeley, USA
203 *>
204 * =====================================================================
205  SUBROUTINE sbdsdc( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
206  $ work, iwork, info )
207 *
208 * -- LAPACK computational routine (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  CHARACTER compq, uplo
215  INTEGER info, ldu, ldvt, n
216 * ..
217 * .. Array Arguments ..
218  INTEGER iq( * ), iwork( * )
219  REAL d( * ), e( * ), q( * ), u( ldu, * ),
220  $ vt( ldvt, * ), work( * )
221 * ..
222 *
223 * =====================================================================
224 * Changed dimension statement in comment describing E from (N) to
225 * (N-1). Sven, 17 Feb 05.
226 * =====================================================================
227 *
228 * .. Parameters ..
229  REAL zero, one, two
230  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
231 * ..
232 * .. Local Scalars ..
233  INTEGER difl, difr, givcol, givnum, givptr, i, ic,
234  $ icompq, ierr, ii, is, iu, iuplo, ivt, j, k, kk,
235  $ mlvl, nm1, nsize, perm, poles, qstart, smlsiz,
236  $ smlszp, sqre, start, wstart, z
237  REAL cs, eps, orgnrm, p, r, sn
238 * ..
239 * .. External Functions ..
240  LOGICAL lsame
241  INTEGER ilaenv
242  REAL slamch, slanst
243  EXTERNAL slamch, slanst, ilaenv, lsame
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL scopy, slartg, slascl, slasd0, slasda, slasdq,
247  $ slaset, slasr, sswap, xerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC REAL, abs, int, log, sign
251 * ..
252 * .. Executable Statements ..
253 *
254 * Test the input parameters.
255 *
256  info = 0
257 *
258  iuplo = 0
259  IF( lsame( uplo, 'U' ) )
260  $ iuplo = 1
261  IF( lsame( uplo, 'L' ) )
262  $ iuplo = 2
263  IF( lsame( compq, 'N' ) ) THEN
264  icompq = 0
265  ELSE IF( lsame( compq, 'P' ) ) THEN
266  icompq = 1
267  ELSE IF( lsame( compq, 'I' ) ) THEN
268  icompq = 2
269  ELSE
270  icompq = -1
271  END IF
272  IF( iuplo.EQ.0 ) THEN
273  info = -1
274  ELSE IF( icompq.LT.0 ) THEN
275  info = -2
276  ELSE IF( n.LT.0 ) THEN
277  info = -3
278  ELSE IF( ( ldu.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldu.LT.
279  $ n ) ) ) THEN
280  info = -7
281  ELSE IF( ( ldvt.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldvt.LT.
282  $ n ) ) ) THEN
283  info = -9
284  END IF
285  IF( info.NE.0 ) THEN
286  CALL xerbla( 'SBDSDC', -info )
287  return
288  END IF
289 *
290 * Quick return if possible
291 *
292  IF( n.EQ.0 )
293  $ return
294  smlsiz = ilaenv( 9, 'SBDSDC', ' ', 0, 0, 0, 0 )
295  IF( n.EQ.1 ) THEN
296  IF( icompq.EQ.1 ) THEN
297  q( 1 ) = sign( one, d( 1 ) )
298  q( 1+smlsiz*n ) = one
299  ELSE IF( icompq.EQ.2 ) THEN
300  u( 1, 1 ) = sign( one, d( 1 ) )
301  vt( 1, 1 ) = one
302  END IF
303  d( 1 ) = abs( d( 1 ) )
304  return
305  END IF
306  nm1 = n - 1
307 *
308 * If matrix lower bidiagonal, rotate to be upper bidiagonal
309 * by applying Givens rotations on the left
310 *
311  wstart = 1
312  qstart = 3
313  IF( icompq.EQ.1 ) THEN
314  CALL scopy( n, d, 1, q( 1 ), 1 )
315  CALL scopy( n-1, e, 1, q( n+1 ), 1 )
316  END IF
317  IF( iuplo.EQ.2 ) THEN
318  qstart = 5
319  wstart = 2*n - 1
320  DO 10 i = 1, n - 1
321  CALL slartg( d( i ), e( i ), cs, sn, r )
322  d( i ) = r
323  e( i ) = sn*d( i+1 )
324  d( i+1 ) = cs*d( i+1 )
325  IF( icompq.EQ.1 ) THEN
326  q( i+2*n ) = cs
327  q( i+3*n ) = sn
328  ELSE IF( icompq.EQ.2 ) THEN
329  work( i ) = cs
330  work( nm1+i ) = -sn
331  END IF
332  10 continue
333  END IF
334 *
335 * If ICOMPQ = 0, use SLASDQ to compute the singular values.
336 *
337  IF( icompq.EQ.0 ) THEN
338  CALL slasdq( 'U', 0, n, 0, 0, 0, d, e, vt, ldvt, u, ldu, u,
339  $ ldu, work( wstart ), 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 slaset( 'A', n, n, zero, one, u, ldu )
349  CALL slaset( 'A', n, n, zero, one, vt, ldvt )
350  CALL slasdq( '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 slaset( 'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
356  $ n )
357  CALL slaset( 'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
358  $ n )
359  CALL slasdq( '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 slaset( 'A', n, n, zero, one, u, ldu )
370  CALL slaset( 'A', n, n, zero, one, vt, ldvt )
371  END IF
372 *
373 * Scale.
374 *
375  orgnrm = slanst( 'M', n, d, e )
376  IF( orgnrm.EQ.zero )
377  $ return
378  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
379  CALL slascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
380 *
381  eps = slamch( 'Epsilon' )
382 *
383  mlvl = int( log( REAL( N ) / REAL( 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 slasd0( 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 slasda( 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 slascl( '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 sswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
494  CALL sswap( 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 slasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u, ldu )
516 *
517  return
518 *
519 * End of SBDSDC
520 *
521  END