LAPACK 3.12.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.

 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 196 of file dbdsdc.f.

198*
199* -- LAPACK computational routine --
200* -- LAPACK is a software package provided by Univ. of Tennessee, --
201* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
202*
203* .. Scalar Arguments ..
204 CHARACTER COMPQ, UPLO
205 INTEGER INFO, LDU, LDVT, N
206* ..
207* .. Array Arguments ..
208 INTEGER IQ( * ), IWORK( * )
209 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
210 $ VT( LDVT, * ), WORK( * )
211* ..
212*
213* =====================================================================
214* Changed dimension statement in comment describing E from (N) to
215* (N-1). Sven, 17 Feb 05.
216* =====================================================================
217*
218* .. Parameters ..
219 DOUBLE PRECISION ZERO, ONE, TWO
220 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
221* ..
222* .. Local Scalars ..
223 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
224 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
225 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
226 $ SMLSZP, SQRE, START, WSTART, Z
227 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
228* ..
229* .. External Functions ..
230 LOGICAL LSAME
231 INTEGER ILAENV
232 DOUBLE PRECISION DLAMCH, DLANST
233 EXTERNAL lsame, ilaenv, dlamch, dlanst
234* ..
235* .. External Subroutines ..
236 EXTERNAL dcopy, dlartg, dlascl, dlasd0, dlasda, dlasdq,
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC abs, dble, int, log, sign
241* ..
242* .. Executable Statements ..
243*
244* Test the input parameters.
245*
246 info = 0
247*
248 iuplo = 0
249 IF( lsame( uplo, 'U' ) )
250 $ iuplo = 1
251 IF( lsame( uplo, 'L' ) )
252 $ iuplo = 2
253 IF( lsame( compq, 'N' ) ) THEN
254 icompq = 0
255 ELSE IF( lsame( compq, 'P' ) ) THEN
256 icompq = 1
257 ELSE IF( lsame( compq, 'I' ) ) THEN
258 icompq = 2
259 ELSE
260 icompq = -1
261 END IF
262 IF( iuplo.EQ.0 ) THEN
263 info = -1
264 ELSE IF( icompq.LT.0 ) THEN
265 info = -2
266 ELSE IF( n.LT.0 ) THEN
267 info = -3
268 ELSE IF( ( ldu.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldu.LT.
269 $ n ) ) ) THEN
270 info = -7
271 ELSE IF( ( ldvt.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldvt.LT.
272 $ n ) ) ) THEN
273 info = -9
274 END IF
275 IF( info.NE.0 ) THEN
276 CALL xerbla( 'DBDSDC', -info )
277 RETURN
278 END IF
279*
280* Quick return if possible
281*
282 IF( n.EQ.0 )
283 $ RETURN
284 smlsiz = ilaenv( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
285 IF( n.EQ.1 ) THEN
286 IF( icompq.EQ.1 ) THEN
287 q( 1 ) = sign( one, d( 1 ) )
288 q( 1+smlsiz*n ) = one
289 ELSE IF( icompq.EQ.2 ) THEN
290 u( 1, 1 ) = sign( one, d( 1 ) )
291 vt( 1, 1 ) = one
292 END IF
293 d( 1 ) = abs( d( 1 ) )
294 RETURN
295 END IF
296 nm1 = n - 1
297*
298* If matrix lower bidiagonal, rotate to be upper bidiagonal
299* by applying Givens rotations on the left
300*
301 wstart = 1
302 qstart = 3
303 IF( icompq.EQ.1 ) THEN
304 CALL dcopy( n, d, 1, q( 1 ), 1 )
305 CALL dcopy( n-1, e, 1, q( n+1 ), 1 )
306 END IF
307 IF( iuplo.EQ.2 ) THEN
308 qstart = 5
309 IF( icompq .EQ. 2 ) wstart = 2*n - 1
310 DO 10 i = 1, n - 1
311 CALL dlartg( d( i ), e( i ), cs, sn, r )
312 d( i ) = r
313 e( i ) = sn*d( i+1 )
314 d( i+1 ) = cs*d( i+1 )
315 IF( icompq.EQ.1 ) THEN
316 q( i+2*n ) = cs
317 q( i+3*n ) = sn
318 ELSE IF( icompq.EQ.2 ) THEN
319 work( i ) = cs
320 work( nm1+i ) = -sn
321 END IF
322 10 CONTINUE
323 END IF
324*
325* If ICOMPQ = 0, use DLASDQ to compute the singular values.
326*
327 IF( icompq.EQ.0 ) THEN
328* Ignore WSTART, instead using WORK( 1 ), since the two vectors
329* for CS and -SN above are added only if ICOMPQ == 2,
330* and adding them exceeds documented WORK size of 4*n.
331 CALL dlasdq( 'U', 0, n, 0, 0, 0, d, e, vt, ldvt, u, ldu, u,
332 $ ldu, work( 1 ), info )
333 GO TO 40
334 END IF
335*
336* If N is smaller than the minimum divide size SMLSIZ, then solve
337* the problem with another solver.
338*
339 IF( n.LE.smlsiz ) THEN
340 IF( icompq.EQ.2 ) THEN
341 CALL dlaset( 'A', n, n, zero, one, u, ldu )
342 CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
343 CALL dlasdq( 'U', 0, n, n, n, 0, d, e, vt, ldvt, u, ldu, u,
344 $ ldu, work( wstart ), info )
345 ELSE IF( icompq.EQ.1 ) THEN
346 iu = 1
347 ivt = iu + n
348 CALL dlaset( 'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
349 $ n )
350 CALL dlaset( 'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
351 $ n )
352 CALL dlasdq( 'U', 0, n, n, n, 0, d, e,
353 $ q( ivt+( qstart-1 )*n ), n,
354 $ q( iu+( qstart-1 )*n ), n,
355 $ q( iu+( qstart-1 )*n ), n, work( wstart ),
356 $ info )
357 END IF
358 GO TO 40
359 END IF
360*
361 IF( icompq.EQ.2 ) THEN
362 CALL dlaset( 'A', n, n, zero, one, u, ldu )
363 CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
364 END IF
365*
366* Scale.
367*
368 orgnrm = dlanst( 'M', n, d, e )
369 IF( orgnrm.EQ.zero )
370 $ RETURN
371 CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
372 CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
373*
374 eps = (0.9d+0)*dlamch( 'Epsilon' )
375*
376 mlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
377 smlszp = smlsiz + 1
378*
379 IF( icompq.EQ.1 ) THEN
380 iu = 1
381 ivt = 1 + smlsiz
382 difl = ivt + smlszp
383 difr = difl + mlvl
384 z = difr + mlvl*2
385 ic = z + mlvl
386 is = ic + 1
387 poles = is + 1
388 givnum = poles + 2*mlvl
389*
390 k = 1
391 givptr = 2
392 perm = 3
393 givcol = perm + mlvl
394 END IF
395*
396 DO 20 i = 1, n
397 IF( abs( d( i ) ).LT.eps ) THEN
398 d( i ) = sign( eps, d( i ) )
399 END IF
400 20 CONTINUE
401*
402 start = 1
403 sqre = 0
404*
405 DO 30 i = 1, nm1
406 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
407*
408* Subproblem found. First determine its size and then
409* apply divide and conquer on it.
410*
411 IF( i.LT.nm1 ) THEN
412*
413* A subproblem with E(I) small for I < NM1.
414*
415 nsize = i - start + 1
416 ELSE IF( abs( e( i ) ).GE.eps ) THEN
417*
418* A subproblem with E(NM1) not too small but I = NM1.
419*
420 nsize = n - start + 1
421 ELSE
422*
423* A subproblem with E(NM1) small. This implies an
424* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
425* first.
426*
427 nsize = i - start + 1
428 IF( icompq.EQ.2 ) THEN
429 u( n, n ) = sign( one, d( n ) )
430 vt( n, n ) = one
431 ELSE IF( icompq.EQ.1 ) THEN
432 q( n+( qstart-1 )*n ) = sign( one, d( n ) )
433 q( n+( smlsiz+qstart-1 )*n ) = one
434 END IF
435 d( n ) = abs( d( n ) )
436 END IF
437 IF( icompq.EQ.2 ) THEN
438 CALL dlasd0( nsize, sqre, d( start ), e( start ),
439 $ u( start, start ), ldu, vt( start, start ),
440 $ ldvt, smlsiz, iwork, work( wstart ), info )
441 ELSE
442 CALL dlasda( icompq, smlsiz, nsize, sqre, d( start ),
443 $ e( start ), q( start+( iu+qstart-2 )*n ), n,
444 $ q( start+( ivt+qstart-2 )*n ),
445 $ iq( start+k*n ), q( start+( difl+qstart-2 )*
446 $ n ), q( start+( difr+qstart-2 )*n ),
447 $ q( start+( z+qstart-2 )*n ),
448 $ q( start+( poles+qstart-2 )*n ),
449 $ iq( start+givptr*n ), iq( start+givcol*n ),
450 $ n, iq( start+perm*n ),
451 $ q( start+( givnum+qstart-2 )*n ),
452 $ q( start+( ic+qstart-2 )*n ),
453 $ q( start+( is+qstart-2 )*n ),
454 $ work( wstart ), iwork, info )
455 END IF
456 IF( info.NE.0 ) THEN
457 RETURN
458 END IF
459 start = i + 1
460 END IF
461 30 CONTINUE
462*
463* Unscale
464*
465 CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, ierr )
466 40 CONTINUE
467*
468* Use Selection Sort to minimize swaps of singular vectors
469*
470 DO 60 ii = 2, n
471 i = ii - 1
472 kk = i
473 p = d( i )
474 DO 50 j = ii, n
475 IF( d( j ).GT.p ) THEN
476 kk = j
477 p = d( j )
478 END IF
479 50 CONTINUE
480 IF( kk.NE.i ) THEN
481 d( kk ) = d( i )
482 d( i ) = p
483 IF( icompq.EQ.1 ) THEN
484 iq( i ) = kk
485 ELSE IF( icompq.EQ.2 ) THEN
486 CALL dswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
487 CALL dswap( n, vt( i, 1 ), ldvt, vt( kk, 1 ), ldvt )
488 END IF
489 ELSE IF( icompq.EQ.1 ) THEN
490 iq( i ) = i
491 END IF
492 60 CONTINUE
493*
494* If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
495*
496 IF( icompq.EQ.1 ) THEN
497 IF( iuplo.EQ.1 ) THEN
498 iq( n ) = 1
499 ELSE
500 iq( n ) = 0
501 END IF
502 END IF
503*
504* If B is lower bidiagonal, update U by those Givens rotations
505* which rotated B to be upper bidiagonal
506*
507 IF( ( iuplo.EQ.2 ) .AND. ( icompq.EQ.2 ) )
508 $ CALL dlasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u, ldu )
509*
510 RETURN
511*
512* End of DBDSDC
513*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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 dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 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:152
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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: