205 SUBROUTINE dbdsdc( 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 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( ldu, * ),
220 $ vt( ldvt, * ), work( * )
229 DOUBLE PRECISION ZERO, ONE, TWO
230 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+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 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
242 DOUBLE PRECISION DLAMCH, DLANST
243 EXTERNAL lsame, ilaenv, dlamch, dlanst
250 INTRINSIC abs, dble, 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(
'DBDSDC', -info )
294 smlsiz = ilaenv( 9,
'DBDSDC',
' ', 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 dcopy( n, d, 1, q( 1 ), 1 )
315 CALL dcopy( n-1, e, 1, q( n+1 ), 1 )
317 IF( iuplo.EQ.2 )
THEN
321 CALL dlartg( 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 dlasdq(
'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 dlaset(
'A', n, n, zero, one, u, ldu )
349 CALL dlaset(
'A', n, n, zero, one, vt, ldvt )
350 CALL dlasdq(
'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 dlaset(
'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
357 CALL dlaset(
'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
359 CALL dlasdq(
'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 dlaset(
'A', n, n, zero, one, u, ldu )
370 CALL dlaset(
'A', n, n, zero, one, vt, ldvt )
375 orgnrm = dlanst(
'M', n, d, e )
378 CALL dlascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
379 CALL dlascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
381 eps = (0.9d+0)*dlamch(
'Epsilon' )
383 mlvl = int( log( dble( n ) / dble( 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 dlasd0( nsize, sqre, d( start ), e( start ),
446 $ u( start, start ), ldu, vt( start, start ),
447 $ ldvt, smlsiz, iwork, work( wstart ), info )
449 CALL dlasda( 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 dlascl(
'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 dswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
494 CALL dswap( 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 dlasr(
'L',
'V',
'B', n, n, work( 1 ), work( n ), u, ldu )
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...
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
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...
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
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.
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
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...
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...