LAPACK 3.12.1
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 258 of file sgesvdx.f.

261*
262* -- LAPACK driver routine --
263* -- LAPACK is a software package provided by Univ. of Tennessee, --
264* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265*
266* .. Scalar Arguments ..
267 CHARACTER JOBU, JOBVT, RANGE
268 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
269 REAL VL, VU
270* ..
271* .. Array Arguments ..
272 INTEGER IWORK( * )
273 REAL A( LDA, * ), S( * ), U( LDU, * ),
274 $ VT( LDVT, * ), WORK( * )
275* ..
276*
277* =====================================================================
278*
279* .. Parameters ..
280 REAL ZERO, ONE
281 parameter( zero = 0.0e0, one = 1.0e0 )
282* ..
283* .. Local Scalars ..
284 CHARACTER JOBZ, RNGTGK
285 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
286 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
287 $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
288 $ J, MAXWRK, MINMN, MINWRK, MNTHR
289 REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
290* ..
291* .. Local Arrays ..
292 REAL DUM( 1 )
293* ..
294* .. External Subroutines ..
295 EXTERNAL sbdsvdx, sgebrd, sgelqf, sgeqrf,
296 $ slacpy,
297 $ slascl, slaset, sormbr, sormlq, sormqr,
298 $ scopy, xerbla
299* ..
300* .. External Functions ..
301 LOGICAL LSAME
302 INTEGER ILAENV
303 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
304 EXTERNAL lsame, ilaenv, slamch, slange,
305 $ 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,
388 $ 0 )
389 IF( m.GE.mnthr ) THEN
390*
391* Path 1 (M much larger than N)
392*
393 maxwrk = n +
394 $ n*ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
395 maxwrk = max( maxwrk, n*(n+5) + 2*n*
396 $ ilaenv( 1, 'SGEBRD', ' ', n, n, -1, -1 ) )
397 IF (wantu) THEN
398 maxwrk = max(maxwrk,n*(n*3+6)+n*
399 $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
400 END IF
401 IF (wantvt) THEN
402 maxwrk = max(maxwrk,n*(n*3+6)+n*
403 $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
404 END IF
405 minwrk = n*(n*3+20)
406 ELSE
407*
408* Path 2 (M at least N, but not much larger)
409*
410 maxwrk = 4*n + ( m+n )*
411 $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
412 IF (wantu) THEN
413 maxwrk = max(maxwrk,n*(n*2+5)+n*
414 $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
415 END IF
416 IF (wantvt) THEN
417 maxwrk = max(maxwrk,n*(n*2+5)+n*
418 $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
419 END IF
420 minwrk = max(n*(n*2+19),4*n+m)
421 END IF
422 ELSE
423 mnthr = ilaenv( 6, 'SGESVD', jobu // jobvt, m, n, 0,
424 $ 0 )
425 IF( n.GE.mnthr ) THEN
426*
427* Path 1t (N much larger than M)
428*
429 maxwrk = m +
430 $ m*ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
431 maxwrk = max( maxwrk, m*(m+5) + 2*m*
432 $ ilaenv( 1, 'SGEBRD', ' ', m, m, -1, -1 ) )
433 IF (wantu) THEN
434 maxwrk = max(maxwrk,m*(m*3+6)+m*
435 $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
436 END IF
437 IF (wantvt) THEN
438 maxwrk = max(maxwrk,m*(m*3+6)+m*
439 $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
440 END IF
441 minwrk = m*(m*3+20)
442 ELSE
443*
444* Path 2t (N at least M, but not much larger)
445*
446 maxwrk = 4*m + ( m+n )*
447 $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
448 IF (wantu) THEN
449 maxwrk = max(maxwrk,m*(m*2+5)+m*
450 $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
451 END IF
452 IF (wantvt) THEN
453 maxwrk = max(maxwrk,m*(m*2+5)+m*
454 $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
455 END IF
456 minwrk = max(m*(m*2+19),4*m+n)
457 END IF
458 END IF
459 END IF
460 maxwrk = max( maxwrk, minwrk )
461 work( 1 ) = sroundup_lwork( maxwrk )
462*
463 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
464 info = -19
465 END IF
466 END IF
467*
468 IF( info.NE.0 ) THEN
469 CALL xerbla( 'SGESVDX', -info )
470 RETURN
471 ELSE IF( lquery ) THEN
472 RETURN
473 END IF
474*
475* Quick return if possible
476*
477 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
478 RETURN
479 END IF
480*
481* Set singular values indices accord to RANGE.
482*
483 IF( alls ) THEN
484 rngtgk = 'I'
485 iltgk = 1
486 iutgk = min( m, n )
487 ELSE IF( inds ) THEN
488 rngtgk = 'I'
489 iltgk = il
490 iutgk = iu
491 ELSE
492 rngtgk = 'V'
493 iltgk = 0
494 iutgk = 0
495 END IF
496*
497* Get machine constants
498*
499 eps = slamch( 'P' )
500 smlnum = sqrt( slamch( 'S' ) ) / eps
501 bignum = one / smlnum
502*
503* Scale A if max element outside range [SMLNUM,BIGNUM]
504*
505 anrm = slange( 'M', m, n, a, lda, dum )
506 iscl = 0
507 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
508 iscl = 1
509 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
510 ELSE IF( anrm.GT.bignum ) THEN
511 iscl = 1
512 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
513 END IF
514*
515 IF( m.GE.n ) THEN
516*
517* A has at least as many rows as columns. If A has sufficiently
518* more rows than columns, first reduce A using the QR
519* decomposition.
520*
521 IF( m.GE.mnthr ) THEN
522*
523* Path 1 (M much larger than N):
524* A = Q * R = Q * ( QB * B * PB**T )
525* = Q * ( QB * ( UB * S * VB**T ) * PB**T )
526* U = Q * QB * UB; V**T = VB**T * PB**T
527*
528* Compute A=Q*R
529* (Workspace: need 2*N, prefer N+N*NB)
530*
531 itau = 1
532 itemp = itau + n
533 CALL sgeqrf( m, n, a, lda, work( itau ), work( itemp ),
534 $ lwork-itemp+1, info )
535*
536* Copy R into WORK and bidiagonalize it:
537* (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
538*
539 iqrf = itemp
540 id = iqrf + n*n
541 ie = id + n
542 itauq = ie + n
543 itaup = itauq + n
544 itemp = itaup + n
545 CALL slacpy( 'U', n, n, a, lda, work( iqrf ), n )
546 CALL slaset( 'L', n-1, n-1, zero, zero, work( iqrf+1 ),
547 $ n )
548 CALL sgebrd( n, n, work( iqrf ), n, work( id ),
549 $ work( ie ),
550 $ work( itauq ), work( itaup ), work( itemp ),
551 $ lwork-itemp+1, info )
552*
553* Solve eigenvalue problem TGK*Z=Z*S.
554* (Workspace: need 14*N + 2*N*(N+1))
555*
556 itgkz = itemp
557 itemp = itgkz + n*(n*2+1)
558 CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ),
559 $ work( ie ),
560 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
561 $ n*2, work( itemp ), iwork, info)
562*
563* If needed, compute left singular vectors.
564*
565 IF( wantu ) THEN
566 j = itgkz
567 DO i = 1, ns
568 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
569 j = j + n*2
570 END DO
571 CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ),
572 $ ldu )
573*
574* Call SORMBR to compute QB*UB.
575* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
576*
577 CALL sormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
578 $ work( itauq ), u, ldu, work( itemp ),
579 $ lwork-itemp+1, info )
580*
581* Call SORMQR to compute Q*(QB*UB).
582* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
583*
584 CALL sormqr( '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 j = itgkz + n
593 DO i = 1, ns
594 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
595 j = j + n*2
596 END DO
597*
598* Call SORMBR to compute VB**T * PB**T
599* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
600*
601 CALL sormbr( 'P', 'R', 'T', ns, n, n, work( iqrf ), n,
602 $ work( itaup ), vt, ldvt, work( itemp ),
603 $ lwork-itemp+1, info )
604 END IF
605 ELSE
606*
607* Path 2 (M at least N, but not much larger)
608* Reduce A to bidiagonal form without QR decomposition
609* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
610* U = QB * UB; V**T = VB**T * PB**T
611*
612* Bidiagonalize A
613* (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
614*
615 id = 1
616 ie = id + n
617 itauq = ie + n
618 itaup = itauq + n
619 itemp = itaup + n
620 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
621 $ work( itauq ), work( itaup ), work( itemp ),
622 $ lwork-itemp+1, info )
623*
624* Solve eigenvalue problem TGK*Z=Z*S.
625* (Workspace: need 14*N + 2*N*(N+1))
626*
627 itgkz = itemp
628 itemp = itgkz + n*(n*2+1)
629 CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ),
630 $ work( ie ),
631 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
632 $ n*2, work( itemp ), iwork, info)
633*
634* If needed, compute left singular vectors.
635*
636 IF( wantu ) THEN
637 j = itgkz
638 DO i = 1, ns
639 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
640 j = j + n*2
641 END DO
642 CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ),
643 $ ldu )
644*
645* Call SORMBR to compute QB*UB.
646* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
647*
648 CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
649 $ work( itauq ), u, ldu, work( itemp ),
650 $ lwork-itemp+1, ierr )
651 END IF
652*
653* If needed, compute right singular vectors.
654*
655 IF( wantvt) THEN
656 j = itgkz + n
657 DO i = 1, ns
658 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
659 j = j + n*2
660 END DO
661*
662* Call SORMBR to compute VB**T * PB**T
663* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
664*
665 CALL sormbr( 'P', 'R', 'T', ns, n, n, a, lda,
666 $ work( itaup ), vt, ldvt, work( itemp ),
667 $ lwork-itemp+1, ierr )
668 END IF
669 END IF
670 ELSE
671*
672* A has more columns than rows. If A has sufficiently more
673* columns than rows, first reduce A using the LQ decomposition.
674*
675 IF( n.GE.mnthr ) THEN
676*
677* Path 1t (N much larger than M):
678* A = L * Q = ( QB * B * PB**T ) * Q
679* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
680* U = QB * UB ; V**T = VB**T * PB**T * Q
681*
682* Compute A=L*Q
683* (Workspace: need 2*M, prefer M+M*NB)
684*
685 itau = 1
686 itemp = itau + m
687 CALL sgelqf( m, n, a, lda, work( itau ), work( itemp ),
688 $ lwork-itemp+1, info )
689
690* Copy L into WORK and bidiagonalize it:
691* (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
692*
693 ilqf = itemp
694 id = ilqf + m*m
695 ie = id + m
696 itauq = ie + m
697 itaup = itauq + m
698 itemp = itaup + m
699 CALL slacpy( 'L', m, m, a, lda, work( ilqf ), m )
700 CALL slaset( 'U', m-1, m-1, zero, zero, work( ilqf+m ),
701 $ m )
702 CALL sgebrd( m, m, work( ilqf ), m, work( id ),
703 $ work( ie ),
704 $ work( itauq ), work( itaup ), work( itemp ),
705 $ lwork-itemp+1, info )
706*
707* Solve eigenvalue problem TGK*Z=Z*S.
708* (Workspace: need 2*M*M+14*M)
709*
710 itgkz = itemp
711 itemp = itgkz + m*(m*2+1)
712 CALL sbdsvdx( 'U', jobz, rngtgk, m, work( id ),
713 $ work( ie ),
714 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
715 $ m*2, work( itemp ), iwork, info)
716*
717* If needed, compute left singular vectors.
718*
719 IF( wantu ) THEN
720 j = itgkz
721 DO i = 1, ns
722 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
723 j = j + m*2
724 END DO
725*
726* Call SORMBR to compute QB*UB.
727* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
728*
729 CALL sormbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
730 $ work( itauq ), u, ldu, work( itemp ),
731 $ lwork-itemp+1, info )
732 END IF
733*
734* If needed, compute right singular vectors.
735*
736 IF( wantvt) THEN
737 j = itgkz + m
738 DO i = 1, ns
739 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
740 j = j + m*2
741 END DO
742 CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ),
743 $ ldvt)
744*
745* Call SORMBR to compute (VB**T)*(PB**T)
746* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
747*
748 CALL sormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
749 $ work( itaup ), vt, ldvt, work( itemp ),
750 $ lwork-itemp+1, info )
751*
752* Call SORMLQ to compute ((VB**T)*(PB**T))*Q.
753* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
754*
755 CALL sormlq( 'R', 'N', ns, n, m, a, lda,
756 $ work( itau ), vt, ldvt, work( itemp ),
757 $ lwork-itemp+1, info )
758 END IF
759 ELSE
760*
761* Path 2t (N greater than M, but not much larger)
762* Reduce to bidiagonal form without LQ decomposition
763* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
764* U = QB * UB; V**T = VB**T * PB**T
765*
766* Bidiagonalize A
767* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
768*
769 id = 1
770 ie = id + m
771 itauq = ie + m
772 itaup = itauq + m
773 itemp = itaup + m
774 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
775 $ work( itauq ), work( itaup ), work( itemp ),
776 $ lwork-itemp+1, info )
777*
778* Solve eigenvalue problem TGK*Z=Z*S.
779* (Workspace: need 2*M*M+14*M)
780*
781 itgkz = itemp
782 itemp = itgkz + m*(m*2+1)
783 CALL sbdsvdx( 'L', jobz, rngtgk, m, work( id ),
784 $ work( ie ),
785 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
786 $ m*2, work( itemp ), iwork, info)
787*
788* If needed, compute left singular vectors.
789*
790 IF( wantu ) THEN
791 j = itgkz
792 DO i = 1, ns
793 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
794 j = j + m*2
795 END DO
796*
797* Call SORMBR to compute QB*UB.
798* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
799*
800 CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
801 $ work( itauq ), u, ldu, work( itemp ),
802 $ lwork-itemp+1, info )
803 END IF
804*
805* If needed, compute right singular vectors.
806*
807 IF( wantvt) THEN
808 j = itgkz + m
809 DO i = 1, ns
810 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
811 j = j + m*2
812 END DO
813 CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ),
814 $ ldvt)
815*
816* Call SORMBR to compute VB**T * PB**T
817* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
818*
819 CALL sormbr( 'P', 'R', 'T', ns, n, m, a, lda,
820 $ work( itaup ), vt, ldvt, work( itemp ),
821 $ lwork-itemp+1, info )
822 END IF
823 END IF
824 END IF
825*
826* Undo scaling if necessary
827*
828 IF( iscl.EQ.1 ) THEN
829 IF( anrm.GT.bignum )
830 $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
831 $ s, minmn, info )
832 IF( anrm.LT.smlnum )
833 $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
834 $ s, minmn, info )
835 END IF
836*
837* Return optimal workspace in WORK(1)
838*
839 work( 1 ) = sroundup_lwork( maxwrk )
840*
841 RETURN
842*
843* End of SGESVDX
844*
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
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgebrd
subroutine sgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
SGEBRD
Definition sgebrd.f:204
sgelqf
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:142
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:144
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
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:112
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
slaset
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:108
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
sormbr
subroutine sormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMBR
Definition sormbr.f:194
sormlq
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
Definition sormlq.f:166
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:166
Here is the call graph for this function:
Here is the caller graph for this function: