205 SUBROUTINE sbdsdc( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
206 $ work, iwork, info )
214 CHARACTER COMPQ, UPLO
215 INTEGER INFO, LDU, LDVT, N
218 INTEGER IQ( * ), IWORK( * )
219 REAL D( * ), E( * ), Q( * ), U( ldu, * ),
220 $ vt( ldvt, * ), work( * )
230 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
233 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
234 $ icompq, ierr, ii, is, iu, iuplo, ivt, j, k, kk,
235 $ mlvl, nm1, nsize, perm, poles, qstart, smlsiz,
236 $ smlszp, sqre, start, wstart, z
237 REAL CS, EPS, ORGNRM, P, R, SN
243 EXTERNAL slamch, slanst, ilaenv, lsame
250 INTRINSIC REAL, ABS, INT, LOG, SIGN
259 IF( lsame( uplo,
'U' ) )
261 IF( lsame( uplo,
'L' ) )
263 IF( lsame( compq,
'N' ) )
THEN
265 ELSE IF( lsame( compq,
'P' ) )
THEN
267 ELSE IF( lsame( compq,
'I' ) )
THEN
272 IF( iuplo.EQ.0 )
THEN
274 ELSE IF( icompq.LT.0 )
THEN
276 ELSE IF( n.LT.0 )
THEN
278 ELSE IF( ( ldu.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldu.LT.
281 ELSE IF( ( ldvt.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldvt.LT.
286 CALL
xerbla(
'SBDSDC', -info )
294 smlsiz = ilaenv( 9,
'SBDSDC',
' ', 0, 0, 0, 0 )
296 IF( icompq.EQ.1 )
THEN
297 q( 1 ) = sign( one, d( 1 ) )
298 q( 1+smlsiz*n ) = one
299 ELSE IF( icompq.EQ.2 )
THEN
300 u( 1, 1 ) = sign( one, d( 1 ) )
303 d( 1 ) = abs( d( 1 ) )
313 IF( icompq.EQ.1 )
THEN
314 CALL
scopy( n, d, 1, q( 1 ), 1 )
315 CALL
scopy( n-1, e, 1, q( n+1 ), 1 )
317 IF( iuplo.EQ.2 )
THEN
321 CALL
slartg( d( i ), e( i ), cs, sn, r )
324 d( i+1 ) = cs*d( i+1 )
325 IF( icompq.EQ.1 )
THEN
328 ELSE IF( icompq.EQ.2 )
THEN
337 IF( icompq.EQ.0 )
THEN
338 CALL
slasdq(
'U', 0, n, 0, 0, 0, d, e, vt, ldvt, u, ldu, u,
339 $ ldu, work( wstart ), info )
346 IF( n.LE.smlsiz )
THEN
347 IF( icompq.EQ.2 )
THEN
348 CALL
slaset(
'A', n, n, zero, one, u, ldu )
349 CALL
slaset(
'A', n, n, zero, one, vt, ldvt )
350 CALL
slasdq(
'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
355 CALL
slaset(
'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
357 CALL
slaset(
'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
359 CALL
slasdq(
'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 ),
368 IF( icompq.EQ.2 )
THEN
369 CALL
slaset(
'A', n, n, zero, one, u, ldu )
370 CALL
slaset(
'A', n, n, zero, one, vt, ldvt )
375 orgnrm = slanst(
'M', n, d, e )
378 CALL
slascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
379 CALL
slascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
381 eps = slamch(
'Epsilon' )
383 mlvl = int( log(
REAL( N ) /
REAL( SMLSIZ+1 ) ) / log( TWO ) ) + 1
386 IF( icompq.EQ.1 )
THEN
395 givnum = poles + 2*mlvl
404 IF( abs( d( i ) ).LT.eps )
THEN
405 d( i ) = sign( eps, d( i ) )
413 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
422 nsize = i - start + 1
423 ELSE IF( abs( e( i ) ).GE.eps )
THEN
427 nsize = n - start + 1
434 nsize = i - start + 1
435 IF( icompq.EQ.2 )
THEN
436 u( n, n ) = sign( one, d( n ) )
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
442 d( n ) = abs( d( n ) )
444 IF( icompq.EQ.2 )
THEN
445 CALL
slasd0( nsize, sqre, d( start ), e( start ),
446 $ u( start, start ), ldu, vt( start, start ),
447 $ ldvt, smlsiz, iwork, work( wstart ), info )
449 CALL
slasda( 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 )
472 CALL
slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, ierr )
482 IF( d( j ).GT.p )
THEN
490 IF( icompq.EQ.1 )
THEN
492 ELSE IF( icompq.EQ.2 )
THEN
493 CALL
sswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
494 CALL
sswap( n, vt( i, 1 ), ldvt, vt( kk, 1 ), ldvt )
496 ELSE IF( icompq.EQ.1 )
THEN
503 IF( icompq.EQ.1 )
THEN
504 IF( iuplo.EQ.1 )
THEN
514 IF( ( iuplo.EQ.2 ) .AND. ( icompq.EQ.2 ) )
515 $ CALL
slasr(
'L',
'V',
'B', n, n, work( 1 ), work( n ), u, ldu )
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...
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
subroutine xerbla(SRNAME, INFO)
XERBLA
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...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
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...
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.
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...
subroutine slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.