LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ cgesvdx()

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.

Definition at line 267 of file cgesvdx.f.

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