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

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

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

Purpose:
  DGESVDX 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.
 
  DGESVDX uses an eigenvalue problem for obtaining the SVD, which 
  allows for the computation of a subset of singular values and 
  vectors. See DBDSVDX 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION
          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 DOUBLE PRECISION array, dimension (min(M,N))
          The singular values of A, sorted so that S(i) >= S(i+1).
[out]U
          U is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DBDSVDX/DSTEVX.
[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 DBDSVDX/DSTEVX.
                 if INFO = N*2 + 1, an internal error occurred in
                 DBDSVDX
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 265 of file dgesvdx.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  DOUBLE PRECISION vl, vu
275 * ..
276 * .. Array Arguments ..
277  INTEGER iwork( * )
278  DOUBLE PRECISION a( lda, * ), s( * ), u( ldu, * ),
279  $ vt( ldvt, * ), work( * )
280 * ..
281 *
282 * =====================================================================
283 *
284 * .. Parameters ..
285  DOUBLE PRECISION zero, one
286  parameter ( zero = 0.0d0, one = 1.0d0 )
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  DOUBLE PRECISION abstol, anrm, bignum, eps, smlnum
295 * ..
296 * .. Local Arrays ..
297  DOUBLE PRECISION dum( 1 )
298 * ..
299 * .. External Subroutines ..
300  EXTERNAL dbdsvdx, dgebrd, dgelqf, dgeqrf, dlacpy,
302  $ dscal, xerbla
303 * ..
304 * .. External Functions ..
305  LOGICAL lsame
306  INTEGER ilaenv
307  DOUBLE PRECISION dlamch, dlange, dnrm2
308  EXTERNAL lsame, ilaenv, dlamch, dlange, dnrm2
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*dlamch('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, 'DGESVD', 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, 'DGEQRF', ' ', m, n, -1, -1 )
397  maxwrk = max( maxwrk, n*(n+5) + 2*n*
398  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
399  IF (wantu) THEN
400  maxwrk = max(maxwrk,n*(n*3+6)+n*
401  $ ilaenv( 1, 'DORMQR', ' ', n, n, -1, -1 ) )
402  END IF
403  IF (wantvt) THEN
404  maxwrk = max(maxwrk,n*(n*3+6)+n*
405  $ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGEBRD', ' ', m, n, -1, -1 )
414  IF (wantu) THEN
415  maxwrk = max(maxwrk,n*(n*2+5)+n*
416  $ ilaenv( 1, 'DORMQR', ' ', n, n, -1, -1 ) )
417  END IF
418  IF (wantvt) THEN
419  maxwrk = max(maxwrk,n*(n*2+5)+n*
420  $ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGESVD', 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, 'DGELQF', ' ', m, n, -1, -1 )
432  maxwrk = max( maxwrk, m*(m+5) + 2*m*
433  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
434  IF (wantu) THEN
435  maxwrk = max(maxwrk,m*(m*3+6)+m*
436  $ ilaenv( 1, 'DORMQR', ' ', m, m, -1, -1 ) )
437  END IF
438  IF (wantvt) THEN
439  maxwrk = max(maxwrk,m*(m*3+6)+m*
440  $ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGEBRD', ' ', m, n, -1, -1 )
449  IF (wantu) THEN
450  maxwrk = max(maxwrk,m*(m*2+5)+m*
451  $ ilaenv( 1, 'DORMQR', ' ', m, m, -1, -1 ) )
452  END IF
453  IF (wantvt) THEN
454  maxwrk = max(maxwrk,m*(m*2+5)+m*
455  $ ilaenv( 1, 'DORMLQ', ' ', 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 ) = dble( 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( 'DGESVDX', -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 = dlamch( 'P' )
501  smlnum = sqrt( dlamch( 'S' ) ) / eps
502  bignum = one / smlnum
503 *
504 * Scale A if max element outside range [SMLNUM,BIGNUM]
505 *
506  anrm = dlange( 'M', m, n, a, lda, dum )
507  iscl = 0
508  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
509  iscl = 1
510  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
511  ELSE IF( anrm.GT.bignum ) THEN
512  iscl = 1
513  CALL dlascl( '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 dgeqrf( 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 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 )
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 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)
560 *
561 * If needed, compute left singular vectors.
562 *
563  IF( wantu ) THEN
564  j = itgkz
565  DO i = 1, ns
566  CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
567  j = j + n*2
568  END DO
569  CALL dlaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
570 *
571 * Call DORMBR to compute QB*UB.
572 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
573 *
574  CALL dormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
575  $ work( itauq ), u, ldu, work( itemp ),
576  $ lwork-itemp+1, info )
577 *
578 * Call DORMQR to compute Q*(QB*UB).
579 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
580 *
581  CALL dormqr( '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 dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
592  j = j + n*2
593  END DO
594 *
595 * Call DORMBR to compute VB**T * PB**T
596 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
597 *
598  CALL dormbr( '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 dgebrd( 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 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)
629 *
630 * If needed, compute left singular vectors.
631 *
632  IF( wantu ) THEN
633  j = itgkz
634  DO i = 1, ns
635  CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
636  j = j + n*2
637  END DO
638  CALL dlaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
639 *
640 * Call DORMBR to compute QB*UB.
641 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
642 *
643  CALL dormbr( '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 dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
654  j = j + n*2
655  END DO
656 *
657 * Call DORMBR to compute VB**T * PB**T
658 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
659 *
660  CALL dormbr( '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 dgelqf( 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 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 )
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 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)
708 *
709 * If needed, compute left singular vectors.
710 *
711  IF( wantu ) THEN
712  j = itgkz
713  DO i = 1, ns
714  CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
715  j = j + m*2
716  END DO
717 *
718 * Call DORMBR to compute QB*UB.
719 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
720 *
721  CALL dormbr( '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 dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
732  j = j + m*2
733  END DO
734  CALL dlaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
735 *
736 * Call DORMBR to compute (VB**T)*(PB**T)
737 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
738 *
739  CALL dormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
740  $ work( itaup ), vt, ldvt, work( itemp ),
741  $ lwork-itemp+1, info )
742 *
743 * Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
744 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
745 *
746  CALL dormlq( '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 dgebrd( 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 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)
777 *
778 * If needed, compute left singular vectors.
779 *
780  IF( wantu ) THEN
781  j = itgkz
782  DO i = 1, ns
783  CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
784  j = j + m*2
785  END DO
786 *
787 * Call DORMBR to compute QB*UB.
788 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
789 *
790  CALL dormbr( '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 dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
801  j = j + m*2
802  END DO
803  CALL dlaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
804 *
805 * Call DORMBR to compute VB**T * PB**T
806 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
807 *
808  CALL dormbr( '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 dlascl( 'G', 0, 0, bignum, anrm, minmn, 1,
820  $ s, minmn, info )
821  IF( anrm.LT.smlnum )
822  $ CALL dlascl( '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 ) = dble( maxwrk )
829 *
830  RETURN
831 *
832 * End of DGESVDX
833 *
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
Definition: dgebrd.f:207
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...
Definition: dlaset.f:112
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:169
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
Definition: dormbr.f:197
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
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.
Definition: dlascl.f:145
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:169
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:138
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine dbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
DBDSVDX
Definition: dbdsvdx.f:228
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: