LAPACK 3.12.1
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 265 of file zgesvdx.f.

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