 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dgesvdx()

 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```

Definition at line 260 of file dgesvdx.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  DOUBLE PRECISION VL, VU
272 * ..
273 * .. Array Arguments ..
274  INTEGER IWORK( * )
275  DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
276  \$ VT( LDVT, * ), WORK( * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  DOUBLE PRECISION ZERO, ONE
283  parameter( zero = 0.0d0, one = 1.0d0 )
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  DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
292 * ..
293 * .. Local Arrays ..
294  DOUBLE PRECISION DUM( 1 )
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL dbdsvdx, dgebrd, dgelqf, dgeqrf, dlacpy,
299  \$ dcopy, xerbla
300 * ..
301 * .. External Functions ..
302  LOGICAL LSAME
303  INTEGER ILAENV
304  DOUBLE PRECISION DLAMCH, DLANGE
305  EXTERNAL lsame, ilaenv, dlamch, dlange
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*dlamch('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, 'DGESVD', 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, 'DGEQRF', ' ', m, n, -1, -1 )
394  maxwrk = max( maxwrk, n*(n+5) + 2*n*
395  \$ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
396  IF (wantu) THEN
397  maxwrk = max(maxwrk,n*(n*3+6)+n*
398  \$ ilaenv( 1, 'DORMQR', ' ', n, n, -1, -1 ) )
399  END IF
400  IF (wantvt) THEN
401  maxwrk = max(maxwrk,n*(n*3+6)+n*
402  \$ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGEBRD', ' ', m, n, -1, -1 )
411  IF (wantu) THEN
412  maxwrk = max(maxwrk,n*(n*2+5)+n*
413  \$ ilaenv( 1, 'DORMQR', ' ', n, n, -1, -1 ) )
414  END IF
415  IF (wantvt) THEN
416  maxwrk = max(maxwrk,n*(n*2+5)+n*
417  \$ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGESVD', 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, 'DGELQF', ' ', m, n, -1, -1 )
429  maxwrk = max( maxwrk, m*(m+5) + 2*m*
430  \$ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
431  IF (wantu) THEN
432  maxwrk = max(maxwrk,m*(m*3+6)+m*
433  \$ ilaenv( 1, 'DORMQR', ' ', m, m, -1, -1 ) )
434  END IF
435  IF (wantvt) THEN
436  maxwrk = max(maxwrk,m*(m*3+6)+m*
437  \$ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGEBRD', ' ', m, n, -1, -1 )
446  IF (wantu) THEN
447  maxwrk = max(maxwrk,m*(m*2+5)+m*
448  \$ ilaenv( 1, 'DORMQR', ' ', m, m, -1, -1 ) )
449  END IF
450  IF (wantvt) THEN
451  maxwrk = max(maxwrk,m*(m*2+5)+m*
452  \$ ilaenv( 1, 'DORMLQ', ' ', 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 ) = dble( 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( 'DGESVDX', -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 = dlamch( 'P' )
498  smlnum = sqrt( dlamch( 'S' ) ) / eps
499  bignum = one / smlnum
500 *
501 * Scale A if max element outside range [SMLNUM,BIGNUM]
502 *
503  anrm = dlange( 'M', m, n, a, lda, dum )
504  iscl = 0
505  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
506  iscl = 1
507  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
508  ELSE IF( anrm.GT.bignum ) THEN
509  iscl = 1
510  CALL dlascl( '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 dgeqrf( 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 dlacpy( 'U', n, n, a, lda, work( iqrf ), n )
544  CALL dlaset( 'L', n-1, n-1, zero, zero, work( iqrf+1 ), n )
545  CALL dgebrd( 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 dbdsvdx( '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 dcopy( n, work( j ), 1, u( 1,i ), 1 )
564  j = j + n*2
565  END DO
566  CALL dlaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
567 *
568 * Call DORMBR to compute QB*UB.
569 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
570 *
571  CALL dormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
572  \$ work( itauq ), u, ldu, work( itemp ),
573  \$ lwork-itemp+1, info )
574 *
575 * Call DORMQR to compute Q*(QB*UB).
576 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
577 *
578  CALL dormqr( '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 dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
589  j = j + n*2
590  END DO
591 *
592 * Call DORMBR to compute VB**T * PB**T
593 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
594 *
595  CALL dormbr( '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 dgebrd( 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 dbdsvdx( '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 dcopy( n, work( j ), 1, u( 1,i ), 1 )
633  j = j + n*2
634  END DO
635  CALL dlaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
636 *
637 * Call DORMBR to compute QB*UB.
638 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
639 *
640  CALL dormbr( '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 dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
651  j = j + n*2
652  END DO
653 *
654 * Call DORMBR to compute VB**T * PB**T
655 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
656 *
657  CALL dormbr( '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 dgelqf( 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 dlacpy( 'L', m, m, a, lda, work( ilqf ), m )
692  CALL dlaset( 'U', m-1, m-1, zero, zero, work( ilqf+m ), m )
693  CALL dgebrd( 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 dbdsvdx( '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 dcopy( m, work( j ), 1, u( 1,i ), 1 )
712  j = j + m*2
713  END DO
714 *
715 * Call DORMBR to compute QB*UB.
716 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
717 *
718  CALL dormbr( '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 dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
729  j = j + m*2
730  END DO
731  CALL dlaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
732 *
733 * Call DORMBR to compute (VB**T)*(PB**T)
734 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
735 *
736  CALL dormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
737  \$ work( itaup ), vt, ldvt, work( itemp ),
738  \$ lwork-itemp+1, info )
739 *
740 * Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
741 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
742 *
743  CALL dormlq( '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 dgebrd( 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 dbdsvdx( '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 dcopy( m, work( j ), 1, u( 1,i ), 1 )
781  j = j + m*2
782  END DO
783 *
784 * Call DORMBR to compute QB*UB.
785 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
786 *
787  CALL dormbr( '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 dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
798  j = j + m*2
799  END DO
800  CALL dlaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
801 *
802 * Call DORMBR to compute VB**T * PB**T
803 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
804 *
805  CALL dormbr( '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 dlascl( 'G', 0, 0, bignum, anrm, minmn, 1,
817  \$ s, minmn, info )
818  IF( anrm.LT.smlnum )
819  \$ CALL dlascl( '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 ) = dble( maxwrk )
826 *
827  RETURN
828 *
829 * End of DGESVDX
830 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
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:110
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
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
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:114
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:146
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:143
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
Definition: dgebrd.f:205
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:167
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
Definition: dormbr.f:195
subroutine dbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
DBDSVDX
Definition: dbdsvdx.f:226
Here is the call graph for this function:
Here is the caller graph for this function: