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.