LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ sgesvdx()

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.

Definition at line 260 of file sgesvdx.f.

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