258 SUBROUTINE dgesvdx( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
259 $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
260 $ LWORK, IWORK, INFO )
267 CHARACTER JOBU, JOBVT, RANGE
268 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
269 DOUBLE PRECISION VL, VU
273 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
274 $ vt( ldvt, * ), work( * )
280 DOUBLE PRECISION ZERO, ONE
281 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
284 CHARACTER JOBZ, RNGTGK
285 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
286 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
287 $ itau, itaup, itauq, itemp, itgkz, iutgk,
288 $ j, maxwrk, minmn, minwrk, mnthr
289 DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
292 DOUBLE PRECISION DUM( 1 )
303 DOUBLE PRECISION DLAMCH, DLANGE
304 EXTERNAL lsame, ilaenv, dlamch, dlange
307 INTRINSIC max, min, sqrt
315 abstol = 2*dlamch(
'S')
316 lquery = ( lwork.EQ.-1 )
319 wantu = lsame( jobu,
'V' )
320 wantvt = lsame( jobvt,
'V' )
321 IF( wantu .OR. wantvt )
THEN
326 alls = lsame( range,
'A' )
327 vals = lsame( range,
'V' )
328 inds = lsame( range,
'I' )
331 IF( .NOT.lsame( jobu,
'V' ) .AND.
332 $ .NOT.lsame( jobu,
'N' ) )
THEN
334 ELSE IF( .NOT.lsame( jobvt,
'V' ) .AND.
335 $ .NOT.lsame( jobvt,
'N' ) )
THEN
337 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) )
THEN
339 ELSE IF( m.LT.0 )
THEN
341 ELSE IF( n.LT.0 )
THEN
343 ELSE IF( m.GT.lda )
THEN
345 ELSE IF( minmn.GT.0 )
THEN
347 IF( vl.LT.zero )
THEN
349 ELSE IF( vu.LE.vl )
THEN
353 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) )
THEN
355 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn )
THEN
360 IF( wantu .AND. ldu.LT.m )
THEN
362 ELSE IF( wantvt )
THEN
364 IF( ldvt.LT.iu-il+1 )
THEN
367 ELSE IF( ldvt.LT.minmn )
THEN
384 IF( minmn.GT.0 )
THEN
386 mnthr = ilaenv( 6,
'DGESVD', jobu // jobvt, m, n, 0,
388 IF( m.GE.mnthr )
THEN
393 $ n*ilaenv( 1,
'DGEQRF',
' ', m, n, -1, -1 )
394 maxwrk = max( maxwrk, n*(n+5) + 2*n*
395 $ ilaenv( 1,
'DGEBRD',
' ', n, n, -1, -1 ) )
397 maxwrk = max(maxwrk,n*(n*3+6)+n*
398 $ ilaenv( 1,
'DORMQR',
' ', n, n, -1, -1 ) )
401 maxwrk = max(maxwrk,n*(n*3+6)+n*
402 $ ilaenv( 1,
'DORMLQ',
' ', n, n, -1, -1 ) )
409 maxwrk = 4*n + ( m+n )*
410 $ ilaenv( 1,
'DGEBRD',
' ', m, n, -1, -1 )
412 maxwrk = max(maxwrk,n*(n*2+5)+n*
413 $ ilaenv( 1,
'DORMQR',
' ', n, n, -1, -1 ) )
416 maxwrk = max(maxwrk,n*(n*2+5)+n*
417 $ ilaenv( 1,
'DORMLQ',
' ', n, n, -1, -1 ) )
419 minwrk = max(n*(n*2+19),4*n+m)
422 mnthr = ilaenv( 6,
'DGESVD', jobu // jobvt, m, n, 0,
424 IF( n.GE.mnthr )
THEN
429 $ m*ilaenv( 1,
'DGELQF',
' ', m, n, -1, -1 )
430 maxwrk = max( maxwrk, m*(m+5) + 2*m*
431 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
433 maxwrk = max(maxwrk,m*(m*3+6)+m*
434 $ ilaenv( 1,
'DORMQR',
' ', m, m, -1, -1 ) )
437 maxwrk = max(maxwrk,m*(m*3+6)+m*
438 $ ilaenv( 1,
'DORMLQ',
' ', m, m, -1, -1 ) )
445 maxwrk = 4*m + ( m+n )*
446 $ ilaenv( 1,
'DGEBRD',
' ', m, n, -1, -1 )
448 maxwrk = max(maxwrk,m*(m*2+5)+m*
449 $ ilaenv( 1,
'DORMQR',
' ', m, m, -1, -1 ) )
452 maxwrk = max(maxwrk,m*(m*2+5)+m*
453 $ ilaenv( 1,
'DORMLQ',
' ', m, m, -1, -1 ) )
455 minwrk = max(m*(m*2+19),4*m+n)
459 maxwrk = max( maxwrk, minwrk )
460 work( 1 ) = dble( maxwrk )
462 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
468 CALL xerbla(
'DGESVDX', -info )
470 ELSE IF( lquery )
THEN
476 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
499 smlnum = sqrt( dlamch(
'S' ) ) / eps
500 bignum = one / smlnum
504 anrm = dlange(
'M', m, n, a, lda, dum )
506 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
508 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
509 ELSE IF( anrm.GT.bignum )
THEN
511 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
520 IF( m.GE.mnthr )
THEN
532 CALL dgeqrf( m, n, a, lda, work( itau ), work( itemp ),
533 $ lwork-itemp+1, info )
544 CALL dlacpy(
'U', n, n, a, lda, work( iqrf ), n )
545 CALL dlaset(
'L', n-1, n-1, zero, zero, work( iqrf+1 ),
547 CALL dgebrd( n, n, work( iqrf ), n, work( id ),
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 ),
559 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
560 $ n*2, work( itemp ), iwork, info)
567 CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
570 CALL dlaset(
'A', m-n, ns, zero, zero, u( n+1,1 ),
576 CALL dormbr(
'Q',
'L',
'N', n, ns, n, work( iqrf ), n,
577 $ work( itauq ), u, ldu, work( itemp ),
578 $ lwork-itemp+1, info )
583 CALL dormqr(
'L',
'N', m, ns, n, a, lda,
584 $ work( itau ), u, ldu, work( itemp ),
585 $ lwork-itemp+1, info )
593 CALL dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
600 CALL dormbr(
'P',
'R',
'T', ns, n, n, work( iqrf ), n,
601 $ work( itaup ), vt, ldvt, work( itemp ),
602 $ lwork-itemp+1, info )
619 CALL dgebrd( m, n, a, lda, work( id ), work( ie ),
620 $ work( itauq ), work( itaup ), work( itemp ),
621 $ lwork-itemp+1, info )
627 itemp = itgkz + n*(n*2+1)
628 CALL dbdsvdx(
'U', jobz, rngtgk, n, work( id ),
630 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
631 $ n*2, work( itemp ), iwork, info)
638 CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
641 CALL dlaset(
'A', m-n, ns, zero, zero, u( n+1,1 ),
647 CALL dormbr(
'Q',
'L',
'N', m, ns, n, a, lda,
648 $ work( itauq ), u, ldu, work( itemp ),
649 $ lwork-itemp+1, ierr )
657 CALL dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
664 CALL dormbr(
'P',
'R',
'T', ns, n, n, a, lda,
665 $ work( itaup ), vt, ldvt, work( itemp ),
666 $ lwork-itemp+1, ierr )
674 IF( n.GE.mnthr )
THEN
686 CALL dgelqf( m, n, a, lda, work( itau ), work( itemp ),
687 $ lwork-itemp+1, info )
698 CALL dlacpy(
'L', m, m, a, lda, work( ilqf ), m )
699 CALL dlaset(
'U', m-1, m-1, zero, zero, work( ilqf+m ),
701 CALL dgebrd( m, m, work( ilqf ), m, work( id ),
703 $ work( itauq ), work( itaup ), work( itemp ),
704 $ lwork-itemp+1, info )
710 itemp = itgkz + m*(m*2+1)
711 CALL dbdsvdx(
'U', jobz, rngtgk, m, work( id ),
713 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
714 $ m*2, work( itemp ), iwork, info)
721 CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
728 CALL dormbr(
'Q',
'L',
'N', m, ns, m, work( ilqf ), m,
729 $ work( itauq ), u, ldu, work( itemp ),
730 $ lwork-itemp+1, info )
738 CALL dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
741 CALL dlaset(
'A', ns, n-m, zero, zero, vt( 1,m+1 ),
747 CALL dormbr(
'P',
'R',
'T', ns, m, m, work( ilqf ), m,
748 $ work( itaup ), vt, ldvt, work( itemp ),
749 $ lwork-itemp+1, info )
754 CALL dormlq(
'R',
'N', ns, n, m, a, lda,
755 $ work( itau ), vt, ldvt, work( itemp ),
756 $ lwork-itemp+1, info )
773 CALL dgebrd( m, n, a, lda, work( id ), work( ie ),
774 $ work( itauq ), work( itaup ), work( itemp ),
775 $ lwork-itemp+1, info )
781 itemp = itgkz + m*(m*2+1)
782 CALL dbdsvdx(
'L', jobz, rngtgk, m, work( id ),
784 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
785 $ m*2, work( itemp ), iwork, info)
792 CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
799 CALL dormbr(
'Q',
'L',
'N', m, ns, n, a, lda,
800 $ work( itauq ), u, ldu, work( itemp ),
801 $ lwork-itemp+1, info )
809 CALL dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
812 CALL dlaset(
'A', ns, n-m, zero, zero, vt( 1,m+1 ),
818 CALL dormbr(
'P',
'R',
'T', ns, n, m, a, lda,
819 $ work( itaup ), vt, ldvt, work( itemp ),
820 $ lwork-itemp+1, info )
829 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1,
832 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1,
838 work( 1 ) = dble( maxwrk )