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

◆ cgesvdx()

subroutine cgesvdx ( character jobu,
character jobvt,
character range,
integer m,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
integer ns,
real, dimension( * ) s,
complex, dimension( ldu, * ) u,
integer ldu,
complex, dimension( ldvt, * ) vt,
integer ldvt,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer info )

CGESVDX computes the singular value decomposition (SVD) for GE matrices

Download CGESVDX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>  CGESVDX computes the singular value decomposition (SVD) of a complex
!>  M-by-N matrix A, optionally computing the left and/or right singular
!>  vectors. The SVD is written
!>
!>      A = U * SIGMA * transpose(V)
!>
!>  where SIGMA is an M-by-N matrix which is zero except for its
!>  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
!>  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
!>  are the singular values of A; they are real and non-negative, and
!>  are returned in descending order.  The first min(m,n) columns of
!>  U and V are the left and right singular vectors of A.
!>
!>  CGESVDX uses an eigenvalue problem for obtaining the SVD, which
!>  allows for the computation of a subset of singular values and
!>  vectors. See SBDSVDX for details.
!>
!>  Note that the routine returns V**T, not V.
!>
Parameters
[in]JOBU
!>          JOBU is CHARACTER*1
!>          Specifies options for computing all or part of the matrix U:
!>          = 'V':  the first min(m,n) columns of U (the left singular
!>                  vectors) or as specified by RANGE are returned in
!>                  the array U;
!>          = 'N':  no columns of U (no left singular vectors) are
!>                  computed.
!>
[in]JOBVT
!>          JOBVT is CHARACTER*1
!>           Specifies options for computing all or part of the matrix
!>           V**T:
!>           = 'V':  the first min(m,n) rows of V**T (the right singular
!>                   vectors) or as specified by RANGE are returned in
!>                   the array VT;
!>           = 'N':  no rows of V**T (no right singular vectors) are
!>                   computed.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all singular values will be found.
!>          = 'V': all singular values in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th singular values will be found.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the input matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the input matrix A.  N >= 0.
!>
[in,out]A
!>          A is COMPLEX array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, the contents of A are destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for singular values. VU > VL.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for singular values. VU > VL.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest singular value to be returned.
!>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest singular value to be returned.
!>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[out]NS
!>          NS is INTEGER
!>          The total number of singular values found,
!>          0 <= NS <= min(M,N).
!>          If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
!>
[out]S
!>          S is REAL array, dimension (min(M,N))
!>          The singular values of A, sorted so that S(i) >= S(i+1).
!>
[out]U
!>          U is COMPLEX array, dimension (LDU,UCOL)
!>          If JOBU = 'V', U contains columns of U (the left singular
!>          vectors, stored columnwise) as specified by RANGE; if
!>          JOBU = 'N', U is not referenced.
!>          Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
!>          the exact value of NS is not known in advance and an upper
!>          bound must be used.
!>
[in]LDU
!>          LDU is INTEGER
!>          The leading dimension of the array U.  LDU >= 1; if
!>          JOBU = 'V', LDU >= M.
!>
[out]VT
!>          VT is COMPLEX array, dimension (LDVT,N)
!>          If JOBVT = 'V', VT contains the rows of V**T (the right singular
!>          vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
!>          VT is not referenced.
!>          Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
!>          the exact value of NS is not known in advance and an upper
!>          bound must be used.
!>
[in]LDVT
!>          LDVT is INTEGER
!>          The leading dimension of the array VT.  LDVT >= 1; if
!>          JOBVT = 'V', LDVT >= NS (see above).
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
!>          comments inside the code):
!>             - PATH 1  (M much larger than N)
!>             - PATH 1t (N much larger than M)
!>          LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
!>          For good performance, LWORK should generally be larger.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (MAX(1,LRWORK))
!>          LRWORK >= MIN(M,N)*(MIN(M,N)*2+15*MIN(M,N)).
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (12*MIN(M,N))
!>          If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
!>          then IWORK contains the indices of the eigenvectors that failed
!>          to converge in SBDSVDX/SSTEVX.
!>
[out]INFO
!>     INFO is INTEGER
!>           = 0:  successful exit
!>           < 0:  if INFO = -i, the i-th argument had an illegal value
!>           > 0:  if INFO = i, then i eigenvectors failed to converge
!>                 in SBDSVDX/SSTEVX.
!>                 if INFO = N*2 + 1, an internal error occurred in
!>                 SBDSVDX
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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