LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dbdsdc()

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.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 203 of file dbdsdc.f.

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,
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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlanst.f:100
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:111
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: