LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zgesvdx()

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.

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