LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine sgesvdx ( character  JOBU,
character  JOBVT,
character  RANGE,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real  VL,
real  VU,
integer  IL,
integer  IU,
integer  NS,
real, dimension( * )  S,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldvt, * )  VT,
integer  LDVT,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SGESVDX computes the singular value decomposition (SVD) for GE matrices

Download SGESVDX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
  SGESVDX computes the singular value decomposition (SVD) of a real
  M-by-N matrix A, optionally computing the left and/or right singular
  vectors. The SVD is written
 
      A = U * SIGMA * transpose(V)
 
  where SIGMA is an M-by-N matrix which is zero except for its
  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  are the singular values of A; they are real and non-negative, and
  are returned in descending order.  The first min(m,n) columns of
  U and V are the left and right singular vectors of A.
 
  SGESVDX uses an eigenvalue problem for obtaining the SVD, which 
  allows for the computation of a subset of singular values and 
  vectors. See SBDSVDX for details.
 
  Note that the routine returns V**T, not V.
Parameters
[in]JOBU
          JOBU is CHARACTER*1
          Specifies options for computing all or part of the matrix U:
          = 'V':  the first min(m,n) columns of U (the left singular
                  vectors) or as specified by RANGE are returned in 
                  the array U;
          = 'N':  no columns of U (no left singular vectors) are
                  computed.
[in]JOBVT
          JOBVT is CHARACTER*1
           Specifies options for computing all or part of the matrix
           V**T:
           = 'V':  the first min(m,n) rows of V**T (the right singular
                   vectors) or as specified by RANGE are returned in 
                   the array VT;
           = 'N':  no rows of V**T (no right singular vectors) are
                   computed.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all singular values will be found.
          = 'V': all singular values in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th singular values will be found. 
[in]M
          M is INTEGER
          The number of rows of the input matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the input matrix A.  N >= 0.
[in,out]A
          A is REAL array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the contents of A are destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in]VL
          VL is REAL
          If RANGE='V', the lower bound of the interval to
          be searched for singular values. VU > VL.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
          VU is REAL
          If RANGE='V', the upper bound of the interval to
          be searched for singular values. VU > VL.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest singular value to be returned.
          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest singular value to be returned.
          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
          Not referenced if RANGE = 'A' or 'V'.
[out]NS
          NS is INTEGER
          The total number of singular values found,  
          0 <= NS <= min(M,N).
          If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
[out]S
          S is REAL array, dimension (min(M,N))
          The singular values of A, sorted so that S(i) >= S(i+1).
[out]U
          U is REAL array, dimension (LDU,UCOL)
          If JOBU = 'V', U contains columns of U (the left singular 
          vectors, stored columnwise) as specified by RANGE; if 
          JOBU = 'N', U is not referenced.
          Note: The user must ensure that UCOL >= NS; if RANGE = 'V', 
          the exact value of NS is not known in advance and an upper
          bound must be used.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= 1; if
          JOBU = 'V', LDU >= M.
[out]VT
          VT is REAL array, dimension (LDVT,N)
          If JOBVT = 'V', VT contains the rows of V**T (the right singular 
          vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N', 
          VT is not referenced.
          Note: The user must ensure that LDVT >= NS; if RANGE = 'V', 
          the exact value of NS is not known in advance and an upper 
          bound must be used.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.  LDVT >= 1; if
          JOBVT = 'V', LDVT >= NS (see above).
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see 
          comments inside the code):
             - PATH 1  (M much larger than N) 
             - PATH 1t (N much larger than M)
          LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (12*MIN(M,N))
          If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0, 
          then IWORK contains the indices of the eigenvectors that failed 
          to converge in SBDSVDX/SSTEVX.
[out]INFO
     INFO is INTEGER
           = 0:  successful exit
           < 0:  if INFO = -i, the i-th argument had an illegal value
           > 0:  if INFO = i, then i eigenvectors failed to converge
                 in SBDSVDX/SSTEVX.
                 if INFO = N*2 + 1, an internal error occurred in
                 SBDSVDX
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 265 of file sgesvdx.f.

265 *
266 * -- LAPACK driver routine (version 3.6.1) --
267 * -- LAPACK is a software package provided by Univ. of Tennessee, --
268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269 * June 2016
270 *
271 * .. Scalar Arguments ..
272  CHARACTER jobu, jobvt, range
273  INTEGER il, info, iu, lda, ldu, ldvt, lwork, m, n, ns
274  REAL vl, vu
275 * ..
276 * .. Array Arguments ..
277  INTEGER iwork( * )
278  REAL a( lda, * ), s( * ), u( ldu, * ),
279  $ vt( ldvt, * ), work( * )
280 * ..
281 *
282 * =====================================================================
283 *
284 * .. Parameters ..
285  REAL zero, one
286  parameter ( zero = 0.0e0, one = 1.0e0 )
287 * ..
288 * .. Local Scalars ..
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  REAL abstol, anrm, bignum, eps, smlnum
295 * ..
296 * .. Local Arrays ..
297  REAL dum( 1 )
298 * ..
299 * .. External Subroutines ..
300  EXTERNAL sbdsvdx, sgebrd, sgelqf, sgeqrf, slacpy,
302  $ sscal, xerbla
303 * ..
304 * .. External Functions ..
305  LOGICAL lsame
306  INTEGER ilaenv
307  REAL slamch, slange, snrm2
308  EXTERNAL lsame, ilaenv, slamch, slange, snrm2
309 * ..
310 * .. Intrinsic Functions ..
311  INTRINSIC max, min, sqrt
312 * ..
313 * .. Executable Statements ..
314 *
315 * Test the input arguments.
316 *
317  ns = 0
318  info = 0
319  abstol = 2*slamch('S')
320  lquery = ( lwork.EQ.-1 )
321  minmn = min( m, n )
322 
323  wantu = lsame( jobu, 'V' )
324  wantvt = lsame( jobvt, 'V' )
325  IF( wantu .OR. wantvt ) THEN
326  jobz = 'V'
327  ELSE
328  jobz = 'N'
329  END IF
330  alls = lsame( range, 'A' )
331  vals = lsame( range, 'V' )
332  inds = lsame( range, 'I' )
333 *
334  info = 0
335  IF( .NOT.lsame( jobu, 'V' ) .AND.
336  $ .NOT.lsame( jobu, 'N' ) ) THEN
337  info = -1
338  ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
339  $ .NOT.lsame( jobvt, 'N' ) ) THEN
340  info = -2
341  ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
342  info = -3
343  ELSE IF( m.LT.0 ) THEN
344  info = -4
345  ELSE IF( n.LT.0 ) THEN
346  info = -5
347  ELSE IF( m.GT.lda ) THEN
348  info = -7
349  ELSE IF( minmn.GT.0 ) THEN
350  IF( vals ) THEN
351  IF( vl.LT.zero ) THEN
352  info = -8
353  ELSE IF( vu.LE.vl ) THEN
354  info = -9
355  END IF
356  ELSE IF( inds ) THEN
357  IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
358  info = -10
359  ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
360  info = -11
361  END IF
362  END IF
363  IF( info.EQ.0 ) THEN
364  IF( wantu .AND. ldu.LT.m ) THEN
365  info = -15
366  ELSE IF( wantvt ) THEN
367  IF( inds ) THEN
368  IF( ldvt.LT.iu-il+1 ) THEN
369  info = -17
370  END IF
371  ELSE IF( ldvt.LT.minmn ) THEN
372  info = -17
373  END IF
374  END IF
375  END IF
376  END IF
377 *
378 * Compute workspace
379 * (Note: Comments in the code beginning "Workspace:" describe the
380 * minimal amount of workspace needed at that point in the code,
381 * as well as the preferred amount for good performance.
382 * NB refers to the optimal block size for the immediately
383 * following subroutine, as returned by ILAENV.)
384 *
385  IF( info.EQ.0 ) THEN
386  minwrk = 1
387  maxwrk = 1
388  IF( minmn.GT.0 ) THEN
389  IF( m.GE.n ) THEN
390  mnthr = ilaenv( 6, 'SGESVD', jobu // jobvt, m, n, 0, 0 )
391  IF( m.GE.mnthr ) THEN
392 *
393 * Path 1 (M much larger than N)
394 *
395  maxwrk = n +
396  $ n*ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
397  maxwrk = max( maxwrk, n*(n+5) + 2*n*
398  $ ilaenv( 1, 'SGEBRD', ' ', n, n, -1, -1 ) )
399  IF (wantu) THEN
400  maxwrk = max(maxwrk,n*(n*3+6)+n*
401  $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
402  END IF
403  IF (wantvt) THEN
404  maxwrk = max(maxwrk,n*(n*3+6)+n*
405  $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
406  END IF
407  minwrk = n*(n*3+20)
408  ELSE
409 *
410 * Path 2 (M at least N, but not much larger)
411 *
412  maxwrk = 4*n + ( m+n )*
413  $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
414  IF (wantu) THEN
415  maxwrk = max(maxwrk,n*(n*2+5)+n*
416  $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
417  END IF
418  IF (wantvt) THEN
419  maxwrk = max(maxwrk,n*(n*2+5)+n*
420  $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
421  END IF
422  minwrk = max(n*(n*2+19),4*n+m)
423  END IF
424  ELSE
425  mnthr = ilaenv( 6, 'SGESVD', jobu // jobvt, m, n, 0, 0 )
426  IF( n.GE.mnthr ) THEN
427 *
428 * Path 1t (N much larger than M)
429 *
430  maxwrk = m +
431  $ m*ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
432  maxwrk = max( maxwrk, m*(m+5) + 2*m*
433  $ ilaenv( 1, 'SGEBRD', ' ', m, m, -1, -1 ) )
434  IF (wantu) THEN
435  maxwrk = max(maxwrk,m*(m*3+6)+m*
436  $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
437  END IF
438  IF (wantvt) THEN
439  maxwrk = max(maxwrk,m*(m*3+6)+m*
440  $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
441  END IF
442  minwrk = m*(m*3+20)
443  ELSE
444 *
445 * Path 2t (N at least M, but not much larger)
446 *
447  maxwrk = 4*m + ( m+n )*
448  $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
449  IF (wantu) THEN
450  maxwrk = max(maxwrk,m*(m*2+5)+m*
451  $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
452  END IF
453  IF (wantvt) THEN
454  maxwrk = max(maxwrk,m*(m*2+5)+m*
455  $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
456  END IF
457  minwrk = max(m*(m*2+19),4*m+n)
458  END IF
459  END IF
460  END IF
461  maxwrk = max( maxwrk, minwrk )
462  work( 1 ) = REAL( maxwrk )
463 *
464  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
465  info = -19
466  END IF
467  END IF
468 *
469  IF( info.NE.0 ) THEN
470  CALL xerbla( 'SGESVDX', -info )
471  RETURN
472  ELSE IF( lquery ) THEN
473  RETURN
474  END IF
475 *
476 * Quick return if possible
477 *
478  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
479  RETURN
480  END IF
481 *
482 * Set singular values indices accord to RANGE.
483 *
484  IF( alls ) THEN
485  rngtgk = 'I'
486  iltgk = 1
487  iutgk = min( m, n )
488  ELSE IF( inds ) THEN
489  rngtgk = 'I'
490  iltgk = il
491  iutgk = iu
492  ELSE
493  rngtgk = 'V'
494  iltgk = 0
495  iutgk = 0
496  END IF
497 *
498 * Get machine constants
499 *
500  eps = slamch( 'P' )
501  smlnum = sqrt( slamch( 'S' ) ) / eps
502  bignum = one / smlnum
503 *
504 * Scale A if max element outside range [SMLNUM,BIGNUM]
505 *
506  anrm = slange( 'M', m, n, a, lda, dum )
507  iscl = 0
508  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
509  iscl = 1
510  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
511  ELSE IF( anrm.GT.bignum ) THEN
512  iscl = 1
513  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
514  END IF
515 *
516  IF( m.GE.n ) THEN
517 *
518 * A has at least as many rows as columns. If A has sufficiently
519 * more rows than columns, first reduce A using the QR
520 * decomposition.
521 *
522  IF( m.GE.mnthr ) THEN
523 *
524 * Path 1 (M much larger than N):
525 * A = Q * R = Q * ( QB * B * PB**T )
526 * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
527 * U = Q * QB * UB; V**T = VB**T * PB**T
528 *
529 * Compute A=Q*R
530 * (Workspace: need 2*N, prefer N+N*NB)
531 *
532  itau = 1
533  itemp = itau + n
534  CALL sgeqrf( m, n, a, lda, work( itau ), work( itemp ),
535  $ lwork-itemp+1, info )
536 *
537 * Copy R into WORK and bidiagonalize it:
538 * (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
539 *
540  iqrf = itemp
541  id = iqrf + n*n
542  ie = id + n
543  itauq = ie + n
544  itaup = itauq + n
545  itemp = itaup + n
546  CALL slacpy( 'U', n, n, a, lda, work( iqrf ), n )
547  CALL slaset( 'L', n-1, n-1, zero, zero, work( iqrf+1 ), n )
548  CALL sgebrd( n, n, work( iqrf ), n, work( id ), work( ie ),
549  $ work( itauq ), work( itaup ), work( itemp ),
550  $ lwork-itemp+1, info )
551 *
552 * Solve eigenvalue problem TGK*Z=Z*S.
553 * (Workspace: need 14*N + 2*N*(N+1))
554 *
555  itgkz = itemp
556  itemp = itgkz + n*(n*2+1)
557  CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ), work( ie ),
558  $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
559  $ n*2, work( itemp ), iwork, info)
560 *
561 * If needed, compute left singular vectors.
562 *
563  IF( wantu ) THEN
564  j = itgkz
565  DO i = 1, ns
566  CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
567  j = j + n*2
568  END DO
569  CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
570 *
571 * Call SORMBR to compute QB*UB.
572 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
573 *
574  CALL sormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
575  $ work( itauq ), u, ldu, work( itemp ),
576  $ lwork-itemp+1, info )
577 *
578 * Call SORMQR to compute Q*(QB*UB).
579 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
580 *
581  CALL sormqr( 'L', 'N', m, ns, n, a, lda,
582  $ work( itau ), u, ldu, work( itemp ),
583  $ lwork-itemp+1, info )
584  END IF
585 *
586 * If needed, compute right singular vectors.
587 *
588  IF( wantvt) THEN
589  j = itgkz + n
590  DO i = 1, ns
591  CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
592  j = j + n*2
593  END DO
594 *
595 * Call SORMBR to compute VB**T * PB**T
596 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
597 *
598  CALL sormbr( 'P', 'R', 'T', ns, n, n, work( iqrf ), n,
599  $ work( itaup ), vt, ldvt, work( itemp ),
600  $ lwork-itemp+1, info )
601  END IF
602  ELSE
603 *
604 * Path 2 (M at least N, but not much larger)
605 * Reduce A to bidiagonal form without QR decomposition
606 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
607 * U = QB * UB; V**T = VB**T * PB**T
608 *
609 * Bidiagonalize A
610 * (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
611 *
612  id = 1
613  ie = id + n
614  itauq = ie + n
615  itaup = itauq + n
616  itemp = itaup + n
617  CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
618  $ work( itauq ), work( itaup ), work( itemp ),
619  $ lwork-itemp+1, info )
620 *
621 * Solve eigenvalue problem TGK*Z=Z*S.
622 * (Workspace: need 14*N + 2*N*(N+1))
623 *
624  itgkz = itemp
625  itemp = itgkz + n*(n*2+1)
626  CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ), work( ie ),
627  $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
628  $ n*2, work( itemp ), iwork, info)
629 *
630 * If needed, compute left singular vectors.
631 *
632  IF( wantu ) THEN
633  j = itgkz
634  DO i = 1, ns
635  CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
636  j = j + n*2
637  END DO
638  CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
639 *
640 * Call SORMBR to compute QB*UB.
641 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
642 *
643  CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
644  $ work( itauq ), u, ldu, work( itemp ),
645  $ lwork-itemp+1, ierr )
646  END IF
647 *
648 * If needed, compute right singular vectors.
649 *
650  IF( wantvt) THEN
651  j = itgkz + n
652  DO i = 1, ns
653  CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
654  j = j + n*2
655  END DO
656 *
657 * Call SORMBR to compute VB**T * PB**T
658 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
659 *
660  CALL sormbr( 'P', 'R', 'T', ns, n, n, a, lda,
661  $ work( itaup ), vt, ldvt, work( itemp ),
662  $ lwork-itemp+1, ierr )
663  END IF
664  END IF
665  ELSE
666 *
667 * A has more columns than rows. If A has sufficiently more
668 * columns than rows, first reduce A using the LQ decomposition.
669 *
670  IF( n.GE.mnthr ) THEN
671 *
672 * Path 1t (N much larger than M):
673 * A = L * Q = ( QB * B * PB**T ) * Q
674 * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
675 * U = QB * UB ; V**T = VB**T * PB**T * Q
676 *
677 * Compute A=L*Q
678 * (Workspace: need 2*M, prefer M+M*NB)
679 *
680  itau = 1
681  itemp = itau + m
682  CALL sgelqf( m, n, a, lda, work( itau ), work( itemp ),
683  $ lwork-itemp+1, info )
684 
685 * Copy L into WORK and bidiagonalize it:
686 * (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
687 *
688  ilqf = itemp
689  id = ilqf + m*m
690  ie = id + m
691  itauq = ie + m
692  itaup = itauq + m
693  itemp = itaup + m
694  CALL slacpy( 'L', m, m, a, lda, work( ilqf ), m )
695  CALL slaset( 'U', m-1, m-1, zero, zero, work( ilqf+m ), m )
696  CALL sgebrd( m, m, work( ilqf ), m, work( id ), work( ie ),
697  $ work( itauq ), work( itaup ), work( itemp ),
698  $ lwork-itemp+1, info )
699 *
700 * Solve eigenvalue problem TGK*Z=Z*S.
701 * (Workspace: need 2*M*M+14*M)
702 *
703  itgkz = itemp
704  itemp = itgkz + m*(m*2+1)
705  CALL sbdsvdx( 'U', jobz, rngtgk, m, work( id ), work( ie ),
706  $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
707  $ m*2, work( itemp ), iwork, info)
708 *
709 * If needed, compute left singular vectors.
710 *
711  IF( wantu ) THEN
712  j = itgkz
713  DO i = 1, ns
714  CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
715  j = j + m*2
716  END DO
717 *
718 * Call SORMBR to compute QB*UB.
719 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
720 *
721  CALL sormbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
722  $ work( itauq ), u, ldu, work( itemp ),
723  $ lwork-itemp+1, info )
724  END IF
725 *
726 * If needed, compute right singular vectors.
727 *
728  IF( wantvt) THEN
729  j = itgkz + m
730  DO i = 1, ns
731  CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
732  j = j + m*2
733  END DO
734  CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
735 *
736 * Call SORMBR to compute (VB**T)*(PB**T)
737 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
738 *
739  CALL sormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
740  $ work( itaup ), vt, ldvt, work( itemp ),
741  $ lwork-itemp+1, info )
742 *
743 * Call SORMLQ to compute ((VB**T)*(PB**T))*Q.
744 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
745 *
746  CALL sormlq( 'R', 'N', ns, n, m, a, lda,
747  $ work( itau ), vt, ldvt, work( itemp ),
748  $ lwork-itemp+1, info )
749  END IF
750  ELSE
751 *
752 * Path 2t (N greater than M, but not much larger)
753 * Reduce to bidiagonal form without LQ decomposition
754 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
755 * U = QB * UB; V**T = VB**T * PB**T
756 *
757 * Bidiagonalize A
758 * (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
759 *
760  id = 1
761  ie = id + m
762  itauq = ie + m
763  itaup = itauq + m
764  itemp = itaup + m
765  CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
766  $ work( itauq ), work( itaup ), work( itemp ),
767  $ lwork-itemp+1, info )
768 *
769 * Solve eigenvalue problem TGK*Z=Z*S.
770 * (Workspace: need 2*M*M+14*M)
771 *
772  itgkz = itemp
773  itemp = itgkz + m*(m*2+1)
774  CALL sbdsvdx( 'L', jobz, rngtgk, m, work( id ), work( ie ),
775  $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
776  $ m*2, work( itemp ), iwork, info)
777 *
778 * If needed, compute left singular vectors.
779 *
780  IF( wantu ) THEN
781  j = itgkz
782  DO i = 1, ns
783  CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
784  j = j + m*2
785  END DO
786 *
787 * Call SORMBR to compute QB*UB.
788 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
789 *
790  CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
791  $ work( itauq ), u, ldu, work( itemp ),
792  $ lwork-itemp+1, info )
793  END IF
794 *
795 * If needed, compute right singular vectors.
796 *
797  IF( wantvt) THEN
798  j = itgkz + m
799  DO i = 1, ns
800  CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
801  j = j + m*2
802  END DO
803  CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
804 *
805 * Call SORMBR to compute VB**T * PB**T
806 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
807 *
808  CALL sormbr( 'P', 'R', 'T', ns, n, m, a, lda,
809  $ work( itaup ), vt, ldvt, work( itemp ),
810  $ lwork-itemp+1, info )
811  END IF
812  END IF
813  END IF
814 *
815 * Undo scaling if necessary
816 *
817  IF( iscl.EQ.1 ) THEN
818  IF( anrm.GT.bignum )
819  $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
820  $ s, minmn, info )
821  IF( anrm.LT.smlnum )
822  $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
823  $ s, minmn, info )
824  END IF
825 *
826 * Return optimal workspace in WORK(1)
827 *
828  work( 1 ) = REAL( maxwrk )
829 *
830  RETURN
831 *
832 * End of SGESVDX
833 *
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:170
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.
Definition: slascl.f:145
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:138
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
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...
Definition: slaset.f:112
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
Definition: sormlq.f:170
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
Definition: sormbr.f:198
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:116
subroutine sbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
SBDSVDX
Definition: sbdsvdx.f:228
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
Definition: sgebrd.f:207
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
Definition: sgelqf.f:137

Here is the call graph for this function:

Here is the caller graph for this function: