LAPACK 3.12.1
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 194 of file dbdsdc.f.

197*
198* -- LAPACK computational routine --
199* -- LAPACK is a software package provided by Univ. of Tennessee, --
200* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201*
202* .. Scalar Arguments ..
203 CHARACTER COMPQ, UPLO
204 INTEGER INFO, LDU, LDVT, N
205* ..
206* .. Array Arguments ..
207 INTEGER IQ( * ), IWORK( * )
208 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
209 $ VT( LDVT, * ), WORK( * )
210* ..
211*
212* =====================================================================
213* Changed dimension statement in comment describing E from (N) to
214* (N-1). Sven, 17 Feb 05.
215* =====================================================================
216*
217* .. Parameters ..
218 DOUBLE PRECISION ZERO, ONE, TWO
219 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
220* ..
221* .. Local Scalars ..
222 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
223 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
224 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
225 $ SMLSZP, SQRE, START, WSTART, Z
226 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
227* ..
228* .. External Functions ..
229 LOGICAL LSAME
230 INTEGER ILAENV
231 DOUBLE PRECISION DLAMCH, DLANST
232 EXTERNAL lsame, ilaenv, dlamch, dlanst
233* ..
234* .. External Subroutines ..
235 EXTERNAL dcopy, dlartg, dlascl, dlasd0, dlasda,
236 $ dlasdq,
237 $ dlaset, dlasr, dswap, xerbla
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,
344 $ u,
345 $ ldu, work( wstart ), info )
346 ELSE IF( icompq.EQ.1 ) THEN
347 iu = 1
348 ivt = iu + n
349 CALL dlaset( 'A', n, n, zero, one,
350 $ q( iu+( qstart-1 )*n ),
351 $ n )
352 CALL dlaset( 'A', n, n, zero, one,
353 $ q( ivt+( qstart-1 )*n ),
354 $ n )
355 CALL dlasdq( 'U', 0, n, n, n, 0, d, e,
356 $ q( ivt+( qstart-1 )*n ), n,
357 $ q( iu+( qstart-1 )*n ), n,
358 $ q( iu+( qstart-1 )*n ), n, work( wstart ),
359 $ info )
360 END IF
361 GO TO 40
362 END IF
363*
364 IF( icompq.EQ.2 ) THEN
365 CALL dlaset( 'A', n, n, zero, one, u, ldu )
366 CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
367 END IF
368*
369* Scale.
370*
371 orgnrm = dlanst( 'M', n, d, e )
372 IF( orgnrm.EQ.zero )
373 $ RETURN
374 CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
375 CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
376*
377 eps = (0.9d+0)*dlamch( 'Epsilon' )
378*
379 mlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
380 smlszp = smlsiz + 1
381*
382 IF( icompq.EQ.1 ) THEN
383 iu = 1
384 ivt = 1 + smlsiz
385 difl = ivt + smlszp
386 difr = difl + mlvl
387 z = difr + mlvl*2
388 ic = z + mlvl
389 is = ic + 1
390 poles = is + 1
391 givnum = poles + 2*mlvl
392*
393 k = 1
394 givptr = 2
395 perm = 3
396 givcol = perm + mlvl
397 END IF
398*
399 DO 20 i = 1, n
400 IF( abs( d( i ) ).LT.eps ) THEN
401 d( i ) = sign( eps, d( i ) )
402 END IF
403 20 CONTINUE
404*
405 start = 1
406 sqre = 0
407*
408 DO 30 i = 1, nm1
409 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
410*
411* Subproblem found. First determine its size and then
412* apply divide and conquer on it.
413*
414 IF( i.LT.nm1 ) THEN
415*
416* A subproblem with E(I) small for I < NM1.
417*
418 nsize = i - start + 1
419 ELSE IF( abs( e( i ) ).GE.eps ) THEN
420*
421* A subproblem with E(NM1) not too small but I = NM1.
422*
423 nsize = n - start + 1
424 ELSE
425*
426* A subproblem with E(NM1) small. This implies an
427* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
428* first.
429*
430 nsize = i - start + 1
431 IF( icompq.EQ.2 ) THEN
432 u( n, n ) = sign( one, d( n ) )
433 vt( n, n ) = one
434 ELSE IF( icompq.EQ.1 ) THEN
435 q( n+( qstart-1 )*n ) = sign( one, d( n ) )
436 q( n+( smlsiz+qstart-1 )*n ) = one
437 END IF
438 d( n ) = abs( d( n ) )
439 END IF
440 IF( icompq.EQ.2 ) THEN
441 CALL dlasd0( nsize, sqre, d( start ), e( start ),
442 $ u( start, start ), ldu, vt( start, start ),
443 $ ldvt, smlsiz, iwork, work( wstart ), info )
444 ELSE
445 CALL dlasda( icompq, smlsiz, nsize, sqre, d( start ),
446 $ e( start ), q( start+( iu+qstart-2 )*n ), n,
447 $ q( start+( ivt+qstart-2 )*n ),
448 $ iq( start+k*n ), q( start+( difl+qstart-2 )*
449 $ n ), q( start+( difr+qstart-2 )*n ),
450 $ q( start+( z+qstart-2 )*n ),
451 $ q( start+( poles+qstart-2 )*n ),
452 $ iq( start+givptr*n ), iq( start+givcol*n ),
453 $ n, iq( start+perm*n ),
454 $ q( start+( givnum+qstart-2 )*n ),
455 $ q( start+( ic+qstart-2 )*n ),
456 $ q( start+( is+qstart-2 )*n ),
457 $ work( wstart ), iwork, info )
458 END IF
459 IF( info.NE.0 ) THEN
460 RETURN
461 END IF
462 start = i + 1
463 END IF
464 30 CONTINUE
465*
466* Unscale
467*
468 CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, ierr )
469 40 CONTINUE
470*
471* Use Selection Sort to minimize swaps of singular vectors
472*
473 DO 60 ii = 2, n
474 i = ii - 1
475 kk = i
476 p = d( i )
477 DO 50 j = ii, n
478 IF( d( j ).GT.p ) THEN
479 kk = j
480 p = d( j )
481 END IF
482 50 CONTINUE
483 IF( kk.NE.i ) THEN
484 d( kk ) = d( i )
485 d( i ) = p
486 IF( icompq.EQ.1 ) THEN
487 iq( i ) = kk
488 ELSE IF( icompq.EQ.2 ) THEN
489 CALL dswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
490 CALL dswap( n, vt( i, 1 ), ldvt, vt( kk, 1 ), ldvt )
491 END IF
492 ELSE IF( icompq.EQ.1 ) THEN
493 iq( i ) = i
494 END IF
495 60 CONTINUE
496*
497* If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
498*
499 IF( icompq.EQ.1 ) THEN
500 IF( iuplo.EQ.1 ) THEN
501 iq( n ) = 1
502 ELSE
503 iq( n ) = 0
504 END IF
505 END IF
506*
507* If B is lower bidiagonal, update U by those Givens rotations
508* which rotated B to be upper bidiagonal
509*
510 IF( ( iuplo.EQ.2 ) .AND. ( icompq.EQ.2 ) )
511 $ CALL dlasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u,
512 $ ldu )
513*
514 RETURN
515*
516* End of DBDSDC
517*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
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:98
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlascl
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:142
dlasd0
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:151
dlasda
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:272
dlasdq
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:210
dlaset
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:108
dlasr
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:197
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dswap
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: