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

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

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

Purpose:
  ZGESVDX computes the singular value decomposition (SVD) of a complex
  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 unitary matrix, and
  V is an N-by-N unitary 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.

  ZGESVDX 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
          LRWORK >= MIN(M,N)*(MIN(M,N)*2+15*MIN(M,N)).
[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 272 of file zgesvdx.f.

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