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

◆ sgesvdx()

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

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

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

Purpose:
  SGESVDX computes the singular value decomposition (SVD) of a real
  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 orthogonal matrix, and
  V is an N-by-N orthogonal 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.

  SGESVDX 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 REAL 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 REAL 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 REAL 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 REAL 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]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 260 of file sgesvdx.f.

263*
264* -- LAPACK driver routine --
265* -- LAPACK is a software package provided by Univ. of Tennessee, --
266* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
267*
268* .. Scalar Arguments ..
269 CHARACTER JOBU, JOBVT, RANGE
270 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
271 REAL VL, VU
272* ..
273* .. Array Arguments ..
274 INTEGER IWORK( * )
275 REAL A( LDA, * ), S( * ), U( LDU, * ),
276 $ VT( LDVT, * ), WORK( * )
277* ..
278*
279* =====================================================================
280*
281* .. Parameters ..
282 REAL ZERO, ONE
283 parameter( zero = 0.0e0, one = 1.0e0 )
284* ..
285* .. Local Scalars ..
286 CHARACTER JOBZ, RNGTGK
287 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
288 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
289 $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
290 $ J, MAXWRK, MINMN, MINWRK, MNTHR
291 REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
292* ..
293* .. Local Arrays ..
294 REAL DUM( 1 )
295* ..
296* .. External Subroutines ..
297 EXTERNAL sbdsvdx, sgebrd, sgelqf, sgeqrf, slacpy,
299 $ scopy, xerbla
300* ..
301* .. External Functions ..
302 LOGICAL LSAME
303 INTEGER ILAENV
304 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
306* ..
307* .. Intrinsic Functions ..
308 INTRINSIC max, min, sqrt
309* ..
310* .. Executable Statements ..
311*
312* Test the input arguments.
313*
314 ns = 0
315 info = 0
316 abstol = 2*slamch('S')
317 lquery = ( lwork.EQ.-1 )
318 minmn = min( m, n )
319
320 wantu = lsame( jobu, 'V' )
321 wantvt = lsame( jobvt, 'V' )
322 IF( wantu .OR. wantvt ) THEN
323 jobz = 'V'
324 ELSE
325 jobz = 'N'
326 END IF
327 alls = lsame( range, 'A' )
328 vals = lsame( range, 'V' )
329 inds = lsame( range, 'I' )
330*
331 info = 0
332 IF( .NOT.lsame( jobu, 'V' ) .AND.
333 $ .NOT.lsame( jobu, 'N' ) ) THEN
334 info = -1
335 ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
336 $ .NOT.lsame( jobvt, 'N' ) ) THEN
337 info = -2
338 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
339 info = -3
340 ELSE IF( m.LT.0 ) THEN
341 info = -4
342 ELSE IF( n.LT.0 ) THEN
343 info = -5
344 ELSE IF( m.GT.lda ) THEN
345 info = -7
346 ELSE IF( minmn.GT.0 ) THEN
347 IF( vals ) THEN
348 IF( vl.LT.zero ) THEN
349 info = -8
350 ELSE IF( vu.LE.vl ) THEN
351 info = -9
352 END IF
353 ELSE IF( inds ) THEN
354 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
355 info = -10
356 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
357 info = -11
358 END IF
359 END IF
360 IF( info.EQ.0 ) THEN
361 IF( wantu .AND. ldu.LT.m ) THEN
362 info = -15
363 ELSE IF( wantvt ) THEN
364 IF( inds ) THEN
365 IF( ldvt.LT.iu-il+1 ) THEN
366 info = -17
367 END IF
368 ELSE IF( ldvt.LT.minmn ) THEN
369 info = -17
370 END IF
371 END IF
372 END IF
373 END IF
374*
375* Compute workspace
376* (Note: Comments in the code beginning "Workspace:" describe the
377* minimal amount of workspace needed at that point in the code,
378* as well as the preferred amount for good performance.
379* NB refers to the optimal block size for the immediately
380* following subroutine, as returned by ILAENV.)
381*
382 IF( info.EQ.0 ) THEN
383 minwrk = 1
384 maxwrk = 1
385 IF( minmn.GT.0 ) THEN
386 IF( m.GE.n ) THEN
387 mnthr = ilaenv( 6, 'SGESVD', jobu // jobvt, m, n, 0, 0 )
388 IF( m.GE.mnthr ) THEN
389*
390* Path 1 (M much larger than N)
391*
392 maxwrk = n +
393 $ n*ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
394 maxwrk = max( maxwrk, n*(n+5) + 2*n*
395 $ ilaenv( 1, 'SGEBRD', ' ', n, n, -1, -1 ) )
396 IF (wantu) THEN
397 maxwrk = max(maxwrk,n*(n*3+6)+n*
398 $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
399 END IF
400 IF (wantvt) THEN
401 maxwrk = max(maxwrk,n*(n*3+6)+n*
402 $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
403 END IF
404 minwrk = n*(n*3+20)
405 ELSE
406*
407* Path 2 (M at least N, but not much larger)
408*
409 maxwrk = 4*n + ( m+n )*
410 $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
411 IF (wantu) THEN
412 maxwrk = max(maxwrk,n*(n*2+5)+n*
413 $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
414 END IF
415 IF (wantvt) THEN
416 maxwrk = max(maxwrk,n*(n*2+5)+n*
417 $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
418 END IF
419 minwrk = max(n*(n*2+19),4*n+m)
420 END IF
421 ELSE
422 mnthr = ilaenv( 6, 'SGESVD', jobu // jobvt, m, n, 0, 0 )
423 IF( n.GE.mnthr ) THEN
424*
425* Path 1t (N much larger than M)
426*
427 maxwrk = m +
428 $ m*ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
429 maxwrk = max( maxwrk, m*(m+5) + 2*m*
430 $ ilaenv( 1, 'SGEBRD', ' ', m, m, -1, -1 ) )
431 IF (wantu) THEN
432 maxwrk = max(maxwrk,m*(m*3+6)+m*
433 $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
434 END IF
435 IF (wantvt) THEN
436 maxwrk = max(maxwrk,m*(m*3+6)+m*
437 $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
438 END IF
439 minwrk = m*(m*3+20)
440 ELSE
441*
442* Path 2t (N at least M, but not much larger)
443*
444 maxwrk = 4*m + ( m+n )*
445 $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
446 IF (wantu) THEN
447 maxwrk = max(maxwrk,m*(m*2+5)+m*
448 $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
449 END IF
450 IF (wantvt) THEN
451 maxwrk = max(maxwrk,m*(m*2+5)+m*
452 $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
453 END IF
454 minwrk = max(m*(m*2+19),4*m+n)
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( 'SGESVDX', -info )
468 RETURN
469 ELSE IF( lquery ) THEN
470 RETURN
471 END IF
472*
473* Quick return if possible
474*
475 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
476 RETURN
477 END IF
478*
479* Set singular values indices accord to RANGE.
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 = slange( 'M', m, n, a, lda, dum )
504 iscl = 0
505 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
506 iscl = 1
507 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
508 ELSE IF( anrm.GT.bignum ) THEN
509 iscl = 1
510 CALL slascl( '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 sgeqrf( 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+5*N, prefer N*N+4*N+2*N*NB)
536*
537 iqrf = itemp
538 id = iqrf + n*n
539 ie = id + n
540 itauq = ie + n
541 itaup = itauq + n
542 itemp = itaup + n
543 CALL slacpy( 'U', n, n, a, lda, work( iqrf ), n )
544 CALL slaset( 'L', n-1, n-1, zero, zero, work( iqrf+1 ), n )
545 CALL sgebrd( n, n, work( iqrf ), n, work( id ), work( ie ),
546 $ work( itauq ), work( itaup ), work( itemp ),
547 $ lwork-itemp+1, info )
548*
549* Solve eigenvalue problem TGK*Z=Z*S.
550* (Workspace: need 14*N + 2*N*(N+1))
551*
552 itgkz = itemp
553 itemp = itgkz + n*(n*2+1)
554 CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ), work( ie ),
555 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
556 $ n*2, work( itemp ), iwork, info)
557*
558* If needed, compute left singular vectors.
559*
560 IF( wantu ) THEN
561 j = itgkz
562 DO i = 1, ns
563 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
564 j = j + n*2
565 END DO
566 CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
567*
568* Call SORMBR to compute QB*UB.
569* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
570*
571 CALL sormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
572 $ work( itauq ), u, ldu, work( itemp ),
573 $ lwork-itemp+1, info )
574*
575* Call SORMQR to compute Q*(QB*UB).
576* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
577*
578 CALL sormqr( 'L', 'N', m, ns, n, a, lda,
579 $ work( itau ), u, ldu, work( itemp ),
580 $ lwork-itemp+1, info )
581 END IF
582*
583* If needed, compute right singular vectors.
584*
585 IF( wantvt) THEN
586 j = itgkz + n
587 DO i = 1, ns
588 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
589 j = j + n*2
590 END DO
591*
592* Call SORMBR to compute VB**T * PB**T
593* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
594*
595 CALL sormbr( 'P', 'R', 'T', 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 4*N+M, prefer 4*N+(M+N)*NB)
608*
609 id = 1
610 ie = id + n
611 itauq = ie + n
612 itaup = itauq + n
613 itemp = itaup + n
614 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
615 $ work( itauq ), work( itaup ), work( itemp ),
616 $ lwork-itemp+1, info )
617*
618* Solve eigenvalue problem TGK*Z=Z*S.
619* (Workspace: need 14*N + 2*N*(N+1))
620*
621 itgkz = itemp
622 itemp = itgkz + n*(n*2+1)
623 CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ), work( ie ),
624 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
625 $ n*2, work( itemp ), iwork, info)
626*
627* If needed, compute left singular vectors.
628*
629 IF( wantu ) THEN
630 j = itgkz
631 DO i = 1, ns
632 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
633 j = j + n*2
634 END DO
635 CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ), ldu )
636*
637* Call SORMBR to compute QB*UB.
638* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
639*
640 CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
641 $ work( itauq ), u, ldu, work( itemp ),
642 $ lwork-itemp+1, ierr )
643 END IF
644*
645* If needed, compute right singular vectors.
646*
647 IF( wantvt) THEN
648 j = itgkz + n
649 DO i = 1, ns
650 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
651 j = j + n*2
652 END DO
653*
654* Call SORMBR to compute VB**T * PB**T
655* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
656*
657 CALL sormbr( 'P', 'R', 'T', ns, n, n, a, lda,
658 $ work( itaup ), vt, ldvt, work( itemp ),
659 $ lwork-itemp+1, ierr )
660 END IF
661 END IF
662 ELSE
663*
664* A has more columns than rows. If A has sufficiently more
665* columns than rows, first reduce A using the LQ decomposition.
666*
667 IF( n.GE.mnthr ) THEN
668*
669* Path 1t (N much larger than M):
670* A = L * Q = ( QB * B * PB**T ) * Q
671* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
672* U = QB * UB ; V**T = VB**T * PB**T * Q
673*
674* Compute A=L*Q
675* (Workspace: need 2*M, prefer M+M*NB)
676*
677 itau = 1
678 itemp = itau + m
679 CALL sgelqf( m, n, a, lda, work( itau ), work( itemp ),
680 $ lwork-itemp+1, info )
681
682* Copy L into WORK and bidiagonalize it:
683* (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
684*
685 ilqf = itemp
686 id = ilqf + m*m
687 ie = id + m
688 itauq = ie + m
689 itaup = itauq + m
690 itemp = itaup + m
691 CALL slacpy( 'L', m, m, a, lda, work( ilqf ), m )
692 CALL slaset( 'U', m-1, m-1, zero, zero, work( ilqf+m ), m )
693 CALL sgebrd( m, m, work( ilqf ), m, work( id ), work( ie ),
694 $ work( itauq ), work( itaup ), work( itemp ),
695 $ lwork-itemp+1, info )
696*
697* Solve eigenvalue problem TGK*Z=Z*S.
698* (Workspace: need 2*M*M+14*M)
699*
700 itgkz = itemp
701 itemp = itgkz + m*(m*2+1)
702 CALL sbdsvdx( 'U', jobz, rngtgk, m, work( id ), work( ie ),
703 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
704 $ m*2, work( itemp ), iwork, info)
705*
706* If needed, compute left singular vectors.
707*
708 IF( wantu ) THEN
709 j = itgkz
710 DO i = 1, ns
711 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
712 j = j + m*2
713 END DO
714*
715* Call SORMBR to compute QB*UB.
716* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
717*
718 CALL sormbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
719 $ work( itauq ), u, ldu, work( itemp ),
720 $ lwork-itemp+1, info )
721 END IF
722*
723* If needed, compute right singular vectors.
724*
725 IF( wantvt) THEN
726 j = itgkz + m
727 DO i = 1, ns
728 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
729 j = j + m*2
730 END DO
731 CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
732*
733* Call SORMBR to compute (VB**T)*(PB**T)
734* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
735*
736 CALL sormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
737 $ work( itaup ), vt, ldvt, work( itemp ),
738 $ lwork-itemp+1, info )
739*
740* Call SORMLQ to compute ((VB**T)*(PB**T))*Q.
741* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
742*
743 CALL sormlq( 'R', 'N', ns, n, m, a, lda,
744 $ work( itau ), vt, ldvt, work( itemp ),
745 $ lwork-itemp+1, info )
746 END IF
747 ELSE
748*
749* Path 2t (N greater than M, but not much larger)
750* Reduce to bidiagonal form without LQ decomposition
751* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
752* U = QB * UB; V**T = VB**T * PB**T
753*
754* Bidiagonalize A
755* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
756*
757 id = 1
758 ie = id + m
759 itauq = ie + m
760 itaup = itauq + m
761 itemp = itaup + m
762 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
763 $ work( itauq ), work( itaup ), work( itemp ),
764 $ lwork-itemp+1, info )
765*
766* Solve eigenvalue problem TGK*Z=Z*S.
767* (Workspace: need 2*M*M+14*M)
768*
769 itgkz = itemp
770 itemp = itgkz + m*(m*2+1)
771 CALL sbdsvdx( 'L', jobz, rngtgk, m, work( id ), work( ie ),
772 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
773 $ m*2, work( itemp ), iwork, info)
774*
775* If needed, compute left singular vectors.
776*
777 IF( wantu ) THEN
778 j = itgkz
779 DO i = 1, ns
780 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
781 j = j + m*2
782 END DO
783*
784* Call SORMBR to compute QB*UB.
785* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
786*
787 CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
788 $ work( itauq ), u, ldu, work( itemp ),
789 $ lwork-itemp+1, info )
790 END IF
791*
792* If needed, compute right singular vectors.
793*
794 IF( wantvt) THEN
795 j = itgkz + m
796 DO i = 1, ns
797 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
798 j = j + m*2
799 END DO
800 CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ), ldvt)
801*
802* Call SORMBR to compute VB**T * PB**T
803* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
804*
805 CALL sormbr( 'P', 'R', 'T', ns, n, m, a, lda,
806 $ work( itaup ), vt, ldvt, work( itemp ),
807 $ lwork-itemp+1, info )
808 END IF
809 END IF
810 END IF
811*
812* Undo scaling if necessary
813*
814 IF( iscl.EQ.1 ) THEN
815 IF( anrm.GT.bignum )
816 $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
817 $ s, minmn, info )
818 IF( anrm.LT.smlnum )
819 $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
820 $ s, minmn, info )
821 END IF
822*
823* Return optimal workspace in WORK(1)
824*
825 work( 1 ) = sroundup_lwork( maxwrk )
826*
827 RETURN
828*
829* End of SGESVDX
830*
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 scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
SGEBRD
Definition sgebrd.f:205
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:143
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
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 slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMBR
Definition sormbr.f:196
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
Definition sormlq.f:168
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
Here is the call graph for this function:
Here is the caller graph for this function: