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

◆ dgesvdx()

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

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

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

Purpose:
!>
!>  DGESVDX 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.
!>
!>  DGESVDX uses an eigenvalue problem for obtaining the SVD, which
!>  allows for the computation of a subset of singular values and
!>  vectors. See DBDSVDX for details.
!>
!>  Note that the routine returns V**T, not V.
!>
Parameters
[in]JOBU
!>          JOBU is CHARACTER*1
!>          Specifies options for computing all or part of the matrix U:
!>          = 'V':  the first min(m,n) columns of U (the left singular
!>                  vectors) or as specified by RANGE are returned in
!>                  the array U;
!>          = 'N':  no columns of U (no left singular vectors) are
!>                  computed.
!>
[in]JOBVT
!>          JOBVT is CHARACTER*1
!>           Specifies options for computing all or part of the matrix
!>           V**T:
!>           = 'V':  the first min(m,n) rows of V**T (the right singular
!>                   vectors) or as specified by RANGE are returned in
!>                   the array VT;
!>           = 'N':  no rows of V**T (no right singular vectors) are
!>                   computed.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all singular values will be found.
!>          = 'V': all singular values in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th singular values will be found.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the input matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the input matrix A.  N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, the contents of A are destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in]VL
!>          VL is DOUBLE PRECISION
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for singular values. VU > VL.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is DOUBLE PRECISION
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for singular values. VU > VL.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest singular value to be returned.
!>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest singular value to be returned.
!>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[out]NS
!>          NS is INTEGER
!>          The total number of singular values found,
!>          0 <= NS <= min(M,N).
!>          If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
!>
[out]S
!>          S is DOUBLE PRECISION array, dimension (min(M,N))
!>          The singular values of A, sorted so that S(i) >= S(i+1).
!>
[out]U
!>          U is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DBDSVDX/DSTEVX.
!>
[out]INFO
!>     INFO is INTEGER
!>           = 0:  successful exit
!>           < 0:  if INFO = -i, the i-th argument had an illegal value
!>           > 0:  if INFO = i, then i eigenvectors failed to converge
!>                 in DBDSVDX/DSTEVX.
!>                 if INFO = N*2 + 1, an internal error occurred in
!>                 DBDSVDX
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 258 of file dgesvdx.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 DOUBLE PRECISION VL, VU
270* ..
271* .. Array Arguments ..
272 INTEGER IWORK( * )
273 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
274 $ VT( LDVT, * ), WORK( * )
275* ..
276*
277* =====================================================================
278*
279* .. Parameters ..
280 DOUBLE PRECISION ZERO, ONE
281 parameter( zero = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
290* ..
291* .. Local Arrays ..
292 DOUBLE PRECISION DUM( 1 )
293* ..
294* .. External Subroutines ..
295 EXTERNAL dbdsvdx, dgebrd, dgelqf, dgeqrf,
296 $ dlacpy,
297 $ dlascl, dlaset, dormbr, dormlq, dormqr,
298 $ dcopy, xerbla
299* ..
300* .. External Functions ..
301 LOGICAL LSAME
302 INTEGER ILAENV
303 DOUBLE PRECISION DLAMCH, DLANGE
304 EXTERNAL lsame, ilaenv, dlamch, dlange
305* ..
306* .. Intrinsic Functions ..
307 INTRINSIC max, min, sqrt
308* ..
309* .. Executable Statements ..
310*
311* Test the input arguments.
312*
313 ns = 0
314 info = 0
315 abstol = 2*dlamch('S')
316 lquery = ( lwork.EQ.-1 )
317 minmn = min( m, n )
318
319 wantu = lsame( jobu, 'V' )
320 wantvt = lsame( jobvt, 'V' )
321 IF( wantu .OR. wantvt ) THEN
322 jobz = 'V'
323 ELSE
324 jobz = 'N'
325 END IF
326 alls = lsame( range, 'A' )
327 vals = lsame( range, 'V' )
328 inds = lsame( range, 'I' )
329*
330 info = 0
331 IF( .NOT.lsame( jobu, 'V' ) .AND.
332 $ .NOT.lsame( jobu, 'N' ) ) THEN
333 info = -1
334 ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
335 $ .NOT.lsame( jobvt, 'N' ) ) THEN
336 info = -2
337 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
338 info = -3
339 ELSE IF( m.LT.0 ) THEN
340 info = -4
341 ELSE IF( n.LT.0 ) THEN
342 info = -5
343 ELSE IF( m.GT.lda ) THEN
344 info = -7
345 ELSE IF( minmn.GT.0 ) THEN
346 IF( vals ) THEN
347 IF( vl.LT.zero ) THEN
348 info = -8
349 ELSE IF( vu.LE.vl ) THEN
350 info = -9
351 END IF
352 ELSE IF( inds ) THEN
353 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
354 info = -10
355 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
356 info = -11
357 END IF
358 END IF
359 IF( info.EQ.0 ) THEN
360 IF( wantu .AND. ldu.LT.m ) THEN
361 info = -15
362 ELSE IF( wantvt ) THEN
363 IF( inds ) THEN
364 IF( ldvt.LT.iu-il+1 ) THEN
365 info = -17
366 END IF
367 ELSE IF( ldvt.LT.minmn ) THEN
368 info = -17
369 END IF
370 END IF
371 END IF
372 END IF
373*
374* Compute workspace
375* (Note: Comments in the code beginning "Workspace:" describe the
376* minimal amount of workspace needed at that point in the code,
377* as well as the preferred amount for good performance.
378* NB refers to the optimal block size for the immediately
379* following subroutine, as returned by ILAENV.)
380*
381 IF( info.EQ.0 ) THEN
382 minwrk = 1
383 maxwrk = 1
384 IF( minmn.GT.0 ) THEN
385 IF( m.GE.n ) THEN
386 mnthr = ilaenv( 6, 'DGESVD', jobu // jobvt, m, n, 0,
387 $ 0 )
388 IF( m.GE.mnthr ) THEN
389*
390* Path 1 (M much larger than N)
391*
392 maxwrk = n +
393 $ n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
394 maxwrk = max( maxwrk, n*(n+5) + 2*n*
395 $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
396 IF (wantu) THEN
397 maxwrk = max(maxwrk,n*(n*3+6)+n*
398 $ ilaenv( 1, 'DORMQR', ' ', n, n, -1, -1 ) )
399 END IF
400 IF (wantvt) THEN
401 maxwrk = max(maxwrk,n*(n*3+6)+n*
402 $ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGEBRD', ' ', m, n, -1, -1 )
411 IF (wantu) THEN
412 maxwrk = max(maxwrk,n*(n*2+5)+n*
413 $ ilaenv( 1, 'DORMQR', ' ', n, n, -1, -1 ) )
414 END IF
415 IF (wantvt) THEN
416 maxwrk = max(maxwrk,n*(n*2+5)+n*
417 $ ilaenv( 1, 'DORMLQ', ' ', 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, 'DGESVD', jobu // jobvt, m, n, 0,
423 $ 0 )
424 IF( n.GE.mnthr ) THEN
425*
426* Path 1t (N much larger than M)
427*
428 maxwrk = m +
429 $ m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
430 maxwrk = max( maxwrk, m*(m+5) + 2*m*
431 $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
432 IF (wantu) THEN
433 maxwrk = max(maxwrk,m*(m*3+6)+m*
434 $ ilaenv( 1, 'DORMQR', ' ', m, m, -1, -1 ) )
435 END IF
436 IF (wantvt) THEN
437 maxwrk = max(maxwrk,m*(m*3+6)+m*
438 $ ilaenv( 1, 'DORMLQ', ' ', m, m, -1, -1 ) )
439 END IF
440 minwrk = m*(m*3+20)
441 ELSE
442*
443* Path 2t (N at least M, but not much larger)
444*
445 maxwrk = 4*m + ( m+n )*
446 $ ilaenv( 1, 'DGEBRD', ' ', m, n, -1, -1 )
447 IF (wantu) THEN
448 maxwrk = max(maxwrk,m*(m*2+5)+m*
449 $ ilaenv( 1, 'DORMQR', ' ', m, m, -1, -1 ) )
450 END IF
451 IF (wantvt) THEN
452 maxwrk = max(maxwrk,m*(m*2+5)+m*
453 $ ilaenv( 1, 'DORMLQ', ' ', m, m, -1, -1 ) )
454 END IF
455 minwrk = max(m*(m*2+19),4*m+n)
456 END IF
457 END IF
458 END IF
459 maxwrk = max( maxwrk, minwrk )
460 work( 1 ) = dble( maxwrk )
461*
462 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
463 info = -19
464 END IF
465 END IF
466*
467 IF( info.NE.0 ) THEN
468 CALL xerbla( 'DGESVDX', -info )
469 RETURN
470 ELSE IF( lquery ) THEN
471 RETURN
472 END IF
473*
474* Quick return if possible
475*
476 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
477 RETURN
478 END IF
479*
480* Set singular values indices accord to RANGE.
481*
482 IF( alls ) THEN
483 rngtgk = 'I'
484 iltgk = 1
485 iutgk = min( m, n )
486 ELSE IF( inds ) THEN
487 rngtgk = 'I'
488 iltgk = il
489 iutgk = iu
490 ELSE
491 rngtgk = 'V'
492 iltgk = 0
493 iutgk = 0
494 END IF
495*
496* Get machine constants
497*
498 eps = dlamch( 'P' )
499 smlnum = sqrt( dlamch( 'S' ) ) / eps
500 bignum = one / smlnum
501*
502* Scale A if max element outside range [SMLNUM,BIGNUM]
503*
504 anrm = dlange( 'M', m, n, a, lda, dum )
505 iscl = 0
506 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
507 iscl = 1
508 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
509 ELSE IF( anrm.GT.bignum ) THEN
510 iscl = 1
511 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
512 END IF
513*
514 IF( m.GE.n ) THEN
515*
516* A has at least as many rows as columns. If A has sufficiently
517* more rows than columns, first reduce A using the QR
518* decomposition.
519*
520 IF( m.GE.mnthr ) THEN
521*
522* Path 1 (M much larger than N):
523* A = Q * R = Q * ( QB * B * PB**T )
524* = Q * ( QB * ( UB * S * VB**T ) * PB**T )
525* U = Q * QB * UB; V**T = VB**T * PB**T
526*
527* Compute A=Q*R
528* (Workspace: need 2*N, prefer N+N*NB)
529*
530 itau = 1
531 itemp = itau + n
532 CALL dgeqrf( m, n, a, lda, work( itau ), work( itemp ),
533 $ lwork-itemp+1, info )
534*
535* Copy R into WORK and bidiagonalize it:
536* (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
537*
538 iqrf = itemp
539 id = iqrf + n*n
540 ie = id + n
541 itauq = ie + n
542 itaup = itauq + n
543 itemp = itaup + n
544 CALL dlacpy( 'U', n, n, a, lda, work( iqrf ), n )
545 CALL dlaset( 'L', n-1, n-1, zero, zero, work( iqrf+1 ),
546 $ n )
547 CALL dgebrd( n, n, work( iqrf ), n, work( id ),
548 $ work( ie ),
549 $ work( itauq ), work( itaup ), work( itemp ),
550 $ lwork-itemp+1, info )
551*
552* Solve eigenvalue problem TGK*Z=Z*S.
553* (Workspace: need 14*N + 2*N*(N+1))
554*
555 itgkz = itemp
556 itemp = itgkz + n*(n*2+1)
557 CALL dbdsvdx( 'U', jobz, rngtgk, n, work( id ),
558 $ work( ie ),
559 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
560 $ n*2, work( itemp ), iwork, info)
561*
562* If needed, compute left singular vectors.
563*
564 IF( wantu ) THEN
565 j = itgkz
566 DO i = 1, ns
567 CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
568 j = j + n*2
569 END DO
570 CALL dlaset( 'A', m-n, ns, zero, zero, u( n+1,1 ),
571 $ ldu )
572*
573* Call DORMBR to compute QB*UB.
574* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
575*
576 CALL dormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
577 $ work( itauq ), u, ldu, work( itemp ),
578 $ lwork-itemp+1, info )
579*
580* Call DORMQR to compute Q*(QB*UB).
581* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
582*
583 CALL dormqr( 'L', 'N', m, ns, n, a, lda,
584 $ work( itau ), u, ldu, work( itemp ),
585 $ lwork-itemp+1, info )
586 END IF
587*
588* If needed, compute right singular vectors.
589*
590 IF( wantvt) THEN
591 j = itgkz + n
592 DO i = 1, ns
593 CALL dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
594 j = j + n*2
595 END DO
596*
597* Call DORMBR to compute VB**T * PB**T
598* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
599*
600 CALL dormbr( 'P', 'R', 'T', ns, n, n, work( iqrf ), n,
601 $ work( itaup ), vt, ldvt, work( itemp ),
602 $ lwork-itemp+1, info )
603 END IF
604 ELSE
605*
606* Path 2 (M at least N, but not much larger)
607* Reduce A to bidiagonal form without QR decomposition
608* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
609* U = QB * UB; V**T = VB**T * PB**T
610*
611* Bidiagonalize A
612* (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
613*
614 id = 1
615 ie = id + n
616 itauq = ie + n
617 itaup = itauq + n
618 itemp = itaup + n
619 CALL dgebrd( m, n, a, lda, work( id ), work( ie ),
620 $ work( itauq ), work( itaup ), work( itemp ),
621 $ lwork-itemp+1, info )
622*
623* Solve eigenvalue problem TGK*Z=Z*S.
624* (Workspace: need 14*N + 2*N*(N+1))
625*
626 itgkz = itemp
627 itemp = itgkz + n*(n*2+1)
628 CALL dbdsvdx( 'U', jobz, rngtgk, n, work( id ),
629 $ work( ie ),
630 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
631 $ n*2, work( itemp ), iwork, info)
632*
633* If needed, compute left singular vectors.
634*
635 IF( wantu ) THEN
636 j = itgkz
637 DO i = 1, ns
638 CALL dcopy( n, work( j ), 1, u( 1,i ), 1 )
639 j = j + n*2
640 END DO
641 CALL dlaset( 'A', m-n, ns, zero, zero, u( n+1,1 ),
642 $ ldu )
643*
644* Call DORMBR to compute QB*UB.
645* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
646*
647 CALL dormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
648 $ work( itauq ), u, ldu, work( itemp ),
649 $ lwork-itemp+1, ierr )
650 END IF
651*
652* If needed, compute right singular vectors.
653*
654 IF( wantvt) THEN
655 j = itgkz + n
656 DO i = 1, ns
657 CALL dcopy( n, work( j ), 1, vt( i,1 ), ldvt )
658 j = j + n*2
659 END DO
660*
661* Call DORMBR to compute VB**T * PB**T
662* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
663*
664 CALL dormbr( 'P', 'R', 'T', 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 dgelqf( 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+5*N, prefer M*M+4*M+2*M*NB)
691*
692 ilqf = itemp
693 id = ilqf + m*m
694 ie = id + m
695 itauq = ie + m
696 itaup = itauq + m
697 itemp = itaup + m
698 CALL dlacpy( 'L', m, m, a, lda, work( ilqf ), m )
699 CALL dlaset( 'U', m-1, m-1, zero, zero, work( ilqf+m ),
700 $ m )
701 CALL dgebrd( m, m, work( ilqf ), m, work( id ),
702 $ work( ie ),
703 $ work( itauq ), work( itaup ), work( itemp ),
704 $ lwork-itemp+1, info )
705*
706* Solve eigenvalue problem TGK*Z=Z*S.
707* (Workspace: need 2*M*M+14*M)
708*
709 itgkz = itemp
710 itemp = itgkz + m*(m*2+1)
711 CALL dbdsvdx( 'U', jobz, rngtgk, m, work( id ),
712 $ work( ie ),
713 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
714 $ m*2, work( itemp ), iwork, info)
715*
716* If needed, compute left singular vectors.
717*
718 IF( wantu ) THEN
719 j = itgkz
720 DO i = 1, ns
721 CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
722 j = j + m*2
723 END DO
724*
725* Call DORMBR to compute QB*UB.
726* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
727*
728 CALL dormbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
729 $ work( itauq ), u, ldu, work( itemp ),
730 $ lwork-itemp+1, info )
731 END IF
732*
733* If needed, compute right singular vectors.
734*
735 IF( wantvt) THEN
736 j = itgkz + m
737 DO i = 1, ns
738 CALL dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
739 j = j + m*2
740 END DO
741 CALL dlaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ),
742 $ ldvt)
743*
744* Call DORMBR to compute (VB**T)*(PB**T)
745* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
746*
747 CALL dormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
748 $ work( itaup ), vt, ldvt, work( itemp ),
749 $ lwork-itemp+1, info )
750*
751* Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
752* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
753*
754 CALL dormlq( 'R', 'N', ns, n, m, a, lda,
755 $ work( itau ), vt, ldvt, work( itemp ),
756 $ lwork-itemp+1, info )
757 END IF
758 ELSE
759*
760* Path 2t (N greater than M, but not much larger)
761* Reduce to bidiagonal form without LQ decomposition
762* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
763* U = QB * UB; V**T = VB**T * PB**T
764*
765* Bidiagonalize A
766* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
767*
768 id = 1
769 ie = id + m
770 itauq = ie + m
771 itaup = itauq + m
772 itemp = itaup + m
773 CALL dgebrd( m, n, a, lda, work( id ), work( ie ),
774 $ work( itauq ), work( itaup ), work( itemp ),
775 $ lwork-itemp+1, info )
776*
777* Solve eigenvalue problem TGK*Z=Z*S.
778* (Workspace: need 2*M*M+14*M)
779*
780 itgkz = itemp
781 itemp = itgkz + m*(m*2+1)
782 CALL dbdsvdx( 'L', jobz, rngtgk, m, work( id ),
783 $ work( ie ),
784 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
785 $ m*2, work( itemp ), iwork, info)
786*
787* If needed, compute left singular vectors.
788*
789 IF( wantu ) THEN
790 j = itgkz
791 DO i = 1, ns
792 CALL dcopy( m, work( j ), 1, u( 1,i ), 1 )
793 j = j + m*2
794 END DO
795*
796* Call DORMBR to compute QB*UB.
797* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
798*
799 CALL dormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
800 $ work( itauq ), u, ldu, work( itemp ),
801 $ lwork-itemp+1, info )
802 END IF
803*
804* If needed, compute right singular vectors.
805*
806 IF( wantvt) THEN
807 j = itgkz + m
808 DO i = 1, ns
809 CALL dcopy( m, work( j ), 1, vt( i,1 ), ldvt )
810 j = j + m*2
811 END DO
812 CALL dlaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ),
813 $ ldvt)
814*
815* Call DORMBR to compute VB**T * PB**T
816* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
817*
818 CALL dormbr( 'P', 'R', 'T', ns, n, m, a, lda,
819 $ work( itaup ), vt, ldvt, work( itemp ),
820 $ lwork-itemp+1, info )
821 END IF
822 END IF
823 END IF
824*
825* Undo scaling if necessary
826*
827 IF( iscl.EQ.1 ) THEN
828 IF( anrm.GT.bignum )
829 $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1,
830 $ s, minmn, info )
831 IF( anrm.LT.smlnum )
832 $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
833 $ s, minmn, info )
834 END IF
835*
836* Return optimal workspace in WORK(1)
837*
838 work( 1 ) = dble( maxwrk )
839*
840 RETURN
841*
842* End of DGESVDX
843*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dbdsvdx
subroutine dbdsvdx(uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
DBDSVDX
Definition dbdsvdx.f:224
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgebrd
subroutine dgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
DGEBRD
Definition dgebrd.f:204
dgelqf
subroutine dgelqf(m, n, a, lda, tau, work, lwork, info)
DGELQF
Definition dgelqf.f:142
dgeqrf
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:144
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:112
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:142
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dormbr
subroutine dormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMBR
Definition dormbr.f:193
dormlq
subroutine dormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMLQ
Definition dormlq.f:165
dormqr
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:165
Here is the call graph for this function:
Here is the caller graph for this function: