SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, \$ LWORK, IWORK, INFO ) * * -- LAPACK driver routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * 8-15-00: Improve consistency of WS calculations (eca) * * .. Scalar Arguments .. CHARACTER JOBZ INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N * .. * .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), \$ VT( LDVT, * ), WORK( * ) * .. * * Purpose * ======= * * DGESDD computes the singular value decomposition (SVD) of a real * M-by-N matrix A, optionally computing the left and right singular * vectors. If singular vectors are desired, it uses a * divide-and-conquer algorithm. * * 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. * * Note that the routine returns VT = V**T, not V. * * The divide and conquer algorithm makes very mild assumptions about * floating point arithmetic. It will work on machines with a guard * digit in add/subtract, or on those binary machines without guard * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or * Cray-2. It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. * * Arguments * ========= * * JOBZ (input) CHARACTER*1 * Specifies options for computing all or part of the matrix U: * = 'A': all M columns of U and all N rows of V**T are * returned in the arrays U and VT; * = 'S': the first min(M,N) columns of U and the first * min(M,N) rows of V**T are returned in the arrays U * and VT; * = 'O': If M >= N, the first N columns of U are overwritten * on the array A and all rows of V**T are returned in * the array VT; * otherwise, all columns of U are returned in the * array U and the first M rows of V**T are overwritten * in the array VT; * = 'N': no columns of U or rows of V**T are computed. * * M (input) INTEGER * The number of rows of the input matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, * if JOBZ = 'O', A is overwritten with the first N columns * of U (the left singular vectors, stored * columnwise) if M >= N; * A is overwritten with the first M rows * of V**T (the right singular vectors, stored * rowwise) otherwise. * if JOBZ .ne. 'O', the contents of A are destroyed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * S (output) DOUBLE PRECISION array, dimension (min(M,N)) * The singular values of A, sorted so that S(i) >= S(i+1). * * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) * UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; * UCOL = min(M,N) if JOBZ = 'S'. * If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M * orthogonal matrix U; * if JOBZ = 'S', U contains the first min(M,N) columns of U * (the left singular vectors, stored columnwise); * if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. * * LDU (input) INTEGER * The leading dimension of the array U. LDU >= 1; if * JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. * * VT (output) DOUBLE PRECISION array, dimension (LDVT,N) * If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the * N-by-N orthogonal matrix V**T; * if JOBZ = 'S', VT contains the first min(M,N) rows of * V**T (the right singular vectors, stored rowwise); * if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. * * LDVT (input) INTEGER * The leading dimension of the array VT. LDVT >= 1; if * JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; * if JOBZ = 'S', LDVT >= min(M,N). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK; * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= 1. * If JOBZ = 'N', * LWORK >= max(14*min(M,N)+4, 10*min(M,N)+2+ * SMLSIZ*(SMLSIZ+8)) + max(M,N) * where SMLSIZ is returned by ILAENV and is equal to the * maximum size of the subproblems at the bottom of the * computation tree (usually about 25). * If JOBZ = 'O', * LWORK >= 5*min(M,N)*min(M,N) + max(M,N) + 9*min(M,N). * If JOBZ = 'S' or 'A' * LWORK >= 4*min(M,N)*min(M,N) + max(M,N) + 9*min(M,N). * For good performance, LWORK should generally be larger. * * If LWORK = -1, a workspace query is assumed. The optimal * size for the WORK array is calculated and stored in WORK(1), * and no other work except argument checking is performed. * * IWORK (workspace) INTEGER array, dimension (8*min(M,N)) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: DBDSDC did not converge, updating process failed. * * Further Details * =============== * * Based on contributions by * Ming Gu and Huan Ren, Computer Science Division, University of * California at Berkeley, USA * * ===================================================================== * * .. Parameters .. INTEGER LQUERV PARAMETER ( LQUERV = -1 ) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS INTEGER BDSPAC, BDSPAN, BLK, CHUNK, I, IE, IERR, IL, \$ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, \$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, \$ MNTHR, NWORK, SMLSIZ, WRKBL DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM * .. * .. Local Arrays .. INTEGER IDUM( 1 ) DOUBLE PRECISION DUM( 1 ) * .. * .. External Subroutines .. EXTERNAL DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY, \$ DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR, \$ XERBLA * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLANGE EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE * .. * .. Intrinsic Functions .. INTRINSIC DBLE, INT, MAX, MIN, SQRT * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 MINMN = MIN( M, N ) MNTHR = INT( MINMN*11.0D0 / 6.0D0 ) WNTQA = LSAME( JOBZ, 'A' ) WNTQS = LSAME( JOBZ, 'S' ) WNTQAS = WNTQA .OR. WNTQS WNTQO = LSAME( JOBZ, 'O' ) WNTQN = LSAME( JOBZ, 'N' ) MINWRK = 1 MAXWRK = 1 * IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -5 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR. \$ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN INFO = -8 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR. \$ ( WNTQS .AND. LDVT.LT.MINMN ) .OR. \$ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN INFO = -10 END IF * SMLSIZ = ILAENV( 9, 'DGESDD', ' ', 0, 0, 0, 0 ) * * Compute workspace * (Note: Comments in the code beginning "Workspace:" describe the * minimal amount of workspace needed at that point in the code, * as well as the preferred amount for good performance. * NB refers to the optimal block size for the immediately * following subroutine, as returned by ILAENV.) * IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN IF( M.GE.N ) THEN * * Compute space needed for DBDSDC * BDSPAC = 3*N*N + 7*N BDSPAN = MAX( 12*N+4, 8*N+2+SMLSIZ*( SMLSIZ+8 ) ) IF( M.GE.MNTHR ) THEN IF( WNTQN ) THEN * * Path 1 (M much larger than N, JOBZ='N') * MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, \$ -1 ) MAXWRK = MAX( MAXWRK, 3*N+2*N* \$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = BDSPAC ELSE IF( WNTQO ) THEN * * Path 2 (M much larger than N, JOBZ='O') * WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, \$ N, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+2*N* \$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*N ) MAXWRK = WRKBL + 2*N*N MINWRK = BDSPAC + 2*N*N + 2*N ELSE IF( WNTQS ) THEN * * Path 3 (M much larger than N, JOBZ='S') * WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, \$ N, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+2*N* \$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*N ) MAXWRK = WRKBL + N*N MINWRK = BDSPAC + N*N + 2*N ELSE IF( WNTQA ) THEN * * Path 4 (M much larger than N, JOBZ='A') * WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, \$ M, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+2*N* \$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*N ) MAXWRK = N*N + WRKBL MINWRK = BDSPAC + N*N + M + N END IF ELSE * * Path 5 (M at least N, but not much larger) * WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1, \$ -1 ) IF( WNTQO ) THEN WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) ) WRKBL = MAX( WRKBL, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*N+M ) MAXWRK = WRKBL + M*N MINWRK = BDSPAC + N*N + 2*N + M ELSE IF( WNTQS ) THEN MAXWRK = MAX( MAXWRK, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) ) MAXWRK = MAX( MAXWRK, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) MAXWRK = MAX( MAXWRK, BDSPAC+2*N+M ) MINWRK = BDSPAC + 2*N + M ELSE IF( WNTQA ) THEN MAXWRK = MAX( MAXWRK, 3*N+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) MAXWRK = MAX( MAXWRK, 3*N+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) ) MAXWRK = MAX( MAXWRK, BDSPAC+2*N+M ) MINWRK = BDSPAC + 2*N + M END IF END IF ELSE * * Compute space needed for DBDSDC * BDSPAC = 3*M*M + 7*M BDSPAN = MAX( 12*M+4, 8*M+2+SMLSIZ*( SMLSIZ+8 ) ) IF( N.GE.MNTHR ) THEN IF( WNTQN ) THEN * * Path 1t (N much larger than M, JOBZ='N') * MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, \$ -1 ) MAXWRK = MAX( MAXWRK, 3*M+2*M* \$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = BDSPAC ELSE IF( WNTQO ) THEN * * Path 2t (N much larger than M, JOBZ='O') * WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, \$ N, M, -1 ) ) WRKBL = MAX( WRKBL, 3*M+2*M* \$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*M ) MAXWRK = WRKBL + 2*M*M MINWRK = BDSPAC + 2*M*M + 2*M ELSE IF( WNTQS ) THEN * * Path 3t (N much larger than M, JOBZ='S') * WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, \$ N, M, -1 ) ) WRKBL = MAX( WRKBL, 3*M+2*M* \$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*M ) MAXWRK = WRKBL + M*M MINWRK = BDSPAC + M*M + 2*M ELSE IF( WNTQA ) THEN * * Path 4t (N much larger than M, JOBZ='A') * WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, \$ N, M, -1 ) ) WRKBL = MAX( WRKBL, 3*M+2*M* \$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*M ) MAXWRK = WRKBL + M*M MINWRK = BDSPAC + M*M + M + N END IF ELSE * * Path 5t (N greater than M, but not much larger) * WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1, \$ -1 ) IF( WNTQO ) THEN WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) WRKBL = MAX( WRKBL, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) ) WRKBL = MAX( WRKBL, BDSPAC+2*M ) MAXWRK = WRKBL + M*N MINWRK = BDSPAC + M*M + 2*M + N ELSE IF( WNTQS ) THEN MAXWRK = MAX( MAXWRK, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) MAXWRK = MAX( MAXWRK, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) ) MAXWRK = MAX( MAXWRK, BDSPAC+2*M ) MINWRK = BDSPAC + 2*M + N ELSE IF( WNTQA ) THEN MAXWRK = MAX( MAXWRK, 3*M+M* \$ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) ) MAXWRK = MAX( MAXWRK, 3*M+N* \$ ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) ) MAXWRK = MAX( MAXWRK, BDSPAC+2*M ) MINWRK = BDSPAC + 2*M + N END IF END IF END IF END IF IF( INFO.EQ.0 ) THEN WORK( 1 ) = MAXWRK IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV ) \$ INFO = -12 END IF * * Quick returns * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGESDD', -INFO ) RETURN END IF IF( LWORK.EQ.LQUERV ) \$ RETURN IF( M.EQ.0 .OR. N.EQ.0 ) THEN RETURN END IF * * Get machine constants * EPS = DLAMCH( 'P' ) SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS BIGNUM = ONE / SMLNUM * * Scale A if max element outside range [SMLNUM,BIGNUM] * ANRM = DLANGE( 'M', M, N, A, LDA, DUM ) ISCL = 0 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN ISCL = 1 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) ELSE IF( ANRM.GT.BIGNUM ) THEN ISCL = 1 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) END IF * IF( M.GE.N ) THEN * * A has at least as many rows as columns. If A has sufficiently * more rows than columns, first reduce using the QR * decomposition (if sufficient workspace available) * IF( M.GE.MNTHR ) THEN * IF( WNTQN ) THEN * * Path 1 (M much larger than N, JOBZ='N') * No singular vectors to be computed * ITAU = 1 NWORK = ITAU + N * * Compute A=Q*R * (Workspace: need 2*N, prefer N+N*NB) * CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Zero out below R * CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) IE = 1 ITAUQ = IE + N ITAUP = ITAUQ + N NWORK = ITAUP + N * * Bidiagonalize R in A * (Workspace: need 4*N, prefer 3*N+2*N*NB) * CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), \$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, \$ IERR ) NWORK = IE + N * * Perform bidiagonal SVD, computing singular values only * (Workspace: need BDSPAN) * CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1, \$ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) * ELSE IF( WNTQO ) THEN * * Path 2 (M much larger than N, JOBZ = 'O') * N left singular vectors to be overwritten on A and * N right singular vectors to be computed in VT * IR = 1 * * WORK(IR) is LDWRKR by N * IF( LWORK.GE.LDA*N+4*N*N+9*N ) THEN LDWRKR = LDA ELSE LDWRKR = ( LWORK-4*N*N-9*N ) / N END IF ITAU = IR + LDWRKR*N NWORK = ITAU + N * * Compute A=Q*R * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Copy R to WORK(IR), zeroing out below it * CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ), \$ LDWRKR ) * * Generate Q in A * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) IE = ITAU ITAUQ = IE + N ITAUP = ITAUQ + N NWORK = ITAUP + N * * Bidiagonalize R in VT, copying result to WORK(IR) * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) * CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), \$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * WORK(IU) is N by N * IU = NWORK NWORK = IU + N*N * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in WORK(IU) and computing right * singular vectors of bidiagonal matrix in VT * (Workspace: need 2*N*N+BDSPAC) * CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N, \$ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Overwrite WORK(IU) by left singular vectors of R * and VT by right singular vectors of R * (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) * CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, \$ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Multiply Q in A by left singular vectors of R in * WORK(IU), storing result in WORK(IR) and copying to A * (Workspace: need 2*N*N, prefer N*N+M*N) * DO 10 I = 1, M, LDWRKR CHUNK = MIN( M-I+1, LDWRKR ) CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ), \$ LDA, WORK( IU ), N, ZERO, WORK( IR ), \$ LDWRKR ) CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, \$ A( I, 1 ), LDA ) 10 CONTINUE * ELSE IF( WNTQS ) THEN * * Path 3 (M much larger than N, JOBZ='S') * N left singular vectors to be computed in U and * N right singular vectors to be computed in VT * IR = 1 * * WORK(IR) is N by N * LDWRKR = N ITAU = IR + LDWRKR*N NWORK = ITAU + N * * Compute A=Q*R * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Copy R to WORK(IR), zeroing out below it * CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ), \$ LDWRKR ) * * Generate Q in A * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) IE = ITAU ITAUQ = IE + N ITAUP = ITAUQ + N NWORK = ITAUP + N * * Bidiagonalize R in WORK(IR) * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) * CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), \$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Perform bidiagonal SVD, computing left singular vectors * of bidiagoal matrix in U and computing right singular * vectors of bidiagonal matrix in VT * (Workspace: need N*N+BDSPAC) * CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, \$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Overwrite U by left singular vectors of R and VT * by right singular vectors of R * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) * CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Multiply Q in A by left singular vectors of R in * WORK(IR), storing result in U * (Workspace: need N*N) * CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ), \$ LDWRKR, ZERO, U, LDU ) * ELSE IF( WNTQA ) THEN * * Path 4 (M much larger than N, JOBZ='A') * M left singular vectors to be computed in U and * N right singular vectors to be computed in VT * IU = 1 * * WORK(IU) is N by N * LDWRKU = N ITAU = IU + LDWRKU*N NWORK = ITAU + N * * Compute A=Q*R, copying result to U * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) * * Generate Q in U * (Workspace: need N*N+N+M, prefer N*N+N+M*NB) CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Produce R in A, zeroing out other entries * CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) IE = ITAU ITAUQ = IE + N ITAUP = ITAUQ + N NWORK = ITAUP + N * * Bidiagonalize R in A * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) * CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), \$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, \$ IERR ) * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in WORK(IU) and computing right * singular vectors of bidiagonal matrix in VT * (Workspace: need N*N+BDSPAC) * CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N, \$ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Overwrite WORK(IU) by left singular vectors of R and VT * by right singular vectors of R * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB) * CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA, \$ WORK( ITAUQ ), WORK( IU ), LDWRKU, \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Multiply Q in U by left singular vectors of R in * WORK(IU), storing result in A * (Workspace: need N*N) * CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ), \$ LDWRKU, ZERO, A, LDA ) * * Copy left singular vectors of A from A to U * CALL DLACPY( 'F', M, N, A, LDA, U, LDU ) * END IF * ELSE * * M .LT. MNTHR * * Path 5 (M at least N, but not much larger) * Reduce to bidiagonal form without QR decomposition * IE = 1 ITAUQ = IE + N ITAUP = ITAUQ + N NWORK = ITAUP + N * * Bidiagonalize A * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) * CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), \$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, \$ IERR ) IF( WNTQN ) THEN * * Perform bidiagonal SVD, only computing singular values * (Workspace: need BDSPAN) * CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1, \$ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) ELSE IF( WNTQO ) THEN IU = NWORK IF( LWORK.GE.M*N+3*N*N+9*N ) THEN * * WORK( IU ) is M by N * LDWRKU = M NWORK = IU + LDWRKU*N CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ), \$ LDWRKU ) ELSE * * WORK( IU ) is N by N * LDWRKU = N NWORK = IU + LDWRKU*N * * WORK(IR) is LDWRKR by N * IR = NWORK LDWRKR = ( LWORK-N*N-3*N ) / N END IF NWORK = IU + LDWRKU*N * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in WORK(IU) and computing right * singular vectors of bidiagonal matrix in VT * (Workspace: need N*N+BDSPAC) * CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), \$ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ), \$ IWORK, INFO ) * * Overwrite VT by right singular vectors of A * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * IF( LWORK.GE.M*N+3*N*N+9*N ) THEN * * Overwrite WORK(IU) by left singular vectors of A * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA, \$ WORK( ITAUQ ), WORK( IU ), LDWRKU, \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Copy left singular vectors of A from WORK(IU) to A * CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) ELSE * * Generate Q in A * (Workspace: need N*N+2*N, prefer N*N+N+N*NB) * CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Multiply Q in A by left singular vectors of * bidiagonal matrix in WORK(IU), storing result in * WORK(IR) and copying to A * (Workspace: need 2*N*N, prefer N*N+M*N) * DO 20 I = 1, M, LDWRKR CHUNK = MIN( M-I+1, LDWRKR ) CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ), \$ LDA, WORK( IU ), LDWRKU, ZERO, \$ WORK( IR ), LDWRKR ) CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, \$ A( I, 1 ), LDA ) 20 CONTINUE END IF * ELSE IF( WNTQS ) THEN * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in VT * (Workspace: need BDSPAC) * CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU ) CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, \$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Overwrite U by left singular vectors of A and VT * by right singular vectors of A * (Workspace: need 3*N, prefer 2*N+N*NB) * CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) ELSE * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in VT * (Workspace: need BDSPAC) * CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU ) CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, \$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Set the right corner of U to identity matrix * CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ), \$ LDU ) * * Overwrite U by left singular vectors of A and VT * by right singular vectors of A * (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) END IF * END IF * ELSE * * A has more columns than rows. If A has sufficiently more * columns than rows, first reduce using the LQ decomposition (if * sufficient workspace available) * IF( N.GE.MNTHR ) THEN * IF( WNTQN ) THEN * * Path 1t (N much larger than M, JOBZ='N') * No singular vectors to be computed * ITAU = 1 NWORK = ITAU + M * * Compute A=L*Q * (Workspace: need 2*M, prefer M+M*NB) * CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Zero out above L * CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA ) IE = 1 ITAUQ = IE + M ITAUP = ITAUQ + M NWORK = ITAUP + M * * Bidiagonalize L in A * (Workspace: need 4*M, prefer 3*M+2*M*NB) * CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), \$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, \$ IERR ) NWORK = IE + M * * Perform bidiagonal SVD, computing singular values only * (Workspace: need BDSPAN) * CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1, \$ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) * ELSE IF( WNTQO ) THEN * * Path 2t (N much larger than M, JOBZ='O') * M right singular vectors to be overwritten on A and * M left singular vectors to be computed in U * IVT = 1 * * IVT is M by M * IL = IVT + M*M IF( LWORK.GE.M*N+4*M*M+9*M ) THEN * * WORK(IL) is M by N * LDWRKL = M CHUNK = N ELSE LDWRKL = M CHUNK = ( LWORK-M*M ) / M END IF ITAU = IL + LDWRKL*M NWORK = ITAU + M * * Compute A=L*Q * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Copy L to WORK(IL), zeroing about above it * CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, \$ WORK( IL+LDWRKL ), LDWRKL ) * * Generate Q in A * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) IE = ITAU ITAUQ = IE + M ITAUP = ITAUQ + M NWORK = ITAUP + M * * Bidiagonalize L in WORK(IL) * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) * CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ), \$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U, and computing right singular * vectors of bidiagonal matrix in WORK(IVT) * (Workspace: need 2*M*M+BDSPAC) * CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, \$ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ), \$ IWORK, INFO ) * * Overwrite U by left singular vectors of L and WORK(IVT) * by right singular vectors of L * (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL, \$ WORK( ITAUP ), WORK( IVT ), M, \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Multiply right singular vectors of L in WORK(IVT) by Q * in A, storing result in WORK(IL) and copying to A * (Workspace: need 2*M*M, prefer M*M+M*N) * DO 30 I = 1, N, CHUNK BLK = MIN( N-I+1, CHUNK ) CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M, \$ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL ) CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL, \$ A( 1, I ), LDA ) 30 CONTINUE * ELSE IF( WNTQS ) THEN * * Path 3t (N much larger than M, JOBZ='S') * M right singular vectors to be computed in VT and * M left singular vectors to be computed in U * IL = 1 * * WORK(IL) is M by M * LDWRKL = M ITAU = IL + LDWRKL*M NWORK = ITAU + M * * Compute A=L*Q * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Copy L to WORK(IL), zeroing out above it * CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, \$ WORK( IL+LDWRKL ), LDWRKL ) * * Generate Q in A * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) IE = ITAU ITAUQ = IE + M ITAUP = ITAUQ + M NWORK = ITAUP + M * * Bidiagonalize L in WORK(IU), copying result to U * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) * CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ), \$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in VT * (Workspace: need M*M+BDSPAC) * CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT, \$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Overwrite U by left singular vectors of L and VT * by right singular vectors of L * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * * Multiply right singular vectors of L in WORK(IL) by * Q in A, storing result in VT * (Workspace: need M*M) * CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL, \$ A, LDA, ZERO, VT, LDVT ) * ELSE IF( WNTQA ) THEN * * Path 4t (N much larger than M, JOBZ='A') * N right singular vectors to be computed in VT and * M left singular vectors to be computed in U * IVT = 1 * * WORK(IVT) is M by M * LDWKVT = M ITAU = IVT + LDWKVT*M NWORK = ITAU + M * * Compute A=L*Q, copying result to VT * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) * * Generate Q in VT * (Workspace: need M*M+M+N, prefer M*M+M+N*NB) * CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Produce L in A, zeroing out other entries * CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA ) IE = ITAU ITAUQ = IE + M ITAUP = ITAUQ + M NWORK = ITAUP + M * * Bidiagonalize L in A * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) * CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), \$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, \$ IERR ) * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in WORK(IVT) * (Workspace: need M*M+BDSPAC) * CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, \$ WORK( IVT ), LDWKVT, DUM, IDUM, \$ WORK( NWORK ), IWORK, INFO ) * * Overwrite U by left singular vectors of L and WORK(IVT) * by right singular vectors of L * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA, \$ WORK( ITAUP ), WORK( IVT ), LDWKVT, \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Multiply right singular vectors of L in WORK(IVT) by * Q in VT, storing result in A * (Workspace: need M*M) * CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT, \$ VT, LDVT, ZERO, A, LDA ) * * Copy right singular vectors of A from A to VT * CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT ) * END IF * ELSE * * N .LT. MNTHR * * Path 5t (N greater than M, but not much larger) * Reduce to bidiagonal form without LQ decomposition * IE = 1 ITAUQ = IE + M ITAUP = ITAUQ + M NWORK = ITAUP + M * * Bidiagonalize A * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) * CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), \$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, \$ IERR ) IF( WNTQN ) THEN * * Perform bidiagonal SVD, only computing singular values * (Workspace: need BDSPAN) * CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1, \$ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) ELSE IF( WNTQO ) THEN LDWKVT = M IVT = NWORK IF( LWORK.GE.M*N+3*M*M+9*M ) THEN * * WORK( IVT ) is M by N * CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ), \$ LDWKVT ) NWORK = IVT + LDWKVT*N ELSE * * WORK( IVT ) is M by M * NWORK = IVT + LDWKVT*M IL = NWORK * * WORK(IL) is M by CHUNK * CHUNK = ( LWORK-M*M-3*M ) / M END IF * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in WORK(IVT) * (Workspace: need M*M+BDSPAC) * CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, \$ WORK( IVT ), LDWKVT, DUM, IDUM, \$ WORK( NWORK ), IWORK, INFO ) * * Overwrite U by left singular vectors of A * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) * IF( LWORK.GE.M*N+3*M*M+9*M ) THEN * * Overwrite WORK(IVT) by left singular vectors of A * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA, \$ WORK( ITAUP ), WORK( IVT ), LDWKVT, \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Copy right singular vectors of A from WORK(IVT) to A * CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) ELSE * * Generate P**T in A * (Workspace: need M*M+2*M, prefer M*M+M+M*NB) * CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), \$ WORK( NWORK ), LWORK-NWORK+1, IERR ) * * Multiply Q in A by right singular vectors of * bidiagonal matrix in WORK(IVT), storing result in * WORK(IL) and copying to A * (Workspace: need 2*M*M, prefer M*M+M*N) * DO 40 I = 1, N, CHUNK BLK = MIN( N-I+1, CHUNK ) CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), \$ LDWKVT, A( 1, I ), LDA, ZERO, \$ WORK( IL ), M ) CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ), \$ LDA ) 40 CONTINUE END IF ELSE IF( WNTQS ) THEN * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in VT * (Workspace: need BDSPAC) * CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT ) CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT, \$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Overwrite U by left singular vectors of A and VT * by right singular vectors of A * (Workspace: need 3*M, prefer 2*M+M*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) ELSE * * Perform bidiagonal SVD, computing left singular vectors * of bidiagonal matrix in U and computing right singular * vectors of bidiagonal matrix in VT * (Workspace: need BDSPAC) * CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT ) CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT, \$ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, \$ INFO ) * * Set the right corner of VT to identity matrix * CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ), \$ LDVT ) * * Overwrite U by left singular vectors of A and VT * by right singular vectors of A * (Workspace: need 2*M+N, prefer 2*M+N*NB) * CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, \$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA, \$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), \$ LWORK-NWORK+1, IERR ) END IF * END IF * END IF * * Undo scaling if necessary * IF( ISCL.EQ.1 ) THEN IF( ANRM.GT.BIGNUM ) \$ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, \$ IERR ) IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM ) \$ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, WORK( 2 ), \$ MINMN, IERR ) IF( ANRM.LT.SMLNUM ) \$ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, \$ IERR ) IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM ) \$ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, WORK( 2 ), \$ MINMN, IERR ) END IF * * Return optimal workspace in WORK(1) * WORK( 1 ) = DBLE( MAXWRK ) * RETURN * * End of DGESDD * END