LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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:152
Here is the call graph for this function:
Here is the caller graph for this function: