260 SUBROUTINE sgesvdx( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
261 $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
262 $ LWORK, IWORK, INFO )
269 CHARACTER JOBU, JOBVT, RANGE
270 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
275 REAL A( LDA, * ), S( * ), U( LDU, * ),
276 $ vt( ldvt, * ), work( * )
283 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
286 CHARACTER JOBZ, RNGTGK
287 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
288 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
289 $ itau, itaup, itauq, itemp, itgkz, iutgk,
290 $ j, maxwrk, minmn, minwrk, mnthr
291 REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
304 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
305 EXTERNAL lsame, ilaenv, slamch, slange, sroundup_lwork
308 INTRINSIC max, min, sqrt
316 abstol = 2*slamch(
'S')
317 lquery = ( lwork.EQ.-1 )
320 wantu = lsame( jobu,
'V' )
321 wantvt = lsame( jobvt,
'V' )
322 IF( wantu .OR. wantvt )
THEN
327 alls = lsame( range,
'A' )
328 vals = lsame( range,
'V' )
329 inds = lsame( range,
'I' )
332 IF( .NOT.lsame( jobu,
'V' ) .AND.
333 $ .NOT.lsame( jobu,
'N' ) )
THEN
335 ELSE IF( .NOT.lsame( jobvt,
'V' ) .AND.
336 $ .NOT.lsame( jobvt,
'N' ) )
THEN
338 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) )
THEN
340 ELSE IF( m.LT.0 )
THEN
342 ELSE IF( n.LT.0 )
THEN
344 ELSE IF( m.GT.lda )
THEN
346 ELSE IF( minmn.GT.0 )
THEN
348 IF( vl.LT.zero )
THEN
350 ELSE IF( vu.LE.vl )
THEN
354 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) )
THEN
356 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn )
THEN
361 IF( wantu .AND. ldu.LT.m )
THEN
363 ELSE IF( wantvt )
THEN
365 IF( ldvt.LT.iu-il+1 )
THEN
368 ELSE IF( ldvt.LT.minmn )
THEN
385 IF( minmn.GT.0 )
THEN
387 mnthr = ilaenv( 6,
'SGESVD', jobu // jobvt, m, n, 0, 0 )
388 IF( m.GE.mnthr )
THEN
393 $ n*ilaenv( 1,
'SGEQRF',
' ', m, n, -1, -1 )
394 maxwrk = max( maxwrk, n*(n+5) + 2*n*
395 $ ilaenv( 1,
'SGEBRD',
' ', n, n, -1, -1 ) )
397 maxwrk = max(maxwrk,n*(n*3+6)+n*
398 $ ilaenv( 1,
'SORMQR',
' ', n, n, -1, -1 ) )
401 maxwrk = max(maxwrk,n*(n*3+6)+n*
402 $ ilaenv( 1,
'SORMLQ',
' ', n, n, -1, -1 ) )
409 maxwrk = 4*n + ( m+n )*
410 $ ilaenv( 1,
'SGEBRD',
' ', m, n, -1, -1 )
412 maxwrk = max(maxwrk,n*(n*2+5)+n*
413 $ ilaenv( 1,
'SORMQR',
' ', n, n, -1, -1 ) )
416 maxwrk = max(maxwrk,n*(n*2+5)+n*
417 $ ilaenv( 1,
'SORMLQ',
' ', n, n, -1, -1 ) )
419 minwrk = max(n*(n*2+19),4*n+m)
422 mnthr = ilaenv( 6,
'SGESVD', jobu // jobvt, m, n, 0, 0 )
423 IF( n.GE.mnthr )
THEN
428 $ m*ilaenv( 1,
'SGELQF',
' ', m, n, -1, -1 )
429 maxwrk = max( maxwrk, m*(m+5) + 2*m*
430 $ ilaenv( 1,
'SGEBRD',
' ', m, m, -1, -1 ) )
432 maxwrk = max(maxwrk,m*(m*3+6)+m*
433 $ ilaenv( 1,
'SORMQR',
' ', m, m, -1, -1 ) )
436 maxwrk = max(maxwrk,m*(m*3+6)+m*
437 $ ilaenv( 1,
'SORMLQ',
' ', m, m, -1, -1 ) )
444 maxwrk = 4*m + ( m+n )*
445 $ ilaenv( 1,
'SGEBRD',
' ', m, n, -1, -1 )
447 maxwrk = max(maxwrk,m*(m*2+5)+m*
448 $ ilaenv( 1,
'SORMQR',
' ', m, m, -1, -1 ) )
451 maxwrk = max(maxwrk,m*(m*2+5)+m*
452 $ ilaenv( 1,
'SORMLQ',
' ', m, m, -1, -1 ) )
454 minwrk = max(m*(m*2+19),4*m+n)
458 maxwrk = max( maxwrk, minwrk )
459 work( 1 ) = sroundup_lwork( maxwrk )
461 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
467 CALL xerbla(
'SGESVDX', -info )
469 ELSE IF( lquery )
THEN
475 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
498 smlnum = sqrt( slamch(
'S' ) ) / eps
499 bignum = one / smlnum
503 anrm = slange(
'M', m, n, a, lda, dum )
505 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
507 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
508 ELSE IF( anrm.GT.bignum )
THEN
510 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
519 IF( m.GE.mnthr )
THEN
531 CALL sgeqrf( m, n, a, lda, work( itau ), work( itemp ),
532 $ lwork-itemp+1, info )
543 CALL slacpy(
'U', n, n, a, lda, work( iqrf ), n )
544 CALL slaset(
'L', n-1, n-1, zero, zero, work( iqrf+1 ), n )
545 CALL sgebrd( n, n, work( iqrf ), n, work( id ), work( ie ),
546 $ work( itauq ), work( itaup ), work( itemp ),
547 $ lwork-itemp+1, info )
553 itemp = itgkz + n*(n*2+1)
554 CALL sbdsvdx(
'U', jobz, rngtgk, n, work( id ), work( ie ),
555 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
556 $ n*2, work( itemp ), iwork, info)
563 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
566 CALL slaset(
'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
571 CALL sormbr(
'Q',
'L',
'N', n, ns, n, work( iqrf ), n,
572 $ work( itauq ), u, ldu, work( itemp ),
573 $ lwork-itemp+1, info )
578 CALL sormqr(
'L',
'N', m, ns, n, a, lda,
579 $ work( itau ), u, ldu, work( itemp ),
580 $ lwork-itemp+1, info )
588 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
595 CALL sormbr(
'P',
'R',
'T', ns, n, n, work( iqrf ), n,
596 $ work( itaup ), vt, ldvt, work( itemp ),
597 $ lwork-itemp+1, info )
614 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
615 $ work( itauq ), work( itaup ), work( itemp ),
616 $ lwork-itemp+1, info )
622 itemp = itgkz + n*(n*2+1)
623 CALL sbdsvdx(
'U', jobz, rngtgk, n, work( id ), work( ie ),
624 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
625 $ n*2, work( itemp ), iwork, info)
632 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
635 CALL slaset(
'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
640 CALL sormbr(
'Q',
'L',
'N', m, ns, n, a, lda,
641 $ work( itauq ), u, ldu, work( itemp ),
642 $ lwork-itemp+1, ierr )
650 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
657 CALL sormbr(
'P',
'R',
'T', ns, n, n, a, lda,
658 $ work( itaup ), vt, ldvt, work( itemp ),
659 $ lwork-itemp+1, ierr )
667 IF( n.GE.mnthr )
THEN
679 CALL sgelqf( m, n, a, lda, work( itau ), work( itemp ),
680 $ lwork-itemp+1, info )
691 CALL slacpy(
'L', m, m, a, lda, work( ilqf ), m )
692 CALL slaset(
'U', m-1, m-1, zero, zero, work( ilqf+m ), m )
693 CALL sgebrd( m, m, work( ilqf ), m, work( id ), work( ie ),
694 $ work( itauq ), work( itaup ), work( itemp ),
695 $ lwork-itemp+1, info )
701 itemp = itgkz + m*(m*2+1)
702 CALL sbdsvdx(
'U', jobz, rngtgk, m, work( id ), work( ie ),
703 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
704 $ m*2, work( itemp ), iwork, info)
711 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
718 CALL sormbr(
'Q',
'L',
'N', m, ns, m, work( ilqf ), m,
719 $ work( itauq ), u, ldu, work( itemp ),
720 $ lwork-itemp+1, info )
728 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
731 CALL slaset(
'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
736 CALL sormbr(
'P',
'R',
'T', ns, m, m, work( ilqf ), m,
737 $ work( itaup ), vt, ldvt, work( itemp ),
738 $ lwork-itemp+1, info )
743 CALL sormlq(
'R',
'N', ns, n, m, a, lda,
744 $ work( itau ), vt, ldvt, work( itemp ),
745 $ lwork-itemp+1, info )
762 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
763 $ work( itauq ), work( itaup ), work( itemp ),
764 $ lwork-itemp+1, info )
770 itemp = itgkz + m*(m*2+1)
771 CALL sbdsvdx(
'L', jobz, rngtgk, m, work( id ), work( ie ),
772 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
773 $ m*2, work( itemp ), iwork, info)
780 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
787 CALL sormbr(
'Q',
'L',
'N', m, ns, n, a, lda,
788 $ work( itauq ), u, ldu, work( itemp ),
789 $ lwork-itemp+1, info )
797 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
800 CALL slaset(
'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
805 CALL sormbr(
'P',
'R',
'T', ns, n, m, a, lda,
806 $ work( itaup ), vt, ldvt, work( itemp ),
807 $ lwork-itemp+1, info )
816 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1,
819 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1,
825 work( 1 ) = sroundup_lwork( maxwrk )
subroutine xerbla(srname, info)
subroutine sbdsvdx(uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
SBDSVDX
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
subroutine sgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
SGEBRD
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
subroutine sgesvdx(jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
SGESVDX computes the singular value decomposition (SVD) for GE matrices
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
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 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 sormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMBR
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR