262 SUBROUTINE dgesvdx( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
263 $ il, iu, ns, s, u, ldu, vt, ldvt, work,
264 $ lwork, iwork, info )
272 CHARACTER JOBU, JOBVT, RANGE
273 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
274 DOUBLE PRECISION VL, VU
278 DOUBLE PRECISION A( lda, * ), S( * ), U( ldu, * ),
279 $ vt( ldvt, * ), work( * )
285 DOUBLE PRECISION ZERO, ONE
286 parameter ( zero = 0.0d0, one = 1.0d0 )
289 CHARACTER JOBZ, RNGTGK
290 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
291 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
292 $ itau, itaup, itauq, itemp, itgkz, iutgk,
293 $ j, maxwrk, minmn, minwrk, mnthr
294 DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
297 DOUBLE PRECISION DUM( 1 )
307 DOUBLE PRECISION DLAMCH, DLANGE
308 EXTERNAL lsame, ilaenv, dlamch, dlange
311 INTRINSIC max, min, sqrt
319 abstol = 2*dlamch(
'S')
320 lquery = ( lwork.EQ.-1 )
323 wantu = lsame( jobu,
'V' )
324 wantvt = lsame( jobvt,
'V' )
325 IF( wantu .OR. wantvt )
THEN
330 alls = lsame( range,
'A' )
331 vals = lsame( range,
'V' )
332 inds = lsame( range,
'I' )
335 IF( .NOT.lsame( jobu,
'V' ) .AND.
336 $ .NOT.lsame( jobu,
'N' ) )
THEN
338 ELSE IF( .NOT.lsame( jobvt,
'V' ) .AND.
339 $ .NOT.lsame( jobvt,
'N' ) )
THEN
341 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) )
THEN
343 ELSE IF( m.LT.0 )
THEN
345 ELSE IF( n.LT.0 )
THEN
347 ELSE IF( m.GT.lda )
THEN
349 ELSE IF( minmn.GT.0 )
THEN
351 IF( vl.LT.zero )
THEN
353 ELSE IF( vu.LE.vl )
THEN
357 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) )
THEN
359 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn )
THEN
364 IF( wantu .AND. ldu.LT.m )
THEN
366 ELSE IF( wantvt )
THEN
368 IF( ldvt.LT.iu-il+1 )
THEN
371 ELSE IF( ldvt.LT.minmn )
THEN
388 IF( minmn.GT.0 )
THEN
390 mnthr = ilaenv( 6,
'DGESVD', jobu // jobvt, m, n, 0, 0 )
391 IF( m.GE.mnthr )
THEN
396 $ n*ilaenv( 1,
'DGEQRF',
' ', m, n, -1, -1 )
397 maxwrk = max( maxwrk, n*(n+5) + 2*n*
398 $ ilaenv( 1,
'DGEBRD',
' ', n, n, -1, -1 ) )
400 maxwrk = max(maxwrk,n*(n*3+6)+n*
401 $ ilaenv( 1,
'DORMQR',
' ', n, n, -1, -1 ) )
404 maxwrk = max(maxwrk,n*(n*3+6)+n*
405 $ ilaenv( 1,
'DORMLQ',
' ', n, n, -1, -1 ) )
412 maxwrk = 4*n + ( m+n )*
413 $ ilaenv( 1,
'DGEBRD',
' ', m, n, -1, -1 )
415 maxwrk = max(maxwrk,n*(n*2+5)+n*
416 $ ilaenv( 1,
'DORMQR',
' ', n, n, -1, -1 ) )
419 maxwrk = max(maxwrk,n*(n*2+5)+n*
420 $ ilaenv( 1,
'DORMLQ',
' ', n, n, -1, -1 ) )
422 minwrk = max(n*(n*2+19),4*n+m)
425 mnthr = ilaenv( 6,
'DGESVD', jobu // jobvt, m, n, 0, 0 )
426 IF( n.GE.mnthr )
THEN
431 $ m*ilaenv( 1,
'DGELQF',
' ', m, n, -1, -1 )
432 maxwrk = max( maxwrk, m*(m+5) + 2*m*
433 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
435 maxwrk = max(maxwrk,m*(m*3+6)+m*
436 $ ilaenv( 1,
'DORMQR',
' ', m, m, -1, -1 ) )
439 maxwrk = max(maxwrk,m*(m*3+6)+m*
440 $ ilaenv( 1,
'DORMLQ',
' ', m, m, -1, -1 ) )
447 maxwrk = 4*m + ( m+n )*
448 $ ilaenv( 1,
'DGEBRD',
' ', m, n, -1, -1 )
450 maxwrk = max(maxwrk,m*(m*2+5)+m*
451 $ ilaenv( 1,
'DORMQR',
' ', m, m, -1, -1 ) )
454 maxwrk = max(maxwrk,m*(m*2+5)+m*
455 $ ilaenv( 1,
'DORMLQ',
' ', m, m, -1, -1 ) )
457 minwrk = max(m*(m*2+19),4*m+n)
461 maxwrk = max( maxwrk, minwrk )
462 work( 1 ) = dble( maxwrk )
464 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
470 CALL xerbla(
'DGESVDX', -info )
472 ELSE IF( lquery )
THEN
478 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
501 smlnum = sqrt( dlamch(
'S' ) ) / eps
502 bignum = one / smlnum
506 anrm = dlange(
'M', m, n, a, lda, dum )
508 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
510 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
511 ELSE IF( anrm.GT.bignum )
THEN
513 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
522 IF( m.GE.mnthr )
THEN
534 CALL dgeqrf( m, n, a, lda, work( itau ), work( itemp ),
535 $ lwork-itemp+1, info )
546 CALL dlacpy(
'U', n, n, a, lda, work( iqrf ), n )
547 CALL dlaset(
'L', n-1, n-1, zero, zero, work( iqrf+1 ), n )
548 CALL dgebrd( n, n, work( iqrf ), n, work( id ), work( ie ),
549 $ work( itauq ), work( itaup ), work( itemp ),
550 $ lwork-itemp+1, info )
556 itemp = itgkz + n*(n*2+1)
557 CALL dbdsvdx(
'U', jobz, rngtgk, n, work( id ), work( ie ),
558 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
559 $ n*2, work( itemp ), iwork, info)
566 CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
569 CALL dlaset(
'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
574 CALL dormbr(
'Q',
'L',
'N', n, ns, n, work( iqrf ), n,
575 $ work( itauq ), u, ldu, work( itemp ),
576 $ lwork-itemp+1, info )
581 CALL dormqr(
'L',
'N', m, ns, n, a, lda,
582 $ work( itau ), u, ldu, work( itemp ),
583 $ lwork-itemp+1, info )
591 CALL dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
598 CALL dormbr(
'P',
'R',
'T', ns, n, n, work( iqrf ), n,
599 $ work( itaup ), vt, ldvt, work( itemp ),
600 $ lwork-itemp+1, info )
617 CALL dgebrd( m, n, a, lda, work( id ), work( ie ),
618 $ work( itauq ), work( itaup ), work( itemp ),
619 $ lwork-itemp+1, info )
625 itemp = itgkz + n*(n*2+1)
626 CALL dbdsvdx(
'U', jobz, rngtgk, n, work( id ), work( ie ),
627 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
628 $ n*2, work( itemp ), iwork, info)
635 CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
638 CALL dlaset(
'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
643 CALL dormbr(
'Q',
'L',
'N', m, ns, n, a, lda,
644 $ work( itauq ), u, ldu, work( itemp ),
645 $ lwork-itemp+1, ierr )
653 CALL dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
660 CALL dormbr(
'P',
'R',
'T', ns, n, n, a, lda,
661 $ work( itaup ), vt, ldvt, work( itemp ),
662 $ lwork-itemp+1, ierr )
670 IF( n.GE.mnthr )
THEN
682 CALL dgelqf( m, n, a, lda, work( itau ), work( itemp ),
683 $ lwork-itemp+1, info )
694 CALL dlacpy(
'L', m, m, a, lda, work( ilqf ), m )
695 CALL dlaset(
'U', m-1, m-1, zero, zero, work( ilqf+m ), m )
696 CALL dgebrd( m, m, work( ilqf ), m, work( id ), work( ie ),
697 $ work( itauq ), work( itaup ), work( itemp ),
698 $ lwork-itemp+1, info )
704 itemp = itgkz + m*(m*2+1)
705 CALL dbdsvdx(
'U', jobz, rngtgk, m, work( id ), work( ie ),
706 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
707 $ m*2, work( itemp ), iwork, info)
714 CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
721 CALL dormbr(
'Q',
'L',
'N', m, ns, m, work( ilqf ), m,
722 $ work( itauq ), u, ldu, work( itemp ),
723 $ lwork-itemp+1, info )
731 CALL dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
734 CALL dlaset(
'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
739 CALL dormbr(
'P',
'R',
'T', ns, m, m, work( ilqf ), m,
740 $ work( itaup ), vt, ldvt, work( itemp ),
741 $ lwork-itemp+1, info )
746 CALL dormlq(
'R',
'N', ns, n, m, a, lda,
747 $ work( itau ), vt, ldvt, work( itemp ),
748 $ lwork-itemp+1, info )
765 CALL dgebrd( m, n, a, lda, work( id ), work( ie ),
766 $ work( itauq ), work( itaup ), work( itemp ),
767 $ lwork-itemp+1, info )
773 itemp = itgkz + m*(m*2+1)
774 CALL dbdsvdx(
'L', jobz, rngtgk, m, work( id ), work( ie ),
775 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
776 $ m*2, work( itemp ), iwork, info)
783 CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
790 CALL dormbr(
'Q',
'L',
'N', m, ns, n, a, lda,
791 $ work( itauq ), u, ldu, work( itemp ),
792 $ lwork-itemp+1, info )
800 CALL dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
803 CALL dlaset(
'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
808 CALL dormbr(
'P',
'R',
'T', ns, n, m, a, lda,
809 $ work( itaup ), vt, ldvt, work( itemp ),
810 $ lwork-itemp+1, info )
819 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1,
822 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1,
828 work( 1 ) = dble( maxwrk )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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 dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
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 dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dgesvdx(JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, IL, IU, NS, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESVDX computes the singular value decomposition (SVD) for GE matrices
subroutine dbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
DBDSVDX