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