Annotation of rpl/lapack/lapack/dgesvdx.f, revision 1.1
1.1 ! bertrand 1: *> \brief <b> DGESVDX computes the singular value decomposition (SVD) for GE matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGESVDX + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdx.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdx.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdx.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
! 22: * $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
! 23: * $ LWORK, IWORK, INFO )
! 24: *
! 25: *
! 26: * .. Scalar Arguments ..
! 27: * CHARACTER JOBU, JOBVT, RANGE
! 28: * INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
! 29: * DOUBLE PRECISION VL, VU
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * INTEGER IWORK( * )
! 33: * DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
! 34: * $ VT( LDVT, * ), WORK( * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: *> DGESVDX computes the singular value decomposition (SVD) of a real
! 44: *> M-by-N matrix A, optionally computing the left and/or right singular
! 45: *> vectors. The SVD is written
! 46: *>
! 47: *> A = U * SIGMA * transpose(V)
! 48: *>
! 49: *> where SIGMA is an M-by-N matrix which is zero except for its
! 50: *> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
! 51: *> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
! 52: *> are the singular values of A; they are real and non-negative, and
! 53: *> are returned in descending order. The first min(m,n) columns of
! 54: *> U and V are the left and right singular vectors of A.
! 55: *>
! 56: *> DGESVDX uses an eigenvalue problem for obtaining the SVD, which
! 57: *> allows for the computation of a subset of singular values and
! 58: *> vectors. See DBDSVDX for details.
! 59: *>
! 60: *> Note that the routine returns V**T, not V.
! 61: *> \endverbatim
! 62: *
! 63: * Arguments:
! 64: * ==========
! 65: *
! 66: *> \param[in] JOBU
! 67: *> \verbatim
! 68: *> JOBU is CHARACTER*1
! 69: *> Specifies options for computing all or part of the matrix U:
! 70: *> = 'V': the first min(m,n) columns of U (the left singular
! 71: *> vectors) or as specified by RANGE are returned in
! 72: *> the array U;
! 73: *> = 'N': no columns of U (no left singular vectors) are
! 74: *> computed.
! 75: *> \endverbatim
! 76: *>
! 77: *> \param[in] JOBVT
! 78: *> \verbatim
! 79: *> JOBVT is CHARACTER*1
! 80: *> Specifies options for computing all or part of the matrix
! 81: *> V**T:
! 82: *> = 'V': the first min(m,n) rows of V**T (the right singular
! 83: *> vectors) or as specified by RANGE are returned in
! 84: *> the array VT;
! 85: *> = 'N': no rows of V**T (no right singular vectors) are
! 86: *> computed.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] RANGE
! 90: *> \verbatim
! 91: *> RANGE is CHARACTER*1
! 92: *> = 'A': all singular values will be found.
! 93: *> = 'V': all singular values in the half-open interval (VL,VU]
! 94: *> will be found.
! 95: *> = 'I': the IL-th through IU-th singular values will be found.
! 96: *> \endverbatim
! 97: *>
! 98: *> \param[in] M
! 99: *> \verbatim
! 100: *> M is INTEGER
! 101: *> The number of rows of the input matrix A. M >= 0.
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in] N
! 105: *> \verbatim
! 106: *> N is INTEGER
! 107: *> The number of columns of the input matrix A. N >= 0.
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[in,out] A
! 111: *> \verbatim
! 112: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 113: *> On entry, the M-by-N matrix A.
! 114: *> On exit, the contents of A are destroyed.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[in] LDA
! 118: *> \verbatim
! 119: *> LDA is INTEGER
! 120: *> The leading dimension of the array A. LDA >= max(1,M).
! 121: *> \endverbatim
! 122: *>
! 123: *> \param[in] VL
! 124: *> \verbatim
! 125: *> VL is DOUBLE PRECISION
! 126: *> VL >=0.
! 127: *> \endverbatim
! 128: *>
! 129: *> \param[in] VU
! 130: *> \verbatim
! 131: *> VU is DOUBLE PRECISION
! 132: *> If RANGE='V', the lower and upper bounds of the interval to
! 133: *> be searched for singular values. VU > VL.
! 134: *> Not referenced if RANGE = 'A' or 'I'.
! 135: *> \endverbatim
! 136: *>
! 137: *> \param[in] IL
! 138: *> \verbatim
! 139: *> IL is INTEGER
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[in] IU
! 143: *> \verbatim
! 144: *> IU is INTEGER
! 145: *> If RANGE='I', the indices (in ascending order) of the
! 146: *> smallest and largest singular values to be returned.
! 147: *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
! 148: *> Not referenced if RANGE = 'A' or 'V'.
! 149: *> \endverbatim
! 150: *>
! 151: *> \param[out] NS
! 152: *> \verbatim
! 153: *> NS is INTEGER
! 154: *> The total number of singular values found,
! 155: *> 0 <= NS <= min(M,N).
! 156: *> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[out] S
! 160: *> \verbatim
! 161: *> S is DOUBLE PRECISION array, dimension (min(M,N))
! 162: *> The singular values of A, sorted so that S(i) >= S(i+1).
! 163: *> \endverbatim
! 164: *>
! 165: *> \param[out] U
! 166: *> \verbatim
! 167: *> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
! 168: *> If JOBU = 'V', U contains columns of U (the left singular
! 169: *> vectors, stored columnwise) as specified by RANGE; if
! 170: *> JOBU = 'N', U is not referenced.
! 171: *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
! 172: *> the exact value of NS is not known ILQFin advance and an upper
! 173: *> bound must be used.
! 174: *> \endverbatim
! 175: *>
! 176: *> \param[in] LDU
! 177: *> \verbatim
! 178: *> LDU is INTEGER
! 179: *> The leading dimension of the array U. LDU >= 1; if
! 180: *> JOBU = 'V', LDU >= M.
! 181: *> \endverbatim
! 182: *>
! 183: *> \param[out] VT
! 184: *> \verbatim
! 185: *> VT is DOUBLE PRECISION array, dimension (LDVT,N)
! 186: *> If JOBVT = 'V', VT contains the rows of V**T (the right singular
! 187: *> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
! 188: *> VT is not referenced.
! 189: *> Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
! 190: *> the exact value of NS is not known in advance and an upper
! 191: *> bound must be used.
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[in] LDVT
! 195: *> \verbatim
! 196: *> LDVT is INTEGER
! 197: *> The leading dimension of the array VT. LDVT >= 1; if
! 198: *> JOBVT = 'V', LDVT >= NS (see above).
! 199: *> \endverbatim
! 200: *>
! 201: *> \param[out] WORK
! 202: *> \verbatim
! 203: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 204: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
! 205: *> \endverbatim
! 206: *>
! 207: *> \param[in] LWORK
! 208: *> \verbatim
! 209: *> LWORK is INTEGER
! 210: *> The dimension of the array WORK.
! 211: *> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
! 212: *> comments inside the code):
! 213: *> - PATH 1 (M much larger than N)
! 214: *> - PATH 1t (N much larger than M)
! 215: *> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
! 216: *> For good performance, LWORK should generally be larger.
! 217: *>
! 218: *> If LWORK = -1, then a workspace query is assumed; the routine
! 219: *> only calculates the optimal size of the WORK array, returns
! 220: *> this value as the first entry of the WORK array, and no error
! 221: *> message related to LWORK is issued by XERBLA.
! 222: *> \endverbatim
! 223: *>
! 224: *> \param[out] IWORK
! 225: *> \verbatim
! 226: *> IWORK is INTEGER array, dimension (12*MIN(M,N))
! 227: *> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
! 228: *> then IWORK contains the indices of the eigenvectors that failed
! 229: *> to converge in DBDSVDX/DSTEVX.
! 230: *> \endverbatim
! 231: *>
! 232: *> \param[out] INFO
! 233: *> \verbatim
! 234: *> INFO is INTEGER
! 235: *> = 0: successful exit
! 236: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 237: *> > 0: if INFO = i, then i eigenvectors failed to converge
! 238: *> in DBDSVDX/DSTEVX.
! 239: *> if INFO = N*2 + 1, an internal error occurred in
! 240: *> DBDSVDX
! 241: *> \endverbatim
! 242: *
! 243: * Authors:
! 244: * ========
! 245: *
! 246: *> \author Univ. of Tennessee
! 247: *> \author Univ. of California Berkeley
! 248: *> \author Univ. of Colorado Denver
! 249: *> \author NAG Ltd.
! 250: *
! 251: *> \date November 2015
! 252: *
! 253: *> \ingroup doubleGEsing
! 254: *
! 255: * =====================================================================
! 256: SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
! 257: $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
! 258: $ LWORK, IWORK, INFO )
! 259: *
! 260: * -- LAPACK driver routine (version 3.6.0) --
! 261: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 262: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 263: * November 2015
! 264: *
! 265: * .. Scalar Arguments ..
! 266: CHARACTER JOBU, JOBVT, RANGE
! 267: INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
! 268: DOUBLE PRECISION VL, VU
! 269: * ..
! 270: * .. Array Arguments ..
! 271: INTEGER IWORK( * )
! 272: DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
! 273: $ VT( LDVT, * ), WORK( * )
! 274: * ..
! 275: *
! 276: * =====================================================================
! 277: *
! 278: * .. Parameters ..
! 279: DOUBLE PRECISION ZERO, ONE
! 280: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 281: * ..
! 282: * .. Local Scalars ..
! 283: CHARACTER JOBZ, RNGTGK
! 284: LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
! 285: INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
! 286: $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
! 287: $ J, MAXWRK, MINMN, MINWRK, MNTHR
! 288: DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
! 289: * ..
! 290: * .. Local Arrays ..
! 291: DOUBLE PRECISION DUM( 1 )
! 292: * ..
! 293: * .. External Subroutines ..
! 294: EXTERNAL DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY,
! 295: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR,
! 296: $ DSCAL, XERBLA
! 297: * ..
! 298: * .. External Functions ..
! 299: LOGICAL LSAME
! 300: INTEGER ILAENV
! 301: DOUBLE PRECISION DLAMCH, DLANGE, DNRM2
! 302: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE, DNRM2
! 303: * ..
! 304: * .. Intrinsic Functions ..
! 305: INTRINSIC MAX, MIN, SQRT
! 306: * ..
! 307: * .. Executable Statements ..
! 308: *
! 309: * Test the input arguments.
! 310: *
! 311: NS = 0
! 312: INFO = 0
! 313: ABSTOL = 2*DLAMCH('S')
! 314: LQUERY = ( LWORK.EQ.-1 )
! 315: MINMN = MIN( M, N )
! 316:
! 317: WANTU = LSAME( JOBU, 'V' )
! 318: WANTVT = LSAME( JOBVT, 'V' )
! 319: IF( WANTU .OR. WANTVT ) THEN
! 320: JOBZ = 'V'
! 321: ELSE
! 322: JOBZ = 'N'
! 323: END IF
! 324: ALLS = LSAME( RANGE, 'A' )
! 325: VALS = LSAME( RANGE, 'V' )
! 326: INDS = LSAME( RANGE, 'I' )
! 327: *
! 328: INFO = 0
! 329: IF( .NOT.LSAME( JOBU, 'V' ) .AND.
! 330: $ .NOT.LSAME( JOBU, 'N' ) ) THEN
! 331: INFO = -1
! 332: ELSE IF( .NOT.LSAME( JOBVT, 'V' ) .AND.
! 333: $ .NOT.LSAME( JOBVT, 'N' ) ) THEN
! 334: INFO = -2
! 335: ELSE IF( .NOT.( ALLS .OR. VALS .OR. INDS ) ) THEN
! 336: INFO = -3
! 337: ELSE IF( M.LT.0 ) THEN
! 338: INFO = -4
! 339: ELSE IF( N.LT.0 ) THEN
! 340: INFO = -5
! 341: ELSE IF( M.GT.LDA ) THEN
! 342: INFO = -7
! 343: ELSE IF( MINMN.GT.0 ) THEN
! 344: IF( VALS ) THEN
! 345: IF( VL.LT.ZERO ) THEN
! 346: INFO = -8
! 347: ELSE IF( VU.LE.VL ) THEN
! 348: INFO = -9
! 349: END IF
! 350: ELSE IF( INDS ) THEN
! 351: IF( IL.LT.1 .OR. IL.GT.MAX( 1, MINMN ) ) THEN
! 352: INFO = -10
! 353: ELSE IF( IU.LT.MIN( MINMN, IL ) .OR. IU.GT.MINMN ) THEN
! 354: INFO = -11
! 355: END IF
! 356: END IF
! 357: IF( INFO.EQ.0 ) THEN
! 358: IF( WANTU .AND. LDU.LT.M ) THEN
! 359: INFO = -15
! 360: ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
! 361: INFO = -16
! 362: END IF
! 363: END IF
! 364: END IF
! 365: *
! 366: * Compute workspace
! 367: * (Note: Comments in the code beginning "Workspace:" describe the
! 368: * minimal amount of workspace needed at that point in the code,
! 369: * as well as the preferred amount for good performance.
! 370: * NB refers to the optimal block size for the immediately
! 371: * following subroutine, as returned by ILAENV.)
! 372: *
! 373: IF( INFO.EQ.0 ) THEN
! 374: MINWRK = 1
! 375: MAXWRK = 1
! 376: IF( MINMN.GT.0 ) THEN
! 377: IF( M.GE.N ) THEN
! 378: MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
! 379: IF( M.GE.MNTHR ) THEN
! 380: *
! 381: * Path 1 (M much larger than N)
! 382: *
! 383: MAXWRK = N*(N*2+16) +
! 384: $ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
! 385: MAXWRK = MAX( MAXWRK, N*(N*2+20) + 2*N*
! 386: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
! 387: MINWRK = N*(N*2+21)
! 388: ELSE
! 389: *
! 390: * Path 2 (M at least N, but not much larger)
! 391: *
! 392: MAXWRK = N*(N*2+19) + ( M+N )*
! 393: $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
! 394: MINWRK = N*(N*2+20) + M
! 395: END IF
! 396: ELSE
! 397: MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
! 398: IF( N.GE.MNTHR ) THEN
! 399: *
! 400: * Path 1t (N much larger than M)
! 401: *
! 402: MAXWRK = M*(M*2+16) +
! 403: $ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
! 404: MAXWRK = MAX( MAXWRK, M*(M*2+20) + 2*M*
! 405: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
! 406: MINWRK = M*(M*2+21)
! 407: ELSE
! 408: *
! 409: * Path 2t (N greater than M, but not much larger)
! 410: *
! 411: MAXWRK = M*(M*2+19) + ( M+N )*
! 412: $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
! 413: MINWRK = M*(M*2+20) + N
! 414: END IF
! 415: END IF
! 416: END IF
! 417: MAXWRK = MAX( MAXWRK, MINWRK )
! 418: WORK( 1 ) = DBLE( MAXWRK )
! 419: *
! 420: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 421: INFO = -19
! 422: END IF
! 423: END IF
! 424: *
! 425: IF( INFO.NE.0 ) THEN
! 426: CALL XERBLA( 'DGESVDX', -INFO )
! 427: RETURN
! 428: ELSE IF( LQUERY ) THEN
! 429: RETURN
! 430: END IF
! 431: *
! 432: * Quick return if possible
! 433: *
! 434: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 435: RETURN
! 436: END IF
! 437: *
! 438: * Set singular values indices accord to RANGE.
! 439: *
! 440: IF( ALLS ) THEN
! 441: RNGTGK = 'I'
! 442: ILTGK = 1
! 443: IUTGK = MIN( M, N )
! 444: ELSE IF( INDS ) THEN
! 445: RNGTGK = 'I'
! 446: ILTGK = IL
! 447: IUTGK = IU
! 448: ELSE
! 449: RNGTGK = 'V'
! 450: ILTGK = 0
! 451: IUTGK = 0
! 452: END IF
! 453: *
! 454: * Get machine constants
! 455: *
! 456: EPS = DLAMCH( 'P' )
! 457: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
! 458: BIGNUM = ONE / SMLNUM
! 459: *
! 460: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 461: *
! 462: ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
! 463: ISCL = 0
! 464: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 465: ISCL = 1
! 466: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 467: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 468: ISCL = 1
! 469: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 470: END IF
! 471: *
! 472: IF( M.GE.N ) THEN
! 473: *
! 474: * A has at least as many rows as columns. If A has sufficiently
! 475: * more rows than columns, first reduce A using the QR
! 476: * decomposition.
! 477: *
! 478: IF( M.GE.MNTHR ) THEN
! 479: *
! 480: * Path 1 (M much larger than N):
! 481: * A = Q * R = Q * ( QB * B * PB**T )
! 482: * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
! 483: * U = Q * QB * UB; V**T = VB**T * PB**T
! 484: *
! 485: * Compute A=Q*R
! 486: * (Workspace: need 2*N, prefer N+N*NB)
! 487: *
! 488: ITAU = 1
! 489: ITEMP = ITAU + N
! 490: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
! 491: $ LWORK-ITEMP+1, INFO )
! 492: *
! 493: * Copy R into WORK and bidiagonalize it:
! 494: * (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
! 495: *
! 496: IQRF = ITEMP
! 497: ID = IQRF + N*N
! 498: IE = ID + N
! 499: ITAUQ = IE + N
! 500: ITAUP = ITAUQ + N
! 501: ITEMP = ITAUP + N
! 502: CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
! 503: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N )
! 504: CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),
! 505: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
! 506: $ LWORK-ITEMP+1, INFO )
! 507: *
! 508: * Solve eigenvalue problem TGK*Z=Z*S.
! 509: * (Workspace: need 14*N + 2*N*(N+1))
! 510: *
! 511: ITGKZ = ITEMP
! 512: ITEMP = ITGKZ + N*(N*2+1)
! 513: CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
! 514: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
! 515: $ N*2, WORK( ITEMP ), IWORK, INFO)
! 516: *
! 517: * If needed, compute left singular vectors.
! 518: *
! 519: IF( WANTU ) THEN
! 520: J = ITGKZ
! 521: DO I = 1, NS
! 522: CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
! 523: J = J + N*2
! 524: END DO
! 525: CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
! 526: *
! 527: * Call DORMBR to compute QB*UB.
! 528: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
! 529: *
! 530: CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
! 531: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
! 532: $ LWORK-ITEMP+1, INFO )
! 533: *
! 534: * Call DORMQR to compute Q*(QB*UB).
! 535: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
! 536: *
! 537: CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,
! 538: $ WORK( ITAU ), U, LDU, WORK( ITEMP ),
! 539: $ LWORK-ITEMP+1, INFO )
! 540: END IF
! 541: *
! 542: * If needed, compute right singular vectors.
! 543: *
! 544: IF( WANTVT) THEN
! 545: J = ITGKZ + N
! 546: DO I = 1, NS
! 547: CALL DCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT )
! 548: J = J + N*2
! 549: END DO
! 550: *
! 551: * Call DORMBR to compute VB**T * PB**T
! 552: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
! 553: *
! 554: CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,
! 555: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
! 556: $ LWORK-ITEMP+1, INFO )
! 557: END IF
! 558: ELSE
! 559: *
! 560: * Path 2 (M at least N, but not much larger)
! 561: * Reduce A to bidiagonal form without QR decomposition
! 562: * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
! 563: * U = QB * UB; V**T = VB**T * PB**T
! 564: *
! 565: * Bidiagonalize A
! 566: * (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
! 567: *
! 568: ID = 1
! 569: IE = ID + N
! 570: ITAUQ = IE + N
! 571: ITAUP = ITAUQ + N
! 572: ITEMP = ITAUP + N
! 573: CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
! 574: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
! 575: $ LWORK-ITEMP+1, INFO )
! 576: *
! 577: * Solve eigenvalue problem TGK*Z=Z*S.
! 578: * (Workspace: need 14*N + 2*N*(N+1))
! 579: *
! 580: ITGKZ = ITEMP
! 581: ITEMP = ITGKZ + N*(N*2+1)
! 582: CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
! 583: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
! 584: $ N*2, WORK( ITEMP ), IWORK, INFO)
! 585: *
! 586: * If needed, compute left singular vectors.
! 587: *
! 588: IF( WANTU ) THEN
! 589: J = ITGKZ
! 590: DO I = 1, NS
! 591: CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
! 592: J = J + N*2
! 593: END DO
! 594: CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
! 595: *
! 596: * Call DORMBR to compute QB*UB.
! 597: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
! 598: *
! 599: CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
! 600: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
! 601: $ LWORK-ITEMP+1, IERR )
! 602: END IF
! 603: *
! 604: * If needed, compute right singular vectors.
! 605: *
! 606: IF( WANTVT) THEN
! 607: J = ITGKZ + N
! 608: DO I = 1, NS
! 609: CALL DCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT )
! 610: J = J + N*2
! 611: END DO
! 612: *
! 613: * Call DORMBR to compute VB**T * PB**T
! 614: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
! 615: *
! 616: CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,
! 617: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
! 618: $ LWORK-ITEMP+1, IERR )
! 619: END IF
! 620: END IF
! 621: ELSE
! 622: *
! 623: * A has more columns than rows. If A has sufficiently more
! 624: * columns than rows, first reduce A using the LQ decomposition.
! 625: *
! 626: IF( N.GE.MNTHR ) THEN
! 627: *
! 628: * Path 1t (N much larger than M):
! 629: * A = L * Q = ( QB * B * PB**T ) * Q
! 630: * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
! 631: * U = QB * UB ; V**T = VB**T * PB**T * Q
! 632: *
! 633: * Compute A=L*Q
! 634: * (Workspace: need 2*M, prefer M+M*NB)
! 635: *
! 636: ITAU = 1
! 637: ITEMP = ITAU + M
! 638: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
! 639: $ LWORK-ITEMP+1, INFO )
! 640:
! 641: * Copy L into WORK and bidiagonalize it:
! 642: * (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
! 643: *
! 644: ILQF = ITEMP
! 645: ID = ILQF + M*M
! 646: IE = ID + M
! 647: ITAUQ = IE + M
! 648: ITAUP = ITAUQ + M
! 649: ITEMP = ITAUP + M
! 650: CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
! 651: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M )
! 652: CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),
! 653: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
! 654: $ LWORK-ITEMP+1, INFO )
! 655: *
! 656: * Solve eigenvalue problem TGK*Z=Z*S.
! 657: * (Workspace: need 2*M*M+14*M)
! 658: *
! 659: ITGKZ = ITEMP
! 660: ITEMP = ITGKZ + M*(M*2+1)
! 661: CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
! 662: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
! 663: $ M*2, WORK( ITEMP ), IWORK, INFO)
! 664: *
! 665: * If needed, compute left singular vectors.
! 666: *
! 667: IF( WANTU ) THEN
! 668: J = ITGKZ
! 669: DO I = 1, NS
! 670: CALL DCOPY( M, WORK( J ), 1, U( 1,I ), 1 )
! 671: J = J + M*2
! 672: END DO
! 673: *
! 674: * Call DORMBR to compute QB*UB.
! 675: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
! 676: *
! 677: CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
! 678: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
! 679: $ LWORK-ITEMP+1, INFO )
! 680: END IF
! 681: *
! 682: * If needed, compute right singular vectors.
! 683: *
! 684: IF( WANTVT) THEN
! 685: J = ITGKZ + M
! 686: DO I = 1, NS
! 687: CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
! 688: J = J + M*2
! 689: END DO
! 690: CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
! 691: *
! 692: * Call DORMBR to compute (VB**T)*(PB**T)
! 693: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
! 694: *
! 695: CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,
! 696: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
! 697: $ LWORK-ITEMP+1, INFO )
! 698: *
! 699: * Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
! 700: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
! 701: *
! 702: CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,
! 703: $ WORK( ITAU ), VT, LDVT, WORK( ITEMP ),
! 704: $ LWORK-ITEMP+1, INFO )
! 705: END IF
! 706: ELSE
! 707: *
! 708: * Path 2t (N greater than M, but not much larger)
! 709: * Reduce to bidiagonal form without LQ decomposition
! 710: * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
! 711: * U = QB * UB; V**T = VB**T * PB**T
! 712: *
! 713: * Bidiagonalize A
! 714: * (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
! 715: *
! 716: ID = 1
! 717: IE = ID + M
! 718: ITAUQ = IE + M
! 719: ITAUP = ITAUQ + M
! 720: ITEMP = ITAUP + M
! 721: CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
! 722: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
! 723: $ LWORK-ITEMP+1, INFO )
! 724: *
! 725: * Solve eigenvalue problem TGK*Z=Z*S.
! 726: * (Workspace: need 2*M*M+14*M)
! 727: *
! 728: ITGKZ = ITEMP
! 729: ITEMP = ITGKZ + M*(M*2+1)
! 730: CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
! 731: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
! 732: $ M*2, WORK( ITEMP ), IWORK, INFO)
! 733: *
! 734: * If needed, compute left singular vectors.
! 735: *
! 736: IF( WANTU ) THEN
! 737: J = ITGKZ
! 738: DO I = 1, NS
! 739: CALL DCOPY( M, WORK( J ), 1, U( 1,I ), 1 )
! 740: J = J + M*2
! 741: END DO
! 742: *
! 743: * Call DORMBR to compute QB*UB.
! 744: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
! 745: *
! 746: CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
! 747: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
! 748: $ LWORK-ITEMP+1, INFO )
! 749: END IF
! 750: *
! 751: * If needed, compute right singular vectors.
! 752: *
! 753: IF( WANTVT) THEN
! 754: J = ITGKZ + M
! 755: DO I = 1, NS
! 756: CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
! 757: J = J + M*2
! 758: END DO
! 759: CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
! 760: *
! 761: * Call DORMBR to compute VB**T * PB**T
! 762: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
! 763: *
! 764: CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,
! 765: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
! 766: $ LWORK-ITEMP+1, INFO )
! 767: END IF
! 768: END IF
! 769: END IF
! 770: *
! 771: * Undo scaling if necessary
! 772: *
! 773: IF( ISCL.EQ.1 ) THEN
! 774: IF( ANRM.GT.BIGNUM )
! 775: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1,
! 776: $ S, MINMN, INFO )
! 777: IF( ANRM.LT.SMLNUM )
! 778: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1,
! 779: $ S, MINMN, INFO )
! 780: END IF
! 781: *
! 782: * Return optimal workspace in WORK(1)
! 783: *
! 784: WORK( 1 ) = DBLE( MAXWRK )
! 785: *
! 786: RETURN
! 787: *
! 788: * End of DGESVDX
! 789: *
! 790: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgesvdx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack