version 1.13, 2012/12/14 14:22:29
|
version 1.17, 2017/06/17 10:53:49
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DGESVD + dependencies |
*> Download DGESVD + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, |
* SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, |
* WORK, LWORK, INFO ) |
* WORK, LWORK, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER JOBU, JOBVT |
* CHARACTER JOBU, JOBVT |
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N |
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N |
Line 29
|
Line 29
|
* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), |
* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), |
* $ VT( LDVT, * ), WORK( * ) |
* $ VT( LDVT, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 173
|
Line 173
|
*> LWORK is INTEGER |
*> LWORK is INTEGER |
*> The dimension of the array WORK. |
*> The dimension of the array WORK. |
*> LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code): |
*> LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code): |
*> - PATH 1 (M much larger than N, JOBU='N') |
*> - PATH 1 (M much larger than N, JOBU='N') |
*> - PATH 1t (N much larger than M, JOBVT='N') |
*> - PATH 1t (N much larger than M, JOBVT='N') |
*> LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) for the other paths |
*> LWORK >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths |
*> For good performance, LWORK should generally be larger. |
*> For good performance, LWORK should generally be larger. |
*> |
*> |
*> If LWORK = -1, then a workspace query is assumed; the routine |
*> If LWORK = -1, then a workspace query is assumed; the routine |
Line 198
|
Line 198
|
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date April 2012 |
*> \date April 2012 |
* |
* |
Line 211
|
Line 211
|
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, |
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, |
$ VT, LDVT, WORK, LWORK, INFO ) |
$ VT, LDVT, WORK, LWORK, INFO ) |
* |
* |
* -- LAPACK driver routine (version 3.4.1) -- |
* -- LAPACK driver routine (version 3.7.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* April 2012 |
* April 2012 |
Line 314
|
Line 314
|
BDSPAC = 5*N |
BDSPAC = 5*N |
* Compute space needed for DGEQRF |
* Compute space needed for DGEQRF |
CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
LWORK_DGEQRF=DUM(1) |
LWORK_DGEQRF = INT( DUM(1) ) |
* Compute space needed for DORGQR |
* Compute space needed for DORGQR |
CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
LWORK_DORGQR_N=DUM(1) |
LWORK_DORGQR_N = INT( DUM(1) ) |
CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
LWORK_DORGQR_M=DUM(1) |
LWORK_DORGQR_M = INT( DUM(1) ) |
* Compute space needed for DGEBRD |
* Compute space needed for DGEBRD |
CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), |
CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), |
$ DUM(1), DUM(1), -1, IERR ) |
$ DUM(1), DUM(1), -1, IERR ) |
LWORK_DGEBRD=DUM(1) |
LWORK_DGEBRD = INT( DUM(1) ) |
* Compute space needed for DORGBR P |
* Compute space needed for DORGBR P |
CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1), |
CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_P=DUM(1) |
LWORK_DORGBR_P = INT( DUM(1) ) |
* Compute space needed for DORGBR Q |
* Compute space needed for DORGBR Q |
CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1), |
CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_Q=DUM(1) |
LWORK_DORGBR_Q = INT( DUM(1) ) |
* |
* |
IF( M.GE.MNTHR ) THEN |
IF( M.GE.MNTHR ) THEN |
IF( WNTUN ) THEN |
IF( WNTUN ) THEN |
Line 339
|
Line 339
|
* Path 1 (M much larger than N, JOBU='N') |
* Path 1 (M much larger than N, JOBU='N') |
* |
* |
MAXWRK = N + LWORK_DGEQRF |
MAXWRK = N + LWORK_DGEQRF |
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DGEBRD ) |
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD ) |
IF( WNTVO .OR. WNTVAS ) |
IF( WNTVO .OR. WNTVAS ) |
$ MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P ) |
$ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 4*N, BDSPAC ) |
MINWRK = MAX( 4*N, BDSPAC ) |
ELSE IF( WNTUO .AND. WNTVN ) THEN |
ELSE IF( WNTUO .AND. WNTVN ) THEN |
Line 349
|
Line 349
|
* Path 2 (M much larger than N, JOBU='O', JOBVT='N') |
* Path 2 (M much larger than N, JOBU='O', JOBVT='N') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) |
MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N ) |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUO .AND. WNTVAS ) THEN |
ELSE IF( WNTUO .AND. WNTVAS ) THEN |
* |
* |
* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or |
* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or |
* 'A') |
* 'A') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) |
MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N ) |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUS .AND. WNTVN ) THEN |
ELSE IF( WNTUS .AND. WNTVN ) THEN |
* |
* |
* Path 4 (M much larger than N, JOBU='S', JOBVT='N') |
* Path 4 (M much larger than N, JOBU='S', JOBVT='N') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUS .AND. WNTVO ) THEN |
ELSE IF( WNTUS .AND. WNTVO ) THEN |
* |
* |
* Path 5 (M much larger than N, JOBU='S', JOBVT='O') |
* Path 5 (M much larger than N, JOBU='S', JOBVT='O') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*N*N + WRKBL |
MAXWRK = 2*N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUS .AND. WNTVAS ) THEN |
ELSE IF( WNTUS .AND. WNTVAS ) THEN |
* |
* |
* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or |
* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or |
* 'A') |
* 'A') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUA .AND. WNTVN ) THEN |
ELSE IF( WNTUA .AND. WNTVN ) THEN |
* |
* |
* Path 7 (M much larger than N, JOBU='A', JOBVT='N') |
* Path 7 (M much larger than N, JOBU='A', JOBVT='N') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUA .AND. WNTVO ) THEN |
ELSE IF( WNTUA .AND. WNTVO ) THEN |
* |
* |
* Path 8 (M much larger than N, JOBU='A', JOBVT='O') |
* Path 8 (M much larger than N, JOBU='A', JOBVT='O') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*N*N + WRKBL |
MAXWRK = 2*N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUA .AND. WNTVAS ) THEN |
ELSE IF( WNTUA .AND. WNTVAS ) THEN |
* |
* |
* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or |
* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or |
* 'A') |
* 'A') |
* |
* |
WRKBL = N + LWORK_DGEQRF |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M ) |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
END IF |
END IF |
ELSE |
ELSE |
* |
* |
Line 447
|
Line 447
|
* |
* |
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), |
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), |
$ DUM(1), DUM(1), -1, IERR ) |
$ DUM(1), DUM(1), -1, IERR ) |
LWORK_DGEBRD=DUM(1) |
LWORK_DGEBRD = INT( DUM(1) ) |
MAXWRK = 3*N + LWORK_DGEBRD |
MAXWRK = 3*N + LWORK_DGEBRD |
IF( WNTUS .OR. WNTUO ) THEN |
IF( WNTUS .OR. WNTUO ) THEN |
CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1), |
CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_Q=DUM(1) |
LWORK_DORGBR_Q = INT( DUM(1) ) |
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q ) |
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q ) |
END IF |
END IF |
IF( WNTUA ) THEN |
IF( WNTUA ) THEN |
CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1), |
CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_Q=DUM(1) |
LWORK_DORGBR_Q = INT( DUM(1) ) |
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q ) |
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q ) |
END IF |
END IF |
IF( .NOT.WNTVN ) THEN |
IF( .NOT.WNTVN ) THEN |
MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P ) |
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P ) |
END IF |
END IF |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
END IF |
END IF |
ELSE IF( MINMN.GT.0 ) THEN |
ELSE IF( MINMN.GT.0 ) THEN |
* |
* |
Line 475
|
Line 475
|
BDSPAC = 5*M |
BDSPAC = 5*M |
* Compute space needed for DGELQF |
* Compute space needed for DGELQF |
CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
LWORK_DGELQF=DUM(1) |
LWORK_DGELQF = INT( DUM(1) ) |
* Compute space needed for DORGLQ |
* Compute space needed for DORGLQ |
CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR ) |
CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR ) |
LWORK_DORGLQ_N=DUM(1) |
LWORK_DORGLQ_N = INT( DUM(1) ) |
CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) |
CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) |
LWORK_DORGLQ_M=DUM(1) |
LWORK_DORGLQ_M = INT( DUM(1) ) |
* Compute space needed for DGEBRD |
* Compute space needed for DGEBRD |
CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), |
CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), |
$ DUM(1), DUM(1), -1, IERR ) |
$ DUM(1), DUM(1), -1, IERR ) |
LWORK_DGEBRD=DUM(1) |
LWORK_DGEBRD = INT( DUM(1) ) |
* Compute space needed for DORGBR P |
* Compute space needed for DORGBR P |
CALL DORGBR( 'P', M, M, M, A, N, DUM(1), |
CALL DORGBR( 'P', M, M, M, A, N, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_P=DUM(1) |
LWORK_DORGBR_P = INT( DUM(1) ) |
* Compute space needed for DORGBR Q |
* Compute space needed for DORGBR Q |
CALL DORGBR( 'Q', M, M, M, A, N, DUM(1), |
CALL DORGBR( 'Q', M, M, M, A, N, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_Q=DUM(1) |
LWORK_DORGBR_Q = INT( DUM(1) ) |
IF( N.GE.MNTHR ) THEN |
IF( N.GE.MNTHR ) THEN |
IF( WNTVN ) THEN |
IF( WNTVN ) THEN |
* |
* |
* Path 1t(N much larger than M, JOBVT='N') |
* Path 1t(N much larger than M, JOBVT='N') |
* |
* |
MAXWRK = M + LWORK_DGELQF |
MAXWRK = M + LWORK_DGELQF |
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DGEBRD ) |
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DGEBRD ) |
IF( WNTUO .OR. WNTUAS ) |
IF( WNTUO .OR. WNTUAS ) |
$ MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q ) |
$ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 4*M, BDSPAC ) |
MINWRK = MAX( 4*M, BDSPAC ) |
ELSE IF( WNTVO .AND. WNTUN ) THEN |
ELSE IF( WNTVO .AND. WNTUN ) THEN |
Line 509
|
Line 509
|
* Path 2t(N much larger than M, JOBU='N', JOBVT='O') |
* Path 2t(N much larger than M, JOBU='N', JOBVT='O') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) |
MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M ) |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVO .AND. WNTUAS ) THEN |
ELSE IF( WNTVO .AND. WNTUAS ) THEN |
* |
* |
* Path 3t(N much larger than M, JOBU='S' or 'A', |
* Path 3t(N much larger than M, JOBU='S' or 'A', |
* JOBVT='O') |
* JOBVT='O') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) |
MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M ) |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVS .AND. WNTUN ) THEN |
ELSE IF( WNTVS .AND. WNTUN ) THEN |
* |
* |
* Path 4t(N much larger than M, JOBU='N', JOBVT='S') |
* Path 4t(N much larger than M, JOBU='N', JOBVT='S') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVS .AND. WNTUO ) THEN |
ELSE IF( WNTVS .AND. WNTUO ) THEN |
* |
* |
* Path 5t(N much larger than M, JOBU='O', JOBVT='S') |
* Path 5t(N much larger than M, JOBU='O', JOBVT='S') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*M*M + WRKBL |
MAXWRK = 2*M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVS .AND. WNTUAS ) THEN |
ELSE IF( WNTVS .AND. WNTUAS ) THEN |
* |
* |
* Path 6t(N much larger than M, JOBU='S' or 'A', |
* Path 6t(N much larger than M, JOBU='S' or 'A', |
* JOBVT='S') |
* JOBVT='S') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVA .AND. WNTUN ) THEN |
ELSE IF( WNTVA .AND. WNTUN ) THEN |
* |
* |
* Path 7t(N much larger than M, JOBU='N', JOBVT='A') |
* Path 7t(N much larger than M, JOBU='N', JOBVT='A') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVA .AND. WNTUO ) THEN |
ELSE IF( WNTVA .AND. WNTUO ) THEN |
* |
* |
* Path 8t(N much larger than M, JOBU='O', JOBVT='A') |
* Path 8t(N much larger than M, JOBU='O', JOBVT='A') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*M*M + WRKBL |
MAXWRK = 2*M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVA .AND. WNTUAS ) THEN |
ELSE IF( WNTVA .AND. WNTUAS ) THEN |
* |
* |
* Path 9t(N much larger than M, JOBU='S' or 'A', |
* Path 9t(N much larger than M, JOBU='S' or 'A', |
* JOBVT='A') |
* JOBVT='A') |
* |
* |
WRKBL = M + LWORK_DGELQF |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N ) |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
END IF |
END IF |
ELSE |
ELSE |
* |
* |
Line 607
|
Line 607
|
* |
* |
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), |
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), |
$ DUM(1), DUM(1), -1, IERR ) |
$ DUM(1), DUM(1), -1, IERR ) |
LWORK_DGEBRD=DUM(1) |
LWORK_DGEBRD = INT( DUM(1) ) |
MAXWRK = 3*M + LWORK_DGEBRD |
MAXWRK = 3*M + LWORK_DGEBRD |
IF( WNTVS .OR. WNTVO ) THEN |
IF( WNTVS .OR. WNTVO ) THEN |
* Compute space needed for DORGBR P |
* Compute space needed for DORGBR P |
CALL DORGBR( 'P', M, N, M, A, N, DUM(1), |
CALL DORGBR( 'P', M, N, M, A, N, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_P=DUM(1) |
LWORK_DORGBR_P = INT( DUM(1) ) |
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P ) |
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P ) |
END IF |
END IF |
IF( WNTVA ) THEN |
IF( WNTVA ) THEN |
CALL DORGBR( 'P', N, N, M, A, N, DUM(1), |
CALL DORGBR( 'P', N, N, M, A, N, DUM(1), |
$ DUM(1), -1, IERR ) |
$ DUM(1), -1, IERR ) |
LWORK_DORGBR_P=DUM(1) |
LWORK_DORGBR_P = INT( DUM(1) ) |
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P ) |
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P ) |
END IF |
END IF |
IF( .NOT.WNTUN ) THEN |
IF( .NOT.WNTUN ) THEN |
MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q ) |
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q ) |
END IF |
END IF |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
END IF |
END IF |
END IF |
END IF |
MAXWRK = MAX( MAXWRK, MINWRK ) |
MAXWRK = MAX( MAXWRK, MINWRK ) |
Line 685
|
Line 685
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Zero out below R |
* Zero out below R |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) |
IF( N .GT. 1 ) THEN |
|
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
|
$ LDA ) |
|
END IF |
IE = 1 |
IE = 1 |
ITAUQ = IE + N |
ITAUQ = IE + N |
ITAUP = ITAUQ + N |
ITAUP = ITAUQ + N |
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 708
|
Line 711
|
IF( WNTVO .OR. WNTVAS ) THEN |
IF( WNTVO .OR. WNTVAS ) THEN |
* |
* |
* If right singular vectors desired, generate P'. |
* If right singular vectors desired, generate P'. |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 739
|
Line 742
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN |
* |
* |
* WORK(IU) is LDA by N, WORK(IR) is LDA by N |
* WORK(IU) is LDA by N, WORK(IR) is LDA by N |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN |
* |
* |
* WORK(IU) is LDA by N, WORK(IR) is N by N |
* WORK(IU) is LDA by N, WORK(IR) is N by N |
* |
* |
Line 762
|
Line 765
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 774
|
Line 777
|
$ LDWRKR ) |
$ LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 784
|
Line 787
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing R |
* Generate left vectors bidiagonalizing R |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 800
|
Line 803
|
* |
* |
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) |
* singular vectors of R in WORK(IR) |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1, |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, |
Line 809
|
Line 812
|
* |
* |
* Multiply Q in A by left singular vectors of R in |
* Multiply Q in A by left singular vectors of R in |
* WORK(IR), storing result in WORK(IU) and copying to A |
* WORK(IR), storing result in WORK(IU) and copying to A |
* (Workspace: need N*N+2*N, prefer N*N+M*N+N) |
* (Workspace: need N*N + 2*N, prefer N*N + M*N + N) |
* |
* |
DO 10 I = 1, M, LDWRKU |
DO 10 I = 1, M, LDWRKU |
CHUNK = MIN( M-I+1, LDWRKU ) |
CHUNK = MIN( M-I+1, LDWRKU ) |
Line 830
|
Line 833
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) |
* (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing A |
* Generate left vectors bidiagonalizing A |
* (Workspace: need 4*N, prefer 3*N+N*NB) |
* (Workspace: need 4*N, prefer 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 863
|
Line 866
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is LDA by N |
* WORK(IU) is LDA by N and WORK(IR) is LDA by N |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* |
* |
Line 886
|
Line 889
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 899
|
Line 902
|
$ VT( 2, 1 ), LDVT ) |
$ VT( 2, 1 ), LDVT ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 909
|
Line 912
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT, copying result to WORK(IR) |
* Bidiagonalize R in VT, copying result to WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 917
|
Line 920
|
CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) |
CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) |
* |
* |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing R in VT |
* Generate right vectors bidiagonalizing R in VT |
* (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB) |
* (Workspace: need N*N + 4*N-1, prefer N*N + 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 933
|
Line 936
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) and computing right |
* singular vectors of R in WORK(IR) and computing right |
* singular vectors of R in VT |
* singular vectors of R in VT |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT, |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT, |
$ WORK( IR ), LDWRKR, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, |
Line 942
|
Line 945
|
* |
* |
* Multiply Q in A by left singular vectors of R in |
* Multiply Q in A by left singular vectors of R in |
* WORK(IR), storing result in WORK(IU) and copying to A |
* WORK(IR), storing result in WORK(IU) and copying to A |
* (Workspace: need N*N+2*N, prefer N*N+M*N+N) |
* (Workspace: need N*N + 2*N, prefer N*N + M*N + N) |
* |
* |
DO 20 I = 1, M, LDWRKU |
DO 20 I = 1, M, LDWRKU |
CHUNK = MIN( M-I+1, LDWRKU ) |
CHUNK = MIN( M-I+1, LDWRKU ) |
Line 961
|
Line 964
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 974
|
Line 977
|
$ VT( 2, 1 ), LDVT ) |
$ VT( 2, 1 ), LDVT ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 984
|
Line 987
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT |
* Bidiagonalize R in VT |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply Q in A by left vectors bidiagonalizing R |
* Multiply Q in A by left vectors bidiagonalizing R |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
$ WORK( ITAUQ ), A, LDA, WORK( IWORK ), |
$ WORK( ITAUQ ), A, LDA, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing R in VT |
* Generate right vectors bidiagonalizing R in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1042
|
Line 1045
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1055
|
Line 1058
|
$ WORK( IR+1 ), LDWRKR ) |
$ WORK( IR+1 ), LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1065
|
Line 1068
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1073
|
Line 1076
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 1082
|
Line 1085
|
* |
* |
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) |
* singular vectors of R in WORK(IR) |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
Line 1103
|
Line 1106
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1121
|
Line 1124
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply Q in U by left vectors bidiagonalizing R |
* Multiply Q in U by left vectors bidiagonalizing R |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
Line 1167
|
Line 1172
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*N |
IR = IU + LDWRKU*N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* |
* |
Line 1186
|
Line 1191
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) |
* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1199
|
Line 1204
|
$ WORK( IU+1 ), LDWRKU ) |
$ WORK( IU+1 ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) |
* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1210
|
Line 1215
|
* |
* |
* Bidiagonalize R in WORK(IU), copying result to |
* Bidiagonalize R in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*N*N+4*N, |
* (Workspace: need 2*N*N + 4*N, |
* prefer 2*N*N+3*N+2*N*NB) |
* prefer 2*N*N+3*N+2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
Line 1221
|
Line 1226
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB) |
* (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IR) |
* Generate right bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*N*N+4*N-1, |
* (Workspace: need 2*N*N + 4*N-1, |
* prefer 2*N*N+3*N+(N-1)*NB) |
* prefer 2*N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
Line 1239
|
Line 1244
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in WORK(IR) |
* right singular vectors of R in WORK(IR) |
* (Workspace: need 2*N*N+BDSPAC) |
* (Workspace: need 2*N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
Line 1266
|
Line 1271
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1284
|
Line 1289
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply Q in U by left vectors bidiagonalizing R |
* Multiply Q in U by left vectors bidiagonalizing R |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing R in A |
* Generate right vectors bidiagonalizing R in A |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1346
|
Line 1353
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1359
|
Line 1366
|
$ WORK( IU+1 ), LDWRKU ) |
$ WORK( IU+1 ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1369
|
Line 1376
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IU), copying result to VT |
* Bidiagonalize R in WORK(IU), copying result to VT |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1379
|
Line 1386
|
$ LDVT ) |
$ LDVT ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need N*N+4*N-1, |
* (Workspace: need N*N + 4*N-1, |
* prefer N*N+3*N+(N-1)*NB) |
* prefer N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
Line 1396
|
Line 1403
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in VT |
* right singular vectors of R in VT |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
Line 1417
|
Line 1424
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1441
|
Line 1448
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT |
* Bidiagonalize R in VT |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1449
|
Line 1456
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in VT |
* in VT |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1503
|
Line 1510
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1517
|
Line 1524
|
$ WORK( IR+1 ), LDWRKR ) |
$ WORK( IR+1 ), LDWRKR ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N*N+N+M, prefer N*N+N+M*NB) |
* (Workspace: need N*N + N + M, prefer N*N + N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1527
|
Line 1534
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1535
|
Line 1542
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IR) |
* Generate left bidiagonalizing vectors in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 1544
|
Line 1551
|
* |
* |
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) |
* singular vectors of R in WORK(IR) |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
Line 1569
|
Line 1576
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N+M, prefer N+M*NB) |
* (Workspace: need N + M, prefer N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1587
|
Line 1594
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1599
|
Line 1608
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in A |
* in A |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
Line 1634
|
Line 1643
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*N |
IR = IU + LDWRKU*N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* |
* |
Line 1653
|
Line 1662
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) |
* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) |
* (Workspace: need 2*N*N + N + M, prefer 2*N*N + N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1678
|
Line 1687
|
* |
* |
* Bidiagonalize R in WORK(IU), copying result to |
* Bidiagonalize R in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*N*N+4*N, |
* (Workspace: need 2*N*N + 4*N, |
* prefer 2*N*N+3*N+2*N*NB) |
* prefer 2*N*N+3*N+2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
Line 1689
|
Line 1698
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB) |
* (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IR) |
* Generate right bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*N*N+4*N-1, |
* (Workspace: need 2*N*N + 4*N-1, |
* prefer 2*N*N+3*N+(N-1)*NB) |
* prefer 2*N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
Line 1707
|
Line 1716
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in WORK(IR) |
* right singular vectors of R in WORK(IR) |
* (Workspace: need 2*N*N+BDSPAC) |
* (Workspace: need 2*N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
Line 1737
|
Line 1746
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N+M, prefer N+M*NB) |
* (Workspace: need N + M, prefer N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1755
|
Line 1764
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1767
|
Line 1778
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in A |
* in A |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in A |
* Generate right bidiagonalizing vectors in A |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1818
|
Line 1829
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N*N+N+M, prefer N*N+N+M*NB) |
* (Workspace: need N*N + N + M, prefer N*N + N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1842
|
Line 1853
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IU), copying result to VT |
* Bidiagonalize R in WORK(IU), copying result to VT |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1852
|
Line 1863
|
$ LDVT ) |
$ LDVT ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need N*N+4*N-1, |
* (Workspace: need N*N + 4*N-1, |
* prefer N*N+3*N+(N-1)*NB) |
* prefer N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
Line 1869
|
Line 1880
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in VT |
* right singular vectors of R in VT |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
Line 1894
|
Line 1905
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N+M, prefer N+M*NB) |
* (Workspace: need N + M, prefer N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1918
|
Line 1929
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT |
* Bidiagonalize R in VT |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1926
|
Line 1937
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in VT |
* in VT |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1967
|
Line 1978
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) |
* (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 1976
|
Line 1987
|
* |
* |
* If left singular vectors desired in U, copy result to U |
* If left singular vectors desired in U, copy result to U |
* and generate left bidiagonalizing vectors in U |
* and generate left bidiagonalizing vectors in U |
* (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB) |
* (Workspace: need 3*N + NCU, prefer 3*N + NCU*NB) |
* |
* |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
IF( WNTUS ) |
IF( WNTUS ) |
Line 1990
|
Line 2001
|
* |
* |
* If right singular vectors desired in VT, copy result to |
* If right singular vectors desired in VT, copy result to |
* VT and generate right bidiagonalizing vectors in VT |
* VT and generate right bidiagonalizing vectors in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
Line 2000
|
Line 2011
|
* |
* |
* If left singular vectors desired in A, generate left |
* If left singular vectors desired in A, generate left |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*N, prefer 3*N+N*NB) |
* (Workspace: need 4*N, prefer 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2009
|
Line 2020
|
* |
* |
* If right singular vectors desired in A, generate right |
* If right singular vectors desired in A, generate right |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2071
|
Line 2082
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
Line 2085
|
Line 2096
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 2093
|
Line 2104
|
IF( WNTUO .OR. WNTUAS ) THEN |
IF( WNTUO .OR. WNTUAS ) THEN |
* |
* |
* If left singular vectors desired, generate Q |
* If left singular vectors desired, generate Q |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2126
|
Line 2137
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
CHUNK = N |
CHUNK = N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* |
* |
Line 2152
|
Line 2163
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2164
|
Line 2175
|
$ WORK( IR+LDWRKR ), LDWRKR ) |
$ WORK( IR+LDWRKR ), LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2174
|
Line 2185
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IR) |
* Bidiagonalize L in WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ), |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing L |
* Generate right vectors bidiagonalizing L |
* (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
Line 2190
|
Line 2201
|
* |
* |
* Perform bidiagonal QR iteration, computing right |
* Perform bidiagonal QR iteration, computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
Line 2199
|
Line 2210
|
* |
* |
* Multiply right singular vectors of L in WORK(IR) by Q |
* Multiply right singular vectors of L in WORK(IR) by Q |
* in A, storing result in WORK(IU) and copying to A |
* in A, storing result in WORK(IU) and copying to A |
* (Workspace: need M*M+2*M, prefer M*M+M*N+M) |
* (Workspace: need M*M + 2*M, prefer M*M + M*N + M) |
* |
* |
DO 30 I = 1, N, CHUNK |
DO 30 I = 1, N, CHUNK |
BLK = MIN( N-I+1, CHUNK ) |
BLK = MIN( N-I+1, CHUNK ) |
Line 2220
|
Line 2231
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) |
* (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing A |
* Generate right vectors bidiagonalizing A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2253
|
Line 2264
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
CHUNK = N |
CHUNK = N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* |
* |
Line 2279
|
Line 2290
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2291
|
Line 2302
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2301
|
Line 2312
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U, copying result to WORK(IR) |
* Bidiagonalize L in U, copying result to WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 2309
|
Line 2320
|
CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) |
CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) |
* |
* |
* Generate right vectors bidiagonalizing L in WORK(IR) |
* Generate right vectors bidiagonalizing L in WORK(IR) |
* (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing L in U |
* Generate left vectors bidiagonalizing L in U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2325
|
Line 2336
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in U, and computing right |
* singular vectors of L in U, and computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, U, LDU, DUM, 1, |
$ WORK( IR ), LDWRKR, U, LDU, DUM, 1, |
Line 2334
|
Line 2345
|
* |
* |
* Multiply right singular vectors of L in WORK(IR) by Q |
* Multiply right singular vectors of L in WORK(IR) by Q |
* in A, storing result in WORK(IU) and copying to A |
* in A, storing result in WORK(IU) and copying to A |
* (Workspace: need M*M+2*M, prefer M*M+M*N+M)) |
* (Workspace: need M*M + 2*M, prefer M*M + M*N + M)) |
* |
* |
DO 40 I = 1, N, CHUNK |
DO 40 I = 1, N, CHUNK |
BLK = MIN( N-I+1, CHUNK ) |
BLK = MIN( N-I+1, CHUNK ) |
Line 2353
|
Line 2364
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2365
|
Line 2376
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2375
|
Line 2386
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U |
* Bidiagonalize L in U |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply right vectors bidiagonalizing L by Q in A |
* Multiply right vectors bidiagonalizing L by Q in A |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
$ WORK( ITAUP ), A, LDA, WORK( IWORK ), |
$ WORK( ITAUP ), A, LDA, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing L in U |
* Generate left vectors bidiagonalizing L in U |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2433
|
Line 2444
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2446
|
Line 2457
|
$ WORK( IR+LDWRKR ), LDWRKR ) |
$ WORK( IR+LDWRKR ), LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2456
|
Line 2467
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IR) |
* Bidiagonalize L in WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 2465
|
Line 2476
|
* |
* |
* Generate right vectors bidiagonalizing L in |
* Generate right vectors bidiagonalizing L in |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
Line 2474
|
Line 2485
|
* |
* |
* Perform bidiagonal QR iteration, computing right |
* Perform bidiagonal QR iteration, computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
Line 2495
|
Line 2506
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2505
|
Line 2516
|
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2520
|
Line 2531
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply right vectors bidiagonalizing L by Q in VT |
* Multiply right vectors bidiagonalizing L by Q in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
Line 2562
|
Line 2573
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*M |
IR = IU + LDWRKU*M |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN |
* |
* |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* |
* |
Line 2581
|
Line 2592
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) |
* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2594
|
Line 2605
|
$ WORK( IU+LDWRKU ), LDWRKU ) |
$ WORK( IU+LDWRKU ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) |
* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2605
|
Line 2616
|
* |
* |
* Bidiagonalize L in WORK(IU), copying result to |
* Bidiagonalize L in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*M*M+4*M, |
* (Workspace: need 2*M*M + 4*M, |
* prefer 2*M*M+3*M+2*M*NB) |
* prefer 2*M*M+3*M+2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
Line 2616
|
Line 2627
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*M*M+4*M-1, |
* (Workspace: need 2*M*M + 4*M-1, |
* prefer 2*M*M+3*M+(M-1)*NB) |
* prefer 2*M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
Line 2624
|
Line 2635
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IR) |
* Generate left bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB) |
* (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 2634
|
Line 2645
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in WORK(IR) and computing |
* singular vectors of L in WORK(IR) and computing |
* right singular vectors of L in WORK(IU) |
* right singular vectors of L in WORK(IU) |
* (Workspace: need 2*M*M+BDSPAC) |
* (Workspace: need 2*M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
Line 2661
|
Line 2672
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2683
|
Line 2694
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply right vectors bidiagonalizing L by Q in VT |
* Multiply right vectors bidiagonalizing L by Q in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors of L in A |
* Generate left bidiagonalizing vectors of L in A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2741
|
Line 2752
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2754
|
Line 2765
|
$ WORK( IU+LDWRKU ), LDWRKU ) |
$ WORK( IU+LDWRKU ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2764
|
Line 2775
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IU), copying result to U |
* Bidiagonalize L in WORK(IU), copying result to U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 2774
|
Line 2785
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need M*M+4*M-1, |
* (Workspace: need M*M + 4*M-1, |
* prefer M*M+3*M+(M-1)*NB) |
* prefer M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
Line 2782
|
Line 2793
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2791
|
Line 2802
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in U and computing right |
* singular vectors of L in U and computing right |
* singular vectors of L in WORK(IU) |
* singular vectors of L in WORK(IU) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
Line 2812
|
Line 2823
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2835
|
Line 2846
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U |
* Bidiagonalize L in U |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 2843
|
Line 2854
|
* |
* |
* Multiply right bidiagonalizing vectors in U by Q |
* Multiply right bidiagonalizing vectors in U by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2877
|
Line 2888
|
* N right singular vectors to be computed in VT and |
* N right singular vectors to be computed in VT and |
* no left singular vectors to be computed |
* no left singular vectors to be computed |
* |
* |
IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN |
IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN |
* |
* |
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
Line 2897
|
Line 2908
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2911
|
Line 2922
|
$ WORK( IR+LDWRKR ), LDWRKR ) |
$ WORK( IR+LDWRKR ), LDWRKR ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M*M+M+N, prefer M*M+M+N*NB) |
* (Workspace: need M*M + M + N, prefer M*M + M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2921
|
Line 2932
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IR) |
* Bidiagonalize L in WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 2929
|
Line 2940
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IR) |
* Generate right bidiagonalizing vectors in WORK(IR) |
* (Workspace: need M*M+4*M-1, |
* (Workspace: need M*M + 4*M-1, |
* prefer M*M+3*M+(M-1)*NB) |
* prefer M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
Line 2939
|
Line 2950
|
* |
* |
* Perform bidiagonal QR iteration, computing right |
* Perform bidiagonal QR iteration, computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
Line 2964
|
Line 2975
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M+N, prefer M+N*NB) |
* (Workspace: need M + N, prefer M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2986
|
Line 2997
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 2994
|
Line 3005
|
* |
* |
* Multiply right bidiagonalizing vectors in A by Q |
* Multiply right bidiagonalizing vectors in A by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
Line 3017
|
Line 3028
|
* N right singular vectors to be computed in VT and |
* N right singular vectors to be computed in VT and |
* M left singular vectors to be overwritten on A |
* M left singular vectors to be overwritten on A |
* |
* |
IF( LWORK.GE.2*M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN |
IF( LWORK.GE.2*M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN |
* |
* |
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
Line 3029
|
Line 3040
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*M |
IR = IU + LDWRKU*M |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN |
* |
* |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* |
* |
Line 3048
|
Line 3059
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) |
* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) |
* (Workspace: need 2*M*M + M + N, prefer 2*M*M + M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3073
|
Line 3084
|
* |
* |
* Bidiagonalize L in WORK(IU), copying result to |
* Bidiagonalize L in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*M*M+4*M, |
* (Workspace: need 2*M*M + 4*M, |
* prefer 2*M*M+3*M+2*M*NB) |
* prefer 2*M*M+3*M+2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
Line 3084
|
Line 3095
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*M*M+4*M-1, |
* (Workspace: need 2*M*M + 4*M-1, |
* prefer 2*M*M+3*M+(M-1)*NB) |
* prefer 2*M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
Line 3092
|
Line 3103
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IR) |
* Generate left bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB) |
* (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 3102
|
Line 3113
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in WORK(IR) and computing |
* singular vectors of L in WORK(IR) and computing |
* right singular vectors of L in WORK(IU) |
* right singular vectors of L in WORK(IU) |
* (Workspace: need 2*M*M+BDSPAC) |
* (Workspace: need 2*M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
Line 3132
|
Line 3143
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M+N, prefer M+N*NB) |
* (Workspace: need M + N, prefer M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3154
|
Line 3165
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 3162
|
Line 3173
|
* |
* |
* Multiply right bidiagonalizing vectors in A by Q |
* Multiply right bidiagonalizing vectors in A by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in A |
* Generate left bidiagonalizing vectors in A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3193
|
Line 3204
|
* N right singular vectors to be computed in VT and |
* N right singular vectors to be computed in VT and |
* M left singular vectors to be computed in U |
* M left singular vectors to be computed in U |
* |
* |
IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN |
IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN |
* |
* |
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
Line 3213
|
Line 3224
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M*M+M+N, prefer M*M+M+N*NB) |
* (Workspace: need M*M + M + N, prefer M*M + M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3237
|
Line 3248
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IU), copying result to U |
* Bidiagonalize L in WORK(IU), copying result to U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 3247
|
Line 3258
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3263
|
Line 3274
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in U and computing right |
* singular vectors of L in U and computing right |
* singular vectors of L in WORK(IU) |
* singular vectors of L in WORK(IU) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
Line 3288
|
Line 3299
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M+N, prefer M+N*NB) |
* (Workspace: need M + N, prefer M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3311
|
Line 3322
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U |
* Bidiagonalize L in U |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 3319
|
Line 3330
|
* |
* |
* Multiply right bidiagonalizing vectors in U by Q |
* Multiply right bidiagonalizing vectors in U by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3360
|
Line 3371
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) |
* (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 3369
|
Line 3380
|
* |
* |
* If left singular vectors desired in U, copy result to U |
* If left singular vectors desired in U, copy result to U |
* and generate left bidiagonalizing vectors in U |
* and generate left bidiagonalizing vectors in U |
* (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB) |
* (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB) |
* |
* |
CALL DLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
Line 3379
|
Line 3390
|
* |
* |
* If right singular vectors desired in VT, copy result to |
* If right singular vectors desired in VT, copy result to |
* VT and generate right bidiagonalizing vectors in VT |
* VT and generate right bidiagonalizing vectors in VT |
* (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB) |
* (Workspace: need 3*M + NRVT, prefer 3*M + NRVT*NB) |
* |
* |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
IF( WNTVA ) |
IF( WNTVA ) |
Line 3393
|
Line 3404
|
* |
* |
* If left singular vectors desired in A, generate left |
* If left singular vectors desired in A, generate left |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB) |
* (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3402
|
Line 3413
|
* |
* |
* If right singular vectors desired in A, generate right |
* If right singular vectors desired in A, generate right |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |