LAPACK  3.9.0
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.
Date
June 2016

Definition at line 272 of file cgesvdx.f.

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