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

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

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

Purpose:
  CGESVDX 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.
 
  CGESVDX 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 COMPLEX 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 COMPLEX 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 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 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 REAL 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 SBDSVDX/SSTEVX.
[out]INFO
     INFO is INTEGER
           = 0:  successful exit
           < 0:  if INFO = -i, the i-th argument had an illegal value
           > 0:  if INFO = i, then i eigenvectors failed to converge
                 in SBDSVDX/SSTEVX.
                 if INFO = N*2 + 1, an internal error occurred in
                 SBDSVDX
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 272 of file cgesvdx.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  REAL vl, vu
282 * ..
283 * .. Array Arguments ..
284  INTEGER iwork( * )
285  REAL s( * ), rwork( * )
286  COMPLEX a( lda, * ), u( ldu, * ), vt( ldvt, * ),
287  $ work( * )
288 * ..
289 *
290 * =====================================================================
291 *
292 * .. Parameters ..
293  COMPLEX czero, cone
294  parameter ( czero = ( 0.0e0, 0.0e0 ),
295  $ cone = ( 1.0e0, 0.0e0 ) )
296  REAL zero, one
297  parameter ( zero = 0.0e0, one = 1.0e0 )
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  REAL abstol, anrm, bignum, eps, smlnum
306 * ..
307 * .. Local Arrays ..
308  REAL dum( 1 )
309 * ..
310 * .. External Subroutines ..
311  EXTERNAL cgebrd, cgelqf, cgeqrf, clascl, claset,
312  $ slascl, xerbla
313 * ..
314 * .. External Functions ..
315  LOGICAL lsame
316  INTEGER ilaenv
317  REAL slamch, clange
318  EXTERNAL lsame, ilaenv, slamch, clange
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*slamch('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, 'CGESVD', 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,'CGEQRF',' ',m,n,-1,-1)
407  maxwrk = max(maxwrk,
408  $ n*n+2*n+2*n*ilaenv(1,'CGEBRD',' ',n,n,-1,-1))
409  IF (wantu .OR. wantvt) THEN
410  maxwrk = max(maxwrk,
411  $ n*n+2*n+n*ilaenv(1,'CUNMQR','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,'CGEBRD',' ',m,n,-1,-1)
419  IF (wantu .OR. wantvt) THEN
420  maxwrk = max(maxwrk,
421  $ 2*n+n*ilaenv(1,'CUNMQR','LN',n,n,n,-1))
422  END IF
423  END IF
424  ELSE
425  mnthr = ilaenv( 6, 'CGESVD', 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,'CGELQF',' ',m,n,-1,-1)
432  maxwrk = max(maxwrk,
433  $ m*m+2*m+2*m*ilaenv(1,'CGEBRD',' ',m,m,-1,-1))
434  IF (wantu .OR. wantvt) THEN
435  maxwrk = max(maxwrk,
436  $ m*m+2*m+m*ilaenv(1,'CUNMQR','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,'CGEBRD',' ',m,n,-1,-1)
445  IF (wantu .OR. wantvt) THEN
446  maxwrk = max(maxwrk,
447  $ 2*m+m*ilaenv(1,'CUNMQR','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 ) = cmplx( REAL( 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( 'CGESVDX', -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 = slamch( 'P' )
492  smlnum = sqrt( slamch( 'S' ) ) / eps
493  bignum = one / smlnum
494 *
495 * Scale A if max element outside range [SMLNUM,BIGNUM]
496 *
497  anrm = clange( 'M', m, n, a, lda, dum )
498  iscl = 0
499  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
500  iscl = 1
501  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
502  ELSE IF( anrm.GT.bignum ) THEN
503  iscl = 1
504  CALL clascl( '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 cgeqrf( 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 clacpy( 'U', n, n, a, lda, work( iqrf ), n )
539  CALL claset( 'L', n-1, n-1, czero, czero,
540  $ work( iqrf+1 ), n )
541  CALL cgebrd( 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 sbdsvdx( '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 ) = cmplx( rwork( k ), zero )
561  k = k + 1
562  END DO
563  k = k + n
564  END DO
565  CALL claset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
566 *
567 * Call CUNMBR to compute QB*UB.
568 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
569 *
570  CALL cunmbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
571  $ work( itauq ), u, ldu, work( itemp ),
572  $ lwork-itemp+1, info )
573 *
574 * Call CUNMQR to compute Q*(QB*UB).
575 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
576 *
577  CALL cunmqr( '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 ) = cmplx( rwork( k ), zero )
589  k = k + 1
590  END DO
591  k = k + n
592  END DO
593 *
594 * Call CUNMBR to compute VB**T * PB**T
595 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
596 *
597  CALL cunmbr( '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 cgebrd( 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 sbdsvdx( '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 ) = cmplx( rwork( k ), zero )
637  k = k + 1
638  END DO
639  k = k + n
640  END DO
641  CALL claset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
642 *
643 * Call CUNMBR to compute QB*UB.
644 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
645 *
646  CALL cunmbr( '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 ) = cmplx( rwork( k ), zero )
658  k = k + 1
659  END DO
660  k = k + n
661  END DO
662 *
663 * Call CUNMBR to compute VB**T * PB**T
664 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
665 *
666  CALL cunmbr( '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 cgelqf( 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 clacpy( 'L', m, m, a, lda, work( ilqf ), m )
702  CALL claset( 'U', m-1, m-1, czero, czero,
703  $ work( ilqf+m ), m )
704  CALL cgebrd( 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 sbdsvdx( '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 ) = cmplx( rwork( k ), zero )
724  k = k + 1
725  END DO
726  k = k + m
727  END DO
728 *
729 * Call CUNMBR to compute QB*UB.
730 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
731 *
732  CALL cunmbr( '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 ) = cmplx( rwork( k ), zero )
744  k = k + 1
745  END DO
746  k = k + m
747  END DO
748  CALL claset( 'A', ns, n-m, czero, czero,
749  $ vt( 1,m+1 ), ldvt )
750 *
751 * Call CUNMBR to compute (VB**T)*(PB**T)
752 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
753 *
754  CALL cunmbr( 'P', 'R', 'C', ns, m, m, work( ilqf ), m,
755  $ work( itaup ), vt, ldvt, work( itemp ),
756  $ lwork-itemp+1, info )
757 *
758 * Call CUNMLQ to compute ((VB**T)*(PB**T))*Q.
759 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
760 *
761  CALL cunmlq( '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 cgebrd( 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 sbdsvdx( '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 ) = cmplx( rwork( k ), zero )
801  k = k + 1
802  END DO
803  k = k + m
804  END DO
805 *
806 * Call CUNMBR to compute QB*UB.
807 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
808 *
809  CALL cunmbr( '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 ) = cmplx( rwork( k ), zero )
821  k = k + 1
822  END DO
823  k = k + m
824  END DO
825  CALL claset( 'A', ns, n-m, czero, czero,
826  $ vt( 1,m+1 ), ldvt )
827 *
828 * Call CUNMBR to compute VB**T * PB**T
829 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
830 *
831  CALL cunmbr( '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 slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
843  $ s, minmn, info )
844  IF( anrm.LT.smlnum )
845  $ CALL slascl( '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 ) = cmplx( REAL( MAXWRK ), zero )
852 *
853  RETURN
854 *
855 * End of CGESVDX
856 *
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
Definition: cunmbr.f:199
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:137
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine sbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
SBDSVDX
Definition: sbdsvdx.f:228
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:138
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
Definition: cgebrd.f:208
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
Definition: cunmlq.f:170
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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: