version 1.13, 2012/12/14 14:22:29
|
version 1.20, 2023/08/07 08:38:50
|
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 |
|
* |
* |
*> \ingroup doubleGEsing |
*> \ingroup doubleGEsing |
* |
* |
Line 211
|
Line 209
|
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 -- |
* -- 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 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU, JOBVT |
CHARACTER JOBU, JOBVT |
Line 314
|
Line 311
|
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 336
|
* 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 346
|
* 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 444
|
* |
* |
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 472
|
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 506
|
* 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 604
|
* |
* |
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 682
|
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 708
|
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 739
|
* 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 762
|
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 774
|
$ 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 784
|
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 800
|
* |
* |
* 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 809
|
* |
* |
* 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 830
|
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 863
|
* 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 886
|
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 899
|
$ 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 909
|
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 917
|
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 933
|
* 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 942
|
* |
* |
* 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 961
|
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 974
|
$ 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 984
|
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 1042
|
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 1055
|
$ 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 1065
|
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 1073
|
$ 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 1082
|
* |
* |
* 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 1103
|
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 1121
|
* |
* |
* 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 1169
|
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 1188
|
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 1201
|
$ 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 1212
|
* |
* |
* 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 1223
|
$ 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 1241
|
* 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 1268
|
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 1286
|
* |
* |
* 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 1350
|
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 1363
|
$ 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 1373
|
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 1383
|
$ 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 1400
|
* 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 1421
|
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 1445
|
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 1453
|
* |
* |
* 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 1507
|
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 1521
|
$ 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 1531
|
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 1539
|
$ 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 1548
|
* |
* |
* 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 1573
|
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 1591
|
* |
* |
* 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 1605
|
* |
* |
* 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 1640
|
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 1659
|
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 1684
|
* |
* |
* 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 1695
|
$ 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 1713
|
* 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 1743
|
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 1761
|
* |
* |
* 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 1775
|
* |
* |
* 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 1826
|
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 1850
|
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 1860
|
$ 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 1877
|
* 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 1902
|
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 1926
|
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 1934
|
* |
* |
* 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 1975
|
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 1984
|
* |
* |
* 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 1998
|
* |
* |
* 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 2008
|
* |
* |
* 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 2017
|
* |
* |
* 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 2079
|
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 2093
|
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 2101
|
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 2134
|
* 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 2160
|
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 2172
|
$ 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 2182
|
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 2198
|
* |
* |
* 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 2207
|
* |
* |
* 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 2228
|
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 2261
|
* 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 2287
|
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 2299
|
$ 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 2309
|
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 2317
|
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 2333
|
* 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 2342
|
* |
* |
* 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 2361
|
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 2373
|
$ 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 2383
|
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 2441
|
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 2454
|
$ 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 2464
|
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 2473
|
* |
* |
* 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 2482
|
* |
* |
* 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 2503
|
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 2513
|
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 2528
|
$ 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 2570
|
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 2589
|
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 2602
|
$ 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 2613
|
* |
* |
* 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 2624
|
$ 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 2632
|
$ 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 2642
|
* 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 2669
|
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 2691
|
$ 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 2749
|
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 2762
|
$ 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 2772
|
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 2782
|
$ 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 2790
|
$ 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 2799
|
* 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 2820
|
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 2843
|
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 2851
|
* |
* |
* 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 2885
|
* 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 2905
|
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 2919
|
$ 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 2929
|
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 2937
|
$ 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 2947
|
* |
* |
* 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 2972
|
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 2994
|
$ 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 3002
|
* |
* |
* 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 3025
|
* 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 3037
|
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 3056
|
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 3081
|
* |
* |
* 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 3092
|
$ 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 3100
|
$ 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 3110
|
* 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 3140
|
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 3162
|
$ 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 3170
|
* |
* |
* 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 3201
|
* 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 3221
|
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 3245
|
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 3255
|
$ 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 3271
|
* 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 3296
|
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 3319
|
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 3327
|
* |
* |
* 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 3368
|
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 3377
|
* |
* |
* 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 3387
|
* |
* |
* 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 3401
|
* |
* |
* 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 3410
|
* |
* |
* 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 ) |