LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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

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

Definition at line 196 of file sbdsdc.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 REAL 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 REAL ZERO, ONE, TWO
220 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+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 REAL CS, EPS, ORGNRM, P, R, SN
228* ..
229* .. External Functions ..
230 LOGICAL LSAME
231 INTEGER ILAENV
232 REAL SLAMCH, SLANST
233 EXTERNAL slamch, slanst, ilaenv, lsame
234* ..
235* .. External Subroutines ..
236 EXTERNAL scopy, slartg, slascl, slasd0, slasda, slasdq,
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, u,
344 \$ ldu, work( wstart ), info )
345 ELSE IF( icompq.EQ.1 ) THEN
346 iu = 1
347 ivt = iu + n
348 CALL slaset( 'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
349 \$ n )
350 CALL slaset( 'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
351 \$ n )
352 CALL slasdq( '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 slaset( 'A', n, n, zero, one, u, ldu )
363 CALL slaset( 'A', n, n, zero, one, vt, ldvt )
364 END IF
365*
366* Scale.
367*
368 orgnrm = slanst( 'M', n, d, e )
369 IF( orgnrm.EQ.zero )
370 \$ RETURN
371 CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
372 CALL slascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
373*
374 eps = slamch( 'Epsilon' )
375*
376 mlvl = int( log( real( n ) / real( 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 slasd0( 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 slasda( 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 slascl( '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 sswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
487 CALL sswap( 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 slasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u, ldu )
509*
510 RETURN
511*
512* End of SBDSDC
513*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:100
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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:143
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 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:273
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:211
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:110
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:199
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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: