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