LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
dbdsdc.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dbdsdc (UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
 DBDSDC More...
 

Function/Subroutine Documentation

subroutine dbdsdc ( character  UPLO,
character  COMPQ,
integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldvt, * )  VT,
integer  LDVT,
double precision, dimension( * )  Q,
integer, dimension( * )  IQ,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DBDSDC

Download DBDSDC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DBDSDC computes the singular value decomposition (SVD) of a real
 N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
 using a divide and conquer method, where S is a diagonal matrix
 with non-negative diagonal elements (the singular values of B), and
 U and VT are orthogonal matrices of left and right singular vectors,
 respectively. DBDSDC can be used to compute all singular values,
 and optionally, singular vectors or singular vectors in compact form.

 This code makes very mild assumptions about floating point
 arithmetic. It will work on machines with a guard digit in
 add/subtract, or on those binary machines without guard digits
 which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.  See DLASD3 for details.

 The code currently calls DLASDQ if singular values only are desired.
 However, it can be slightly modified to compute singular values
 using the divide and conquer method.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  B is upper bidiagonal.
          = 'L':  B is lower bidiagonal.
[in]COMPQ
          COMPQ is CHARACTER*1
          Specifies whether singular vectors are to be computed
          as follows:
          = 'N':  Compute singular values only;
          = 'P':  Compute singular values and compute singular
                  vectors in compact form;
          = 'I':  Compute singular values and singular vectors.
[in]N
          N is INTEGER
          The order of the matrix B.  N >= 0.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the n diagonal elements of the bidiagonal matrix B.
          On exit, if INFO=0, the singular values of B.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N-1)
          On entry, the elements of E contain the offdiagonal
          elements of the bidiagonal matrix whose SVD is desired.
          On exit, E has been destroyed.
[out]U
          U is DOUBLE PRECISION array, dimension (LDU,N)
          If  COMPQ = 'I', then:
             On exit, if INFO = 0, U contains the left singular vectors
             of the bidiagonal matrix.
          For other values of COMPQ, U is not referenced.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= 1.
          If singular vectors are desired, then LDU >= max( 1, N ).
[out]VT
          VT is DOUBLE PRECISION array, dimension (LDVT,N)
          If  COMPQ = 'I', then:
             On exit, if INFO = 0, VT**T contains the right singular
             vectors of the bidiagonal matrix.
          For other values of COMPQ, VT is not referenced.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.  LDVT >= 1.
          If singular vectors are desired, then LDVT >= max( 1, N ).
[out]Q
          Q is DOUBLE PRECISION array, dimension (LDQ)
          If  COMPQ = 'P', then:
             On exit, if INFO = 0, Q and IQ contain the left
             and right singular vectors in a compact form,
             requiring O(N log N) space instead of 2*N**2.
             In particular, Q contains all the DOUBLE PRECISION data in
             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
             words of memory, where SMLSIZ is returned by ILAENV and
             is equal to the maximum size of the subproblems at the
             bottom of the computation tree (usually about 25).
          For other values of COMPQ, Q is not referenced.
[out]IQ
          IQ is INTEGER array, dimension (LDIQ)
          If  COMPQ = 'P', then:
             On exit, if INFO = 0, Q and IQ contain the left
             and right singular vectors in a compact form,
             requiring O(N log N) space instead of 2*N**2.
             In particular, IQ contains all INTEGER data in
             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
             words of memory, where SMLSIZ is returned by ILAENV and
             is equal to the maximum size of the subproblems at the
             bottom of the computation tree (usually about 25).
          For other values of COMPQ, IQ is not referenced.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          If COMPQ = 'N' then LWORK >= (4 * N).
          If COMPQ = 'P' then LWORK >= (6 * N).
          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
[out]IWORK
          IWORK is INTEGER array, dimension (8*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  The algorithm failed to compute a singular value.
                The update process of divide and conquer failed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 207 of file dbdsdc.f.

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  DOUBLE PRECISION 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  DOUBLE PRECISION zero, one, two
230  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+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  DOUBLE PRECISION cs, eps, orgnrm, p, r, sn
238 * ..
239 * .. External Functions ..
240  LOGICAL lsame
241  INTEGER ilaenv
242  DOUBLE PRECISION dlamch, dlanst
243  EXTERNAL lsame, ilaenv, dlamch, dlanst
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL dcopy, dlartg, dlascl, dlasd0, dlasda, dlasdq,
247  $ dlaset, dlasr, dswap, xerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, dble, 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( 'DBDSDC', -info )
287  RETURN
288  END IF
289 *
290 * Quick return if possible
291 *
292  IF( n.EQ.0 )
293  $ RETURN
294  smlsiz = ilaenv( 9, 'DBDSDC', ' ', 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 dcopy( n, d, 1, q( 1 ), 1 )
315  CALL dcopy( 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 dlartg( 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 DLASDQ to compute the singular values.
336 *
337  IF( icompq.EQ.0 ) THEN
338  CALL dlasdq( '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 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 *
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:276
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
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:213
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:141
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:201
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: dlanst.f:102
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
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:112
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:154

Here is the call graph for this function:

Here is the caller graph for this function: