LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sbdsdc()

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

SBDSDC

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

Purpose:
 SBDSDC 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. SBDSDC 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 SLASD3 for details.

 The code currently calls SLASDQ 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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
June 2016
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 207 of file sbdsdc.f.

207 *
208 * -- LAPACK computational routine (version 3.7.1) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * June 2016
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  IF( icompq .EQ. 2 ) 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 * Ignore WSTART, instead using WORK( 1 ), since the two vectors
339 * for CS and -SN above are added only if ICOMPQ == 2,
340 * and adding them exceeds documented WORK size of 4*n.
341  CALL slasdq( 'U', 0, n, 0, 0, 0, d, e, vt, ldvt, u, ldu, u,
342  $ ldu, work( 1 ), info )
343  GO TO 40
344  END IF
345 *
346 * If N is smaller than the minimum divide size SMLSIZ, then solve
347 * the problem with another solver.
348 *
349  IF( n.LE.smlsiz ) THEN
350  IF( icompq.EQ.2 ) THEN
351  CALL slaset( 'A', n, n, zero, one, u, ldu )
352  CALL slaset( 'A', n, n, zero, one, vt, ldvt )
353  CALL slasdq( 'U', 0, n, n, n, 0, d, e, vt, ldvt, u, ldu, u,
354  $ ldu, work( wstart ), info )
355  ELSE IF( icompq.EQ.1 ) THEN
356  iu = 1
357  ivt = iu + n
358  CALL slaset( 'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
359  $ n )
360  CALL slaset( 'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
361  $ n )
362  CALL slasdq( 'U', 0, n, n, n, 0, d, e,
363  $ q( ivt+( qstart-1 )*n ), n,
364  $ q( iu+( qstart-1 )*n ), n,
365  $ q( iu+( qstart-1 )*n ), n, work( wstart ),
366  $ info )
367  END IF
368  GO TO 40
369  END IF
370 *
371  IF( icompq.EQ.2 ) THEN
372  CALL slaset( 'A', n, n, zero, one, u, ldu )
373  CALL slaset( 'A', n, n, zero, one, vt, ldvt )
374  END IF
375 *
376 * Scale.
377 *
378  orgnrm = slanst( 'M', n, d, e )
379  IF( orgnrm.EQ.zero )
380  $ RETURN
381  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
382  CALL slascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
383 *
384  eps = slamch( 'Epsilon' )
385 *
386  mlvl = int( log( REAL( N ) / REAL( SMLSIZ+1 ) ) / log( two ) ) + 1
387  smlszp = smlsiz + 1
388 *
389  IF( icompq.EQ.1 ) THEN
390  iu = 1
391  ivt = 1 + smlsiz
392  difl = ivt + smlszp
393  difr = difl + mlvl
394  z = difr + mlvl*2
395  ic = z + mlvl
396  is = ic + 1
397  poles = is + 1
398  givnum = poles + 2*mlvl
399 *
400  k = 1
401  givptr = 2
402  perm = 3
403  givcol = perm + mlvl
404  END IF
405 *
406  DO 20 i = 1, n
407  IF( abs( d( i ) ).LT.eps ) THEN
408  d( i ) = sign( eps, d( i ) )
409  END IF
410  20 CONTINUE
411 *
412  start = 1
413  sqre = 0
414 *
415  DO 30 i = 1, nm1
416  IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
417 *
418 * Subproblem found. First determine its size and then
419 * apply divide and conquer on it.
420 *
421  IF( i.LT.nm1 ) THEN
422 *
423 * A subproblem with E(I) small for I < NM1.
424 *
425  nsize = i - start + 1
426  ELSE IF( abs( e( i ) ).GE.eps ) THEN
427 *
428 * A subproblem with E(NM1) not too small but I = NM1.
429 *
430  nsize = n - start + 1
431  ELSE
432 *
433 * A subproblem with E(NM1) small. This implies an
434 * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
435 * first.
436 *
437  nsize = i - start + 1
438  IF( icompq.EQ.2 ) THEN
439  u( n, n ) = sign( one, d( n ) )
440  vt( n, n ) = one
441  ELSE IF( icompq.EQ.1 ) THEN
442  q( n+( qstart-1 )*n ) = sign( one, d( n ) )
443  q( n+( smlsiz+qstart-1 )*n ) = one
444  END IF
445  d( n ) = abs( d( n ) )
446  END IF
447  IF( icompq.EQ.2 ) THEN
448  CALL slasd0( nsize, sqre, d( start ), e( start ),
449  $ u( start, start ), ldu, vt( start, start ),
450  $ ldvt, smlsiz, iwork, work( wstart ), info )
451  ELSE
452  CALL slasda( icompq, smlsiz, nsize, sqre, d( start ),
453  $ e( start ), q( start+( iu+qstart-2 )*n ), n,
454  $ q( start+( ivt+qstart-2 )*n ),
455  $ iq( start+k*n ), q( start+( difl+qstart-2 )*
456  $ n ), q( start+( difr+qstart-2 )*n ),
457  $ q( start+( z+qstart-2 )*n ),
458  $ q( start+( poles+qstart-2 )*n ),
459  $ iq( start+givptr*n ), iq( start+givcol*n ),
460  $ n, iq( start+perm*n ),
461  $ q( start+( givnum+qstart-2 )*n ),
462  $ q( start+( ic+qstart-2 )*n ),
463  $ q( start+( is+qstart-2 )*n ),
464  $ work( wstart ), iwork, info )
465  END IF
466  IF( info.NE.0 ) THEN
467  RETURN
468  END IF
469  start = i + 1
470  END IF
471  30 CONTINUE
472 *
473 * Unscale
474 *
475  CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, ierr )
476  40 CONTINUE
477 *
478 * Use Selection Sort to minimize swaps of singular vectors
479 *
480  DO 60 ii = 2, n
481  i = ii - 1
482  kk = i
483  p = d( i )
484  DO 50 j = ii, n
485  IF( d( j ).GT.p ) THEN
486  kk = j
487  p = d( j )
488  END IF
489  50 CONTINUE
490  IF( kk.NE.i ) THEN
491  d( kk ) = d( i )
492  d( i ) = p
493  IF( icompq.EQ.1 ) THEN
494  iq( i ) = kk
495  ELSE IF( icompq.EQ.2 ) THEN
496  CALL sswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
497  CALL sswap( n, vt( i, 1 ), ldvt, vt( kk, 1 ), ldvt )
498  END IF
499  ELSE IF( icompq.EQ.1 ) THEN
500  iq( i ) = i
501  END IF
502  60 CONTINUE
503 *
504 * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
505 *
506  IF( icompq.EQ.1 ) THEN
507  IF( iuplo.EQ.1 ) THEN
508  iq( n ) = 1
509  ELSE
510  iq( n ) = 0
511  END IF
512  END IF
513 *
514 * If B is lower bidiagonal, update U by those Givens rotations
515 * which rotated B to be upper bidiagonal
516 *
517  IF( ( iuplo.EQ.2 ) .AND. ( icompq.EQ.2 ) )
518  $ CALL slasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u, ldu )
519 *
520  RETURN
521 *
522 * End of SBDSDC
523 *
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: slasda.f:275
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: slasdq.f:213
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: slasr.f:201
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slasd0(N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK, WORK, INFO)
SLASD0 computes the singular values of a real upper bidiagonal n-by-m matrix B with diagonal d and of...
Definition: slasd0.f:152
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
real function slanst(NORM, N, D, E)
SLANST 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: slanst.f:102
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
Here is the call graph for this function:
Here is the caller graph for this function: