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

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

Definition at line 194 of file sbdsdc.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 REAL 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 REAL ZERO, ONE, TWO
219 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+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 REAL CS, EPS, ORGNRM, P, R, SN
227* ..
228* .. External Functions ..
229 LOGICAL LSAME
230 INTEGER ILAENV
231 REAL SLAMCH, SLANST
232 EXTERNAL slamch, slanst, ilaenv, lsame
233* ..
234* .. External Subroutines ..
235 EXTERNAL scopy, slartg, slascl, slasd0, slasda,
236 $ slasdq,
237 $ slaset, slasr, sswap, xerbla
238* ..
239* .. Intrinsic Functions ..
240 INTRINSIC real, abs, 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( 'SBDSDC', -info )
277 RETURN
278 END IF
279*
280* Quick return if possible
281*
282 IF( n.EQ.0 )
283 $ RETURN
284 smlsiz = ilaenv( 9, 'SBDSDC', ' ', 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 scopy( n, d, 1, q( 1 ), 1 )
305 CALL scopy( 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 slartg( 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 SLASDQ 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 slasdq( '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 slaset( 'A', n, n, zero, one, u, ldu )
342 CALL slaset( 'A', n, n, zero, one, vt, ldvt )
343 CALL slasdq( '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 slaset( 'A', n, n, zero, one,
350 $ q( iu+( qstart-1 )*n ),
351 $ n )
352 CALL slaset( 'A', n, n, zero, one,
353 $ q( ivt+( qstart-1 )*n ),
354 $ n )
355 CALL slasdq( '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 slaset( 'A', n, n, zero, one, u, ldu )
366 CALL slaset( 'A', n, n, zero, one, vt, ldvt )
367 END IF
368*
369* Scale.
370*
371 orgnrm = slanst( 'M', n, d, e )
372 IF( orgnrm.EQ.zero )
373 $ RETURN
374 CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
375 CALL slascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
376*
377 eps = slamch( 'Epsilon' )
378*
379 mlvl = int( log( real( n ) / real( 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 slasd0( 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 slasda( 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 slascl( '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 sswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
490 CALL sswap( 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 slasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u,
512 $ ldu )
513*
514 RETURN
515*
516* End of SBDSDC
517*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slanst
real function slanst(norm, n, d, e)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slanst.f:98
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
slascl
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:142
slasd0
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:151
slasda
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:272
slasdq
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:210
slaset
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:108
slasr
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:197
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: