diff options
| -rw-r--r-- | SRC/cgelss.f | 95 | ||||
| -rw-r--r-- | SRC/cgesvd.f | 331 | ||||
| -rw-r--r-- | SRC/cungbr.f | 17 | ||||
| -rw-r--r-- | SRC/dgelss.f | 102 | ||||
| -rw-r--r-- | SRC/dgesvd.f | 333 | ||||
| -rw-r--r-- | SRC/dorgbr.f | 17 | ||||
| -rw-r--r-- | SRC/sgelss.f | 94 | ||||
| -rw-r--r-- | SRC/sgesvd.f | 329 | ||||
| -rw-r--r-- | SRC/sorgbr.f | 17 | ||||
| -rw-r--r-- | SRC/zgelss.f | 95 | ||||
| -rw-r--r-- | SRC/zgesvd.f | 332 | ||||
| -rw-r--r-- | SRC/zungbr.f | 17 |
12 files changed, 986 insertions, 793 deletions
diff --git a/SRC/cgelss.f b/SRC/cgelss.f index 85c4f989..6212208c 100644 --- a/SRC/cgelss.f +++ b/SRC/cgelss.f @@ -206,10 +206,13 @@ INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK, $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, $ MAXWRK, MINMN, MINWRK, MM, MNTHR + INTEGER LWORK_CGEQRF, LWORK_CUNMQR, LWORK_CGEBRD, + $ LWORK_CUNMBR, LWORK_CUNGBR, LWORK_CUNMLQ, + $ LWORK_CGELQF REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR * .. * .. Local Arrays .. - COMPLEX VDUM( 1 ) + COMPLEX DUM( 1 ) * .. * .. External Subroutines .. EXTERNAL CBDSQR, CCOPY, CGEBRD, CGELQF, CGEMM, CGEMV, @@ -264,6 +267,13 @@ * Path 1a - overdetermined, with many more rows than * columns * +* Compute space needed for CGEQRF + CALL CGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) + LWORK_CGEQRF=DUM(1) +* Compute space needed for CUNMQR + CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, DUM(1), B, + $ LDB, DUM(1), -1, INFO ) + LWORK_CUNMQR=DUM(1) MM = N MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'CGEQRF', ' ', M, $ N, -1, -1 ) ) @@ -274,12 +284,22 @@ * * Path 1 - overdetermined or exactly determined * - MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, - $ 'CGEBRD', ' ', MM, N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'CUNMBR', - $ 'QLC', MM, NRHS, N, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, - $ 'CUNGBR', 'P', N, N, N, -1 ) ) +* Compute space needed for CGEBRD + CALL CGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_CGEBRD=DUM(1) +* Compute space needed for CUNMBR + CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_CUNMBR=DUM(1) +* Compute space needed for CUNGBR + CALL CUNGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_CUNGBR=DUM(1) +* Compute total workspace needed + MAXWRK = MAX( MAXWRK, 2*N + LWORK_CGEBRD ) + MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR ) + MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR ) MAXWRK = MAX( MAXWRK, N*NRHS ) MINWRK = 2*N + MAX( NRHS, M ) END IF @@ -290,31 +310,56 @@ * Path 2a - underdetermined, with many more columns * than rows * - MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1, - $ 'CGEBRD', ' ', M, M, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1, - $ 'CUNMBR', 'QLC', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1, - $ 'CUNGBR', 'P', M, M, M, -1 ) ) +* Compute space needed for CGELQF + CALL CGELQF( M, N, A, LDA, DUM(1), DUM(1), + $ -1, INFO ) + LWORK_CGELQF=DUM(1) +* Compute space needed for CGEBRD + CALL CGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_CGEBRD=DUM(1) +* Compute space needed for CUNMBR + CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_CUNMBR=DUM(1) +* Compute space needed for CUNGBR + CALL CUNGBR( 'P', M, M, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_CUNGBR=DUM(1) +* Compute space needed for CUNMLQ + CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_CUNMLQ=DUM(1) +* Compute total workspace needed + MAXWRK = M + LWORK_CGELQF + MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_CGEBRD ) + MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_CUNMBR ) + MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_CUNGBR ) IF( NRHS.GT.1 ) THEN MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) ELSE MAXWRK = MAX( MAXWRK, M*M + 2*M ) END IF - MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'CUNMLQ', - $ 'LC', N, NRHS, M, -1 ) ) + MAXWRK = MAX( MAXWRK, M + LWORK_CUNMLQ ) ELSE * * Path 2 - underdetermined * - MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'CGEBRD', ' ', M, - $ N, -1, -1 ) - MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'CUNMBR', - $ 'QLC', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'CUNGBR', - $ 'P', M, N, M, -1 ) ) +* Compute space needed for CGEBRD + CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_CGEBRD=DUM(1) +* Compute space needed for CUNMBR + CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_CUNMBR=DUM(1) +* Compute space needed for CUNGBR + CALL CUNGBR( 'P', M, N, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_CUNGBR=DUM(1) + MAXWRK = 2*M + LWORK_CGEBRD + MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR ) + MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR ) MAXWRK = MAX( MAXWRK, N*NRHS ) END IF END IF @@ -462,7 +507,7 @@ * (CWorkspace: none) * (RWorkspace: need BDSPAC) * - CALL CBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, + CALL CBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM, $ 1, B, LDB, RWORK( IRWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 @@ -659,7 +704,7 @@ * (CWorkspace: none) * (RWorkspace: need BDSPAC) * - CALL CBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, + CALL CBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM, $ 1, B, LDB, RWORK( IRWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 diff --git a/SRC/cgesvd.f b/SRC/cgesvd.f index c1158d67..dc4f8e11 100644 --- a/SRC/cgesvd.f +++ b/SRC/cgesvd.f @@ -245,6 +245,9 @@ $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, $ NRVT, WRKBL + INTEGER LWORK_CGEQRF, LWORK_CUNGQR_N, LWORK_CUNGQR_M, + $ LWORK_CGEBRD, LWORK_CUNGBR_P, LWORK_CUNGBR_Q, + $ LWORK_CGELQF, LWORK_CUNGLQ_N, LWORK_CUNGLQ_M REAL ANRM, BIGNUM, EPS, SMLNUM * .. * .. Local Arrays .. @@ -314,7 +317,28 @@ MAXWRK = 1 IF( M.GE.N .AND. MINMN.GT.0 ) THEN * -* Space needed for CBDSQR is BDSPAC = 5*N +* Space needed for ZBDSQR is BDSPAC = 5*N +* + MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) +* Compute space needed for CGEQRF + CALL CGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_CGEQRF=DUM(1) +* Compute space needed for CUNGQR + CALL CUNGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_CUNGQR_N=DUM(1) + CALL CUNGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_CUNGQR_M=DUM(1) +* Compute space needed for CGEBRD + CALL CGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_CGEBRD=DUM(1) +* Compute space needed for CUNGBR + CALL CUNGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_P=DUM(1) + CALL CUNGBR( 'Q', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_Q=DUM(1) * MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) IF( M.GE.MNTHR ) THEN @@ -322,25 +346,19 @@ * * Path 1 (M much larger than N, JOBU='N') * - MAXWRK = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) + MAXWRK = N + LWORK_CGEQRF + MAXWRK = MAX( MAXWRK, 2*N+LWORK_CGEBRD ) IF( WNTVO .OR. WNTVAS ) - $ MAXWRK = MAX( MAXWRK, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_P ) MINWRK = 3*N ELSE IF( WNTUO .AND. WNTVN ) THEN * * Path 2 (M much larger than N, JOBU='O', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) MINWRK = 2*N + M ELSE IF( WNTUO .AND. WNTVAS ) THEN @@ -348,43 +366,32 @@ * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) MINWRK = 2*N + M ELSE IF( WNTUS .AND. WNTVN ) THEN * * Path 4 (M much larger than N, JOBU='S', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUS .AND. WNTVO ) THEN * * Path 5 (M much larger than N, JOBU='S', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) MAXWRK = 2*N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUS .AND. WNTVAS ) THEN @@ -392,43 +399,32 @@ * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUA .AND. WNTVN ) THEN * * Path 7 (M much larger than N, JOBU='A', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_M ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUA .AND. WNTVO ) THEN * * Path 8 (M much larger than N, JOBU='A', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_M ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) MAXWRK = 2*N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUA .AND. WNTVAS ) THEN @@ -436,15 +432,11 @@ * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_CGEQRF + WRKBL = MAX( WRKBL, N+LWORK_CUNGQR_M ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_CUNGBR_P ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M END IF @@ -452,48 +444,70 @@ * * Path 10 (M at least N, but not much larger) * - MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTUS .OR. WNTUO ) - $ MAXWRK = MAX( MAXWRK, 2*N+N* - $ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) ) - IF( WNTUA ) - $ MAXWRK = MAX( MAXWRK, 2*N+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) - IF( .NOT.WNTVN ) - $ MAXWRK = MAX( MAXWRK, 2*N+( N-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) + CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_CGEBRD=DUM(1) + MAXWRK = 2*N + LWORK_CGEBRD + IF( WNTUS .OR. WNTUO ) THEN + CALL CUNGBR( 'Q', M, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_Q ) + END IF + IF( WNTUA ) THEN + CALL CUNGBR( 'Q', M, M, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_Q ) + END IF + IF( .NOT.WNTVN ) THEN + MAXWRK = MAX( MAXWRK, 2*N+LWORK_CUNGBR_P ) MINWRK = 2*N + M + END IF END IF ELSE IF( MINMN.GT.0 ) THEN * * Space needed for CBDSQR is BDSPAC = 5*M * MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) +* Compute space needed for CGELQF + CALL CGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_CGELQF=DUM(1) +* Compute space needed for CUNGLQ + CALL CUNGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR ) + LWORK_CUNGLQ_N=DUM(1) + CALL CUNGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_CUNGLQ_M=DUM(1) +* Compute space needed for CGEBRD + CALL CGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_CGEBRD=DUM(1) +* Compute space needed for CUNGBR P + CALL CUNGBR( 'P', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_P=DUM(1) +* Compute space needed for CUNGBR Q + CALL CUNGBR( 'Q', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_Q=DUM(1) IF( N.GE.MNTHR ) THEN IF( WNTVN ) THEN * * Path 1t(N much larger than M, JOBVT='N') * - MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) + MAXWRK = M + LWORK_CGELQF + MAXWRK = MAX( MAXWRK, 2*M+LWORK_CGEBRD ) IF( WNTUO .OR. WNTUAS ) - $ MAXWRK = MAX( MAXWRK, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_Q ) MINWRK = 3*M ELSE IF( WNTVO .AND. WNTUN ) THEN * * Path 2t(N much larger than M, JOBU='N', JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) MINWRK = 2*M + N ELSE IF( WNTVO .AND. WNTUAS ) THEN @@ -501,43 +515,32 @@ * Path 3t(N much larger than M, JOBU='S' or 'A', * JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) MINWRK = 2*M + N ELSE IF( WNTVS .AND. WNTUN ) THEN * * Path 4t(N much larger than M, JOBU='N', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVS .AND. WNTUO ) THEN * * Path 5t(N much larger than M, JOBU='O', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) MAXWRK = 2*M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVS .AND. WNTUAS ) THEN @@ -545,43 +548,32 @@ * Path 6t(N much larger than M, JOBU='S' or 'A', * JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVA .AND. WNTUN ) THEN * * Path 7t(N much larger than M, JOBU='N', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_N ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVA .AND. WNTUO ) THEN * * Path 8t(N much larger than M, JOBU='O', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_N ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) MAXWRK = 2*M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVA .AND. WNTUAS ) THEN @@ -589,15 +581,11 @@ * Path 9t(N much larger than M, JOBU='S' or 'A', * JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_CGELQF + WRKBL = MAX( WRKBL, M+LWORK_CUNGLQ_N ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_CUNGBR_Q ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N END IF @@ -605,18 +593,27 @@ * * Path 10t(N greater than M, but not much larger) * - MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTVS .OR. WNTVO ) - $ MAXWRK = MAX( MAXWRK, 2*M+M* - $ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) ) - IF( WNTVA ) - $ MAXWRK = MAX( MAXWRK, 2*M+N* - $ ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) ) - IF( .NOT.WNTUN ) - $ MAXWRK = MAX( MAXWRK, 2*M+( M-1 )* - $ ILAENV( 1, 'CUNGBR', 'Q', M, M, M, -1 ) ) + CALL CGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_CGEBRD=DUM(1) + MAXWRK = 2*M + LWORK_CGEBRD + IF( WNTVS .OR. WNTVO ) THEN +* Compute space needed for CUNGBR P + CALL CUNGBR( 'P', M, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_P ) + END IF + IF( WNTVA ) THEN + CALL CUNGBR( 'P', N, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_CUNGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_P ) + END IF + IF( .NOT.WNTUN ) THEN + MAXWRK = MAX( MAXWRK, 2*M+LWORK_CUNGBR_Q ) MINWRK = 2*M + N + END IF END IF END IF MAXWRK = MAX( MINWRK, MAXWRK ) diff --git a/SRC/cungbr.f b/SRC/cungbr.f index cd76c8e4..b516e943 100644 --- a/SRC/cungbr.f +++ b/SRC/cungbr.f @@ -218,12 +218,21 @@ * IF( INFO.EQ.0 ) THEN IF( WANTQ ) THEN - NB = ILAENV( 1, 'CUNGQR', ' ', M, N, K, -1 ) + IF( M.GE.K ) THEN + CALL CUNGQR( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL CUNGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF ELSE - NB = ILAENV( 1, 'CUNGLQ', ' ', M, N, K, -1 ) + IF( K.LT.N ) THEN + CALL CUNGLQ( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL CUNGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF END IF - LWKOPT = MAX( 1, MN )*NB - WORK( 1 ) = LWKOPT + LWKOPT = WORK( 1 ) END IF * IF( INFO.NE.0 ) THEN diff --git a/SRC/dgelss.f b/SRC/dgelss.f index 60aeb8e4..3e82b5af 100644 --- a/SRC/dgelss.f +++ b/SRC/dgelss.f @@ -196,10 +196,13 @@ INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL, $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, $ MAXWRK, MINMN, MINWRK, MM, MNTHR + INTEGER LWORK_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD, + $ LWORK_DORMBR, LWORK_DORGBR, LWORK_DORMLQ, + $ LWORK_DGELQF DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR * .. * .. Local Arrays .. - DOUBLE PRECISION VDUM( 1 ) + DOUBLE PRECISION DUM( 1 ) * .. * .. External Subroutines .. EXTERNAL DBDSQR, DCOPY, DGEBRD, DGELQF, DGEMM, DGEMV, @@ -252,11 +255,16 @@ * Path 1a - overdetermined, with many more rows than * columns * +* Compute space needed for DGEQRF + CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) + LWORK_DGEQRF=DUM(1) +* Compute space needed for DORMQR + CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B, + $ LDB, DUM(1), -1, INFO ) + LWORK_DORMQR=DUM(1) MM = N - MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'DGEQRF', ' ', M, - $ N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'DORMQR', 'LT', - $ M, NRHS, N, -1 ) ) + MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF ) + MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR ) END IF IF( M.GE.N ) THEN * @@ -265,12 +273,22 @@ * Compute workspace needed for DBDSQR * BDSPAC = MAX( 1, 5*N ) - MAXWRK = MAX( MAXWRK, 3*N + ( MM + N )*ILAENV( 1, - $ 'DGEBRD', ' ', MM, N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*N + NRHS*ILAENV( 1, 'DORMBR', - $ 'QLT', MM, NRHS, N, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*N + ( N - 1 )*ILAENV( 1, - $ 'DORGBR', 'P', N, N, N, -1 ) ) +* Compute space needed for DGEBRD + CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_DGEBRD=DUM(1) +* Compute space needed for DORMBR + CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_DORMBR=DUM(1) +* Compute space needed for DORGBR + CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_DORGBR=DUM(1) +* Compute total workspace needed + MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD ) + MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR ) + MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR ) MAXWRK = MAX( MAXWRK, BDSPAC ) MAXWRK = MAX( MAXWRK, N*NRHS ) MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC ) @@ -287,33 +305,57 @@ * Path 2a - underdetermined, with many more columns * than rows * - MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, - $ 'DGEBRD', ' ', M, M, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, - $ 'DORMBR', 'QLT', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, M*M + 4*M + - $ ( M - 1 )*ILAENV( 1, 'DORGBR', 'P', M, - $ M, M, -1 ) ) +* Compute space needed for DGELQF + CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), + $ -1, INFO ) + LWORK_DGELQF=DUM(1) +* Compute space needed for DGEBRD + CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_DGEBRD=DUM(1) +* Compute space needed for DORMBR + CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_DORMBR=DUM(1) +* Compute space needed for DORGBR + CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_DORGBR=DUM(1) +* Compute space needed for DORMLQ + CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_DORMLQ=DUM(1) +* Compute total workspace needed + MAXWRK = M + LWORK_DGELQF + MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD ) + MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORMBR ) + MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORGBR ) MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC ) IF( NRHS.GT.1 ) THEN MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) ELSE MAXWRK = MAX( MAXWRK, M*M + 2*M ) END IF - MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'DORMLQ', - $ 'LT', N, NRHS, M, -1 ) ) + MAXWRK = MAX( MAXWRK, M + LWORK_DORMLQ ) ELSE * * Path 2 - underdetermined * - MAXWRK = 3*M + ( N + M )*ILAENV( 1, 'DGEBRD', ' ', M, - $ N, -1, -1 ) - MAXWRK = MAX( MAXWRK, 3*M + NRHS*ILAENV( 1, 'DORMBR', - $ 'QLT', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*ILAENV( 1, 'DORGBR', - $ 'P', M, N, M, -1 ) ) +* Compute space needed for DGEBRD + CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_DGEBRD=DUM(1) +* Compute space needed for DORMBR + CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_DORMBR=DUM(1) +* Compute space needed for DORGBR + CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_DORGBR=DUM(1) + MAXWRK = 3*M + LWORK_DGEBRD + MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR ) + MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR ) MAXWRK = MAX( MAXWRK, BDSPAC ) MAXWRK = MAX( MAXWRK, N*NRHS ) END IF @@ -455,7 +497,7 @@ * compute right singular vectors in A * (Workspace: need BDSPAC) * - CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM, + CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, $ 1, B, LDB, WORK( IWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 @@ -638,7 +680,7 @@ * multiplying B by transpose of left singular vectors * (Workspace: need BDSPAC) * - CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM, + CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, $ 1, B, LDB, WORK( IWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 diff --git a/SRC/dgesvd.f b/SRC/dgesvd.f index 3ab77a08..84558d19 100644 --- a/SRC/dgesvd.f +++ b/SRC/dgesvd.f @@ -208,8 +208,8 @@ *> \ingroup doubleGEsing * * ===================================================================== - SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, - $ WORK, LWORK, INFO ) + SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, + $ VT, LDVT, WORK, LWORK, INFO ) * * -- LAPACK driver routine (version 3.3.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -238,6 +238,9 @@ $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, $ NRVT, WRKBL + INTEGER LWORK_DGEQRF, LWORK_DORGQR_N, LWORK_DORGQR_M, + $ LWORK_DGEBRD, LWORK_DORGBR_P, LWORK_DORGBR_Q, + $ LWORK_DGELQF, LWORK_DORGLQ_N, LWORK_DORGLQ_M DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM * .. * .. Local Arrays .. @@ -309,31 +312,46 @@ * MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) BDSPAC = 5*N +* Compute space needed for DGEQRF + CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_DGEQRF=DUM(1) +* Compute space needed for DORGQR + CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_DORGQR_N=DUM(1) + CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_DORGQR_M=DUM(1) +* Compute space needed for DGEBRD + CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_DGEBRD=DUM(1) +* Compute space needed for DORGBR P + CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_P=DUM(1) +* Compute space needed for DORGBR Q + CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_Q=DUM(1) +* IF( M.GE.MNTHR ) THEN IF( WNTUN ) THEN * * Path 1 (M much larger than N, JOBU='N') * - MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) + MAXWRK = N + LWORK_DGEQRF + MAXWRK = MAX( MAXWRK, 3*N+LWORK_DGEBRD ) IF( WNTVO .OR. WNTVAS ) - $ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P ) MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 4*N, BDSPAC ) ELSE IF( WNTUO .AND. WNTVN ) THEN * * Path 2 (M much larger than N, JOBU='O', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) MINWRK = MAX( 3*N+M, BDSPAC ) @@ -342,15 +360,11 @@ * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) MINWRK = MAX( 3*N+M, BDSPAC ) @@ -358,13 +372,10 @@ * * Path 4 (M much larger than N, JOBU='S', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -372,15 +383,11 @@ * * Path 5 (M much larger than N, JOBU='S', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -389,15 +396,11 @@ * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -405,13 +408,10 @@ * * Path 7 (M much larger than N, JOBU='A', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -419,15 +419,11 @@ * * Path 8 (M much larger than N, JOBU='A', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -436,15 +432,11 @@ * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_DGEQRF + WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -453,17 +445,25 @@ * * Path 10 (M at least N, but not much larger) * - MAXWRK = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTUS .OR. WNTUO ) - $ MAXWRK = MAX( MAXWRK, 3*N+N* - $ ILAENV( 1, 'DORGBR', 'Q', M, N, N, -1 ) ) - IF( WNTUA ) - $ MAXWRK = MAX( MAXWRK, 3*N+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, N, -1 ) ) - IF( .NOT.WNTVN ) - $ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* - $ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) + CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_DGEBRD=DUM(1) + MAXWRK = 3*N + LWORK_DGEBRD + IF( WNTUS .OR. WNTUO ) THEN + CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q ) + END IF + IF( WNTUA ) THEN + CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q ) + END IF + IF( .NOT.WNTVN ) THEN + MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P ) + END IF MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 3*N+M, BDSPAC ) END IF @@ -473,31 +473,45 @@ * MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) BDSPAC = 5*M +* Compute space needed for DGELQF + CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_DGELQF=DUM(1) +* Compute space needed for DORGLQ + CALL DORGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR ) + LWORK_DORGLQ_N=DUM(1) + CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_DORGLQ_M=DUM(1) +* Compute space needed for DGEBRD + CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_DGEBRD=DUM(1) +* Compute space needed for DORGBR P + CALL DORGBR( 'P', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_P=DUM(1) +* Compute space needed for DORGBR Q + CALL DORGBR( 'Q', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_Q=DUM(1) IF( N.GE.MNTHR ) THEN IF( WNTVN ) THEN * * Path 1t(N much larger than M, JOBVT='N') * - MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) + MAXWRK = M + LWORK_DGELQF + MAXWRK = MAX( MAXWRK, 3*M+LWORK_DGEBRD ) IF( WNTUO .OR. WNTUAS ) - $ MAXWRK = MAX( MAXWRK, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q ) MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 4*M, BDSPAC ) ELSE IF( WNTVO .AND. WNTUN ) THEN * * Path 2t(N much larger than M, JOBU='N', JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) MINWRK = MAX( 3*M+N, BDSPAC ) @@ -506,15 +520,11 @@ * Path 3t(N much larger than M, JOBU='S' or 'A', * JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) MINWRK = MAX( 3*M+N, BDSPAC ) @@ -522,13 +532,10 @@ * * Path 4t(N much larger than M, JOBU='N', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -536,15 +543,11 @@ * * Path 5t(N much larger than M, JOBU='O', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -553,15 +556,11 @@ * Path 6t(N much larger than M, JOBU='S' or 'A', * JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -569,13 +568,10 @@ * * Path 7t(N much larger than M, JOBU='N', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -583,15 +579,11 @@ * * Path 8t(N much larger than M, JOBU='O', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -600,15 +592,11 @@ * Path 9t(N much larger than M, JOBU='S' or 'A', * JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_DGELQF + WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -617,17 +605,26 @@ * * Path 10t(N greater than M, but not much larger) * - MAXWRK = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTVS .OR. WNTVO ) - $ MAXWRK = MAX( MAXWRK, 3*M+M* - $ ILAENV( 1, 'DORGBR', 'P', M, N, M, -1 ) ) - IF( WNTVA ) - $ MAXWRK = MAX( MAXWRK, 3*M+N* - $ ILAENV( 1, 'DORGBR', 'P', N, N, M, -1 ) ) - IF( .NOT.WNTUN ) - $ MAXWRK = MAX( MAXWRK, 3*M+( M-1 )* - $ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) + CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_DGEBRD=DUM(1) + MAXWRK = 3*M + LWORK_DGEBRD + IF( WNTVS .OR. WNTVO ) THEN +* Compute space needed for DORGBR P + CALL DORGBR( 'P', M, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P ) + END IF + IF( WNTVA ) THEN + CALL DORGBR( 'P', N, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_DORGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P ) + END IF + IF( .NOT.WNTUN ) THEN + MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q ) + END IF MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 3*M+N, BDSPAC ) END IF diff --git a/SRC/dorgbr.f b/SRC/dorgbr.f index 8cb93442..c6390b95 100644 --- a/SRC/dorgbr.f +++ b/SRC/dorgbr.f @@ -217,12 +217,21 @@ * IF( INFO.EQ.0 ) THEN IF( WANTQ ) THEN - NB = ILAENV( 1, 'DORGQR', ' ', M, N, K, -1 ) + IF( M.GE.K ) THEN + CALL DORGQR( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL DORGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF ELSE - NB = ILAENV( 1, 'DORGLQ', ' ', M, N, K, -1 ) + IF( K.LT.N ) THEN + CALL DORGLQ( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL DORGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF END IF - LWKOPT = MAX( 1, MN )*NB - WORK( 1 ) = LWKOPT + LWKOPT = WORK( 1 ) END IF * IF( INFO.NE.0 ) THEN diff --git a/SRC/sgelss.f b/SRC/sgelss.f index 1c2c2416..322bb82a 100644 --- a/SRC/sgelss.f +++ b/SRC/sgelss.f @@ -196,10 +196,12 @@ INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL, $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, $ MAXWRK, MINMN, MINWRK, MM, MNTHR + INTEGER LWORK_SGEQRF, LWORK_SORMQR, LWORK_SGEBRD, + $ LWORK_SORMBR, LWORK_SORGBR, LWORK_SORMLQ REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR * .. * .. Local Arrays .. - REAL VDUM( 1 ) + REAL DUM( 1 ) * .. * .. External Subroutines .. EXTERNAL SBDSQR, SCOPY, SGEBRD, SGELQF, SGEMM, SGEMV, @@ -252,11 +254,16 @@ * Path 1a - overdetermined, with many more rows than * columns * +* Compute space needed for SGEQRF + CALL SGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) + LWORK_SGEQRF=DUM(1) +* Compute space needed for SORMQR + CALL SORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B, + $ LDB, DUM(1), -1, INFO ) + LWORK_SORMQR=DUM(1) MM = N - MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'SGEQRF', ' ', M, - $ N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'SORMQR', 'LT', - $ M, NRHS, N, -1 ) ) + MAXWRK = MAX( MAXWRK, N + LWORK_SGEQRF ) + MAXWRK = MAX( MAXWRK, N + LWORK_SORMQR ) END IF IF( M.GE.N ) THEN * @@ -265,12 +272,22 @@ * Compute workspace needed for SBDSQR * BDSPAC = MAX( 1, 5*N ) - MAXWRK = MAX( MAXWRK, 3*N + ( MM + N )*ILAENV( 1, - $ 'SGEBRD', ' ', MM, N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*N + NRHS*ILAENV( 1, 'SORMBR', - $ 'QLT', MM, NRHS, N, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*N + ( N - 1 )*ILAENV( 1, - $ 'SORGBR', 'P', N, N, N, -1 ) ) +* Compute space needed for SGEBRD + CALL SGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_SGEBRD=DUM(1) +* Compute space needed for SORMBR + CALL SORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_SORMBR=DUM(1) +* Compute space needed for SORGBR + CALL SORGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_SORGBR=DUM(1) +* Compute total workspace needed + MAXWRK = MAX( MAXWRK, 3*N + LWORK_SGEBRD ) + MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORMBR ) + MAXWRK = MAX( MAXWRK, 3*N + LWORK_SORGBR ) MAXWRK = MAX( MAXWRK, BDSPAC ) MAXWRK = MAX( MAXWRK, N*NRHS ) MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC ) @@ -287,33 +304,54 @@ * Path 2a - underdetermined, with many more columns * than rows * +* Compute space needed for SGEBRD + CALL SGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_SGEBRD=DUM(1) +* Compute space needed for SORMBR + CALL SORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_SORMBR=DUM(1) +* Compute space needed for SORGBR + CALL SORGBR( 'P', M, M, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_SORGBR=DUM(1) +* Compute space needed for SORMLQ + CALL SORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_SORMLQ=DUM(1) +* Compute total workspace needed MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, $ -1 ) - MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, - $ 'SGEBRD', ' ', M, M, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, - $ 'SORMBR', 'QLT', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, M*M + 4*M + - $ ( M - 1 )*ILAENV( 1, 'SORGBR', 'P', M, - $ M, M, -1 ) ) + MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SGEBRD ) + MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORMBR ) + MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_SORGBR ) MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC ) IF( NRHS.GT.1 ) THEN MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) ELSE MAXWRK = MAX( MAXWRK, M*M + 2*M ) END IF - MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'SORMLQ', - $ 'LT', N, NRHS, M, -1 ) ) + MAXWRK = MAX( MAXWRK, M + LWORK_SORMLQ ) ELSE * * Path 2 - underdetermined * - MAXWRK = 3*M + ( N + M )*ILAENV( 1, 'SGEBRD', ' ', M, - $ N, -1, -1 ) - MAXWRK = MAX( MAXWRK, 3*M + NRHS*ILAENV( 1, 'SORMBR', - $ 'QLT', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*ILAENV( 1, 'SORGBR', - $ 'P', M, N, M, -1 ) ) +* Compute space needed for SGEBRD + CALL SGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_SGEBRD=DUM(1) +* Compute space needed for SORMBR + CALL SORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_SORMBR=DUM(1) +* Compute space needed for SORGBR + CALL SORGBR( 'P', M, N, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_SORGBR=DUM(1) + MAXWRK = 3*M + LWORK_SGEBRD + MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORMBR ) + MAXWRK = MAX( MAXWRK, 3*M + LWORK_SORGBR ) MAXWRK = MAX( MAXWRK, BDSPAC ) MAXWRK = MAX( MAXWRK, N*NRHS ) END IF @@ -455,7 +493,7 @@ * compute right singular vectors in A * (Workspace: need BDSPAC) * - CALL SBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM, + CALL SBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, $ 1, B, LDB, WORK( IWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 @@ -638,7 +676,7 @@ * multiplying B by transpose of left singular vectors * (Workspace: need BDSPAC) * - CALL SBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM, + CALL SBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM, $ 1, B, LDB, WORK( IWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 diff --git a/SRC/sgesvd.f b/SRC/sgesvd.f index 5ef13c28..8a945afe 100644 --- a/SRC/sgesvd.f +++ b/SRC/sgesvd.f @@ -238,6 +238,9 @@ $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, $ NRVT, WRKBL + INTEGER LWORK_SGEQRF, LWORK_SORGQR_N, LWORK_SORGQR_M, + $ LWORK_SGEBRD, LWORK_SORGBR_P, LWORK_SORGBR_Q, + $ LWORK_SGELQF, LWORK_SORGLQ_N, LWORK_SORGLQ_M REAL ANRM, BIGNUM, EPS, SMLNUM * .. * .. Local Arrays .. @@ -309,31 +312,46 @@ * MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) BDSPAC = 5*N +* Compute space needed for SGEQRF + CALL SGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_SGEQRF=DUM(1) +* Compute space needed for SORGQR + CALL SORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_SORGQR_N=DUM(1) + CALL SORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_SORGQR_M=DUM(1) +* Compute space needed for SGEBRD + CALL SGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_SGEBRD=DUM(1) +* Compute space needed for SORGBR P + CALL SORGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_P=DUM(1) +* Compute space needed for SORGBR Q + CALL SORGBR( 'Q', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_Q=DUM(1) +* IF( M.GE.MNTHR ) THEN IF( WNTUN ) THEN * * Path 1 (M much larger than N, JOBU='N') * - MAXWRK = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) + MAXWRK = N + LWORK_SGEQRF + MAXWRK = MAX( MAXWRK, 3*N+LWORK_SGEBRD ) IF( WNTVO .OR. WNTVAS ) - $ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 3*N+LWORK_SORGBR_P ) MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 4*N, BDSPAC ) ELSE IF( WNTUO .AND. WNTVN ) THEN * * Path 2 (M much larger than N, JOBU='O', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) MINWRK = MAX( 3*N+M, BDSPAC ) @@ -342,15 +360,11 @@ * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) MINWRK = MAX( 3*N+M, BDSPAC ) @@ -358,13 +372,10 @@ * * Path 4 (M much larger than N, JOBU='S', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -372,15 +383,11 @@ * * Path 5 (M much larger than N, JOBU='S', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -389,15 +396,11 @@ * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_N ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -405,13 +408,10 @@ * * Path 7 (M much larger than N, JOBU='A', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_M ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -419,15 +419,11 @@ * * Path 8 (M much larger than N, JOBU='A', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_M ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -436,15 +432,11 @@ * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+2*N* - $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_SGEQRF + WRKBL = MAX( WRKBL, N+LWORK_SORGQR_M ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_Q ) + WRKBL = MAX( WRKBL, 3*N+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = N*N + WRKBL MINWRK = MAX( 3*N+M, BDSPAC ) @@ -453,17 +445,25 @@ * * Path 10 (M at least N, but not much larger) * - MAXWRK = 3*N + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTUS .OR. WNTUO ) - $ MAXWRK = MAX( MAXWRK, 3*N+N* - $ ILAENV( 1, 'SORGBR', 'Q', M, N, N, -1 ) ) - IF( WNTUA ) - $ MAXWRK = MAX( MAXWRK, 3*N+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, N, -1 ) ) - IF( .NOT.WNTVN ) - $ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* - $ ILAENV( 1, 'SORGBR', 'P', N, N, N, -1 ) ) + CALL SGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_SGEBRD=DUM(1) + MAXWRK = 3*N + LWORK_SGEBRD + IF( WNTUS .OR. WNTUO ) THEN + CALL SORGBR( 'Q', M, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 3*N+LWORK_SORGBR_Q ) + END IF + IF( WNTUA ) THEN + CALL SORGBR( 'Q', M, M, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 3*N+LWORK_SORGBR_Q ) + END IF + IF( .NOT.WNTVN ) THEN + MAXWRK = MAX( MAXWRK, 3*N+LWORK_SORGBR_P ) + END IF MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 3*N+M, BDSPAC ) END IF @@ -473,31 +473,45 @@ * MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) BDSPAC = 5*M +* Compute space needed for SGELQF + CALL SGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_SGELQF=DUM(1) +* Compute space needed for SORGLQ + CALL SORGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR ) + LWORK_SORGLQ_N=DUM(1) + CALL SORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_SORGLQ_M=DUM(1) +* Compute space needed for SGEBRD + CALL SGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_SGEBRD=DUM(1) +* Compute space needed for SORGBR P + CALL SORGBR( 'P', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_P=DUM(1) +* Compute space needed for SORGBR Q + CALL SORGBR( 'Q', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_Q=DUM(1) IF( N.GE.MNTHR ) THEN IF( WNTVN ) THEN * * Path 1t(N much larger than M, JOBVT='N') * - MAXWRK = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) + MAXWRK = M + LWORK_SGELQF + MAXWRK = MAX( MAXWRK, 3*M+LWORK_SGEBRD ) IF( WNTUO .OR. WNTUAS ) - $ MAXWRK = MAX( MAXWRK, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 3*M+LWORK_SORGBR_Q ) MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 4*M, BDSPAC ) ELSE IF( WNTVO .AND. WNTUN ) THEN * * Path 2t(N much larger than M, JOBU='N', JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) MINWRK = MAX( 3*M+N, BDSPAC ) @@ -506,15 +520,11 @@ * Path 3t(N much larger than M, JOBU='S' or 'A', * JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) MINWRK = MAX( 3*M+N, BDSPAC ) @@ -522,13 +532,10 @@ * * Path 4t(N much larger than M, JOBU='N', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -536,15 +543,11 @@ * * Path 5t(N much larger than M, JOBU='O', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -554,15 +557,11 @@ * Path 6t(N much larger than M, JOBU='S' or 'A', * JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_M ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -570,13 +569,10 @@ * * Path 7t(N much larger than M, JOBU='N', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_N ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -584,15 +580,11 @@ * * Path 8t(N much larger than M, JOBU='O', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_N ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = 2*M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -601,15 +593,11 @@ * Path 9t(N much larger than M, JOBU='S' or 'A', * JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+2*M* - $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_SGELQF + WRKBL = MAX( WRKBL, M+LWORK_SORGLQ_N ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SGEBRD ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_P ) + WRKBL = MAX( WRKBL, 3*M+LWORK_SORGBR_Q ) WRKBL = MAX( WRKBL, BDSPAC ) MAXWRK = M*M + WRKBL MINWRK = MAX( 3*M+N, BDSPAC ) @@ -618,17 +606,26 @@ * * Path 10t(N greater than M, but not much larger) * - MAXWRK = 3*M + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTVS .OR. WNTVO ) - $ MAXWRK = MAX( MAXWRK, 3*M+M* - $ ILAENV( 1, 'SORGBR', 'P', M, N, M, -1 ) ) - IF( WNTVA ) - $ MAXWRK = MAX( MAXWRK, 3*M+N* - $ ILAENV( 1, 'SORGBR', 'P', N, N, M, -1 ) ) - IF( .NOT.WNTUN ) - $ MAXWRK = MAX( MAXWRK, 3*M+( M-1 )* - $ ILAENV( 1, 'SORGBR', 'Q', M, M, M, -1 ) ) + CALL SGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_SGEBRD=DUM(1) + MAXWRK = 3*M + LWORK_SGEBRD + IF( WNTVS .OR. WNTVO ) THEN +* Compute space needed for SORGBR P + CALL SORGBR( 'P', M, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 3*M+LWORK_SORGBR_P ) + END IF + IF( WNTVA ) THEN + CALL SORGBR( 'P', N, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_SORGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 3*M+LWORK_SORGBR_P ) + END IF + IF( .NOT.WNTUN ) THEN + MAXWRK = MAX( MAXWRK, 3*M+LWORK_SORGBR_Q ) + END IF MAXWRK = MAX( MAXWRK, BDSPAC ) MINWRK = MAX( 3*M+N, BDSPAC ) END IF diff --git a/SRC/sorgbr.f b/SRC/sorgbr.f index 8c66ffc2..8bb36a32 100644 --- a/SRC/sorgbr.f +++ b/SRC/sorgbr.f @@ -217,12 +217,21 @@ * IF( INFO.EQ.0 ) THEN IF( WANTQ ) THEN - NB = ILAENV( 1, 'SORGQR', ' ', M, N, K, -1 ) + IF( M.GE.K ) THEN + CALL SORGQR( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL SORGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF ELSE - NB = ILAENV( 1, 'SORGLQ', ' ', M, N, K, -1 ) + IF( K.LT.N ) THEN + CALL SORGLQ( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL SORGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF END IF - LWKOPT = MAX( 1, MN )*NB - WORK( 1 ) = LWKOPT + LWKOPT = WORK( 1 ) END IF * IF( INFO.NE.0 ) THEN diff --git a/SRC/zgelss.f b/SRC/zgelss.f index ca96cace..fed2fe26 100644 --- a/SRC/zgelss.f +++ b/SRC/zgelss.f @@ -206,10 +206,13 @@ INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK, $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, $ MAXWRK, MINMN, MINWRK, MM, MNTHR + INTEGER LWORK_ZGEQRF, LWORK_ZUNMQR, LWORK_ZGEBRD, + $ LWORK_ZUNMBR, LWORK_ZUNGBR, LWORK_ZUNMLQ, + $ LWORK_ZGELQF DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR * .. * .. Local Arrays .. - COMPLEX*16 VDUM( 1 ) + COMPLEX*16 DUM( 1 ) * .. * .. External Subroutines .. EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY, @@ -264,6 +267,13 @@ * Path 1a - overdetermined, with many more rows than * columns * +* Compute space needed for ZGEQRF + CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) + LWORK_ZGEQRF=DUM(1) +* Compute space needed for ZUNMQR + CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, DUM(1), B, + $ LDB, DUM(1), -1, INFO ) + LWORK_ZUNMQR=DUM(1) MM = N MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M, $ N, -1, -1 ) ) @@ -274,12 +284,22 @@ * * Path 1 - overdetermined or exactly determined * - MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, - $ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', - $ 'QLC', MM, NRHS, N, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, - $ 'ZUNGBR', 'P', N, N, N, -1 ) ) +* Compute space needed for ZGEBRD + CALL ZGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNMBR + CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_ZUNMBR=DUM(1) +* Compute space needed for ZUNGBR + CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_ZUNGBR=DUM(1) +* Compute total workspace needed + MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD ) + MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR ) + MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR ) MAXWRK = MAX( MAXWRK, N*NRHS ) MINWRK = 2*N + MAX( NRHS, M ) END IF @@ -290,31 +310,56 @@ * Path 2a - underdetermined, with many more columns * than rows * - MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1, - $ 'ZGEBRD', ' ', M, M, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1, - $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1, - $ 'ZUNGBR', 'P', M, M, M, -1 ) ) +* Compute space needed for ZGELQF + CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1), + $ -1, INFO ) + LWORK_ZGELQF=DUM(1) +* Compute space needed for ZGEBRD + CALL ZGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNMBR + CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_ZUNMBR=DUM(1) +* Compute space needed for ZUNGBR + CALL ZUNGBR( 'P', M, M, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_ZUNGBR=DUM(1) +* Compute space needed for ZUNMLQ + CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, DUM(1), + $ B, LDB, DUM(1), -1, INFO ) + LWORK_ZUNMLQ=DUM(1) +* Compute total workspace needed + MAXWRK = M + LWORK_ZGELQF + MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZGEBRD ) + MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNMBR ) + MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNGBR ) IF( NRHS.GT.1 ) THEN MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) ELSE MAXWRK = MAX( MAXWRK, M*M + 2*M ) END IF - MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ', - $ 'LC', N, NRHS, M, -1 ) ) + MAXWRK = MAX( MAXWRK, M + LWORK_ZUNMLQ ) ELSE * * Path 2 - underdetermined * - MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M, - $ N, -1, -1 ) - MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR', - $ 'QLC', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR', - $ 'P', M, N, M, -1 ) ) +* Compute space needed for ZGEBRD + CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNMBR + CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, A, LDA, + $ DUM(1), B, LDB, DUM(1), -1, INFO ) + LWORK_ZUNMBR=DUM(1) +* Compute space needed for ZUNGBR + CALL ZUNGBR( 'P', M, N, M, A, LDA, DUM(1), + $ DUM(1), -1, INFO ) + LWORK_ZUNGBR=DUM(1) + MAXWRK = 2*M + LWORK_ZGEBRD + MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR ) + MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR ) MAXWRK = MAX( MAXWRK, N*NRHS ) END IF END IF @@ -462,7 +507,7 @@ * (CWorkspace: none) * (RWorkspace: need BDSPAC) * - CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, + CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM, $ 1, B, LDB, RWORK( IRWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 @@ -659,7 +704,7 @@ * (CWorkspace: none) * (RWorkspace: need BDSPAC) * - CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, + CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM, $ 1, B, LDB, RWORK( IRWORK ), INFO ) IF( INFO.NE.0 ) $ GO TO 70 diff --git a/SRC/zgesvd.f b/SRC/zgesvd.f index 875345d7..7002279f 100644 --- a/SRC/zgesvd.f +++ b/SRC/zgesvd.f @@ -211,8 +211,8 @@ *> \ingroup complex16GEsing * * ===================================================================== - SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, - $ WORK, LWORK, RWORK, INFO ) + SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, + $ VT, LDVT, WORK, LWORK, RWORK, INFO ) * * -- LAPACK driver routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -245,6 +245,9 @@ $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, $ NRVT, WRKBL + INTEGER LWORK_ZGEQRF, LWORK_ZUNGQR_N, LWORK_ZUNGQR_M, + $ LWORK_ZGEBRD, LWORK_ZUNGBR_P, LWORK_ZUNGBR_Q, + $ LWORK_ZGELQF, LWORK_ZUNGLQ_N, LWORK_ZUNGLQ_M DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM * .. * .. Local Arrays .. @@ -317,30 +320,44 @@ * Space needed for ZBDSQR is BDSPAC = 5*N * MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) +* Compute space needed for ZGEQRF + CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_ZGEQRF=DUM(1) +* Compute space needed for ZUNGQR + CALL ZUNGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_ZUNGQR_N=DUM(1) + CALL ZUNGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_ZUNGQR_M=DUM(1) +* Compute space needed for ZGEBRD + CALL ZGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNGBR + CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_P=DUM(1) + CALL ZUNGBR( 'Q', N, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_Q=DUM(1) +* IF( M.GE.MNTHR ) THEN IF( WNTUN ) THEN * * Path 1 (M much larger than N, JOBU='N') * - MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) + MAXWRK = N + LWORK_ZGEQRF + MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZGEBRD ) IF( WNTVO .OR. WNTVAS ) - $ MAXWRK = MAX( MAXWRK, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P ) MINWRK = 3*N ELSE IF( WNTUO .AND. WNTVN ) THEN * * Path 2 (M much larger than N, JOBU='O', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) MINWRK = 2*N + M ELSE IF( WNTUO .AND. WNTVAS ) THEN @@ -348,43 +365,32 @@ * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) MINWRK = 2*N + M ELSE IF( WNTUS .AND. WNTVN ) THEN * * Path 4 (M much larger than N, JOBU='S', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUS .AND. WNTVO ) THEN * * Path 5 (M much larger than N, JOBU='S', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) MAXWRK = 2*N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUS .AND. WNTVAS ) THEN @@ -392,43 +398,32 @@ * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, - $ N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUA .AND. WNTVN ) THEN * * Path 7 (M much larger than N, JOBU='A', JOBVT='N') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUA .AND. WNTVO ) THEN * * Path 8 (M much larger than N, JOBU='A', JOBVT='O') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) MAXWRK = 2*N*N + WRKBL MINWRK = 2*N + M ELSE IF( WNTUA .AND. WNTVAS ) THEN @@ -436,15 +431,11 @@ * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or * 'A') * - WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M, - $ M, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+2*N* - $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) ) - WRKBL = MAX( WRKBL, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + WRKBL = N + LWORK_ZGEQRF + WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) + WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) MAXWRK = N*N + WRKBL MINWRK = 2*N + M END IF @@ -452,48 +443,70 @@ * * Path 10 (M at least N, but not much larger) * - MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTUS .OR. WNTUO ) - $ MAXWRK = MAX( MAXWRK, 2*N+N* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) ) - IF( WNTUA ) - $ MAXWRK = MAX( MAXWRK, 2*N+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) - IF( .NOT.WNTVN ) - $ MAXWRK = MAX( MAXWRK, 2*N+( N-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) + CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_ZGEBRD=DUM(1) + MAXWRK = 2*N + LWORK_ZGEBRD + IF( WNTUS .OR. WNTUO ) THEN + CALL ZUNGBR( 'Q', M, N, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q ) + END IF + IF( WNTUA ) THEN + CALL ZUNGBR( 'Q', M, M, N, A, LDA, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_Q=DUM(1) + MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q ) + END IF + IF( .NOT.WNTVN ) THEN + MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P ) MINWRK = 2*N + M + END IF END IF ELSE IF( MINMN.GT.0 ) THEN * * Space needed for ZBDSQR is BDSPAC = 5*M * MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) +* Compute space needed for ZGELQF + CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_ZGELQF=DUM(1) +* Compute space needed for ZUNGLQ + CALL ZUNGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR ) + LWORK_ZUNGLQ_N=DUM(1) + CALL ZUNGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) + LWORK_ZUNGLQ_M=DUM(1) +* Compute space needed for ZGEBRD + CALL ZGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNGBR P + CALL ZUNGBR( 'P', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_P=DUM(1) +* Compute space needed for ZUNGBR Q + CALL ZUNGBR( 'Q', M, M, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_Q=DUM(1) IF( N.GE.MNTHR ) THEN IF( WNTVN ) THEN * * Path 1t(N much larger than M, JOBVT='N') * - MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) + MAXWRK = M + LWORK_ZGELQF + MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZGEBRD ) IF( WNTUO .OR. WNTUAS ) - $ MAXWRK = MAX( MAXWRK, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + $ MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q ) MINWRK = 3*M ELSE IF( WNTVO .AND. WNTUN ) THEN * * Path 2t(N much larger than M, JOBU='N', JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) MINWRK = 2*M + N ELSE IF( WNTVO .AND. WNTUAS ) THEN @@ -501,43 +514,32 @@ * Path 3t(N much larger than M, JOBU='S' or 'A', * JOBVT='O') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) MINWRK = 2*M + N ELSE IF( WNTVS .AND. WNTUN ) THEN * * Path 4t(N much larger than M, JOBU='N', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVS .AND. WNTUO ) THEN * * Path 5t(N much larger than M, JOBU='O', JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) MAXWRK = 2*M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVS .AND. WNTUAS ) THEN @@ -545,43 +547,32 @@ * Path 6t(N much larger than M, JOBU='S' or 'A', * JOBVT='S') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVA .AND. WNTUN ) THEN * * Path 7t(N much larger than M, JOBU='N', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVA .AND. WNTUO ) THEN * * Path 8t(N much larger than M, JOBU='O', JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) MAXWRK = 2*M*M + WRKBL MINWRK = 2*M + N ELSE IF( WNTVA .AND. WNTUAS ) THEN @@ -589,15 +580,11 @@ * Path 9t(N much larger than M, JOBU='S' or 'A', * JOBVT='A') * - WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) - WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N, - $ N, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+2*M* - $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) ) - WRKBL = MAX( WRKBL, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + WRKBL = M + LWORK_ZGELQF + WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) + WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) MAXWRK = M*M + WRKBL MINWRK = 2*M + N END IF @@ -605,18 +592,27 @@ * * Path 10t(N greater than M, but not much larger) * - MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, - $ -1, -1 ) - IF( WNTVS .OR. WNTVO ) - $ MAXWRK = MAX( MAXWRK, 2*M+M* - $ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) ) - IF( WNTVA ) - $ MAXWRK = MAX( MAXWRK, 2*M+N* - $ ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) ) - IF( .NOT.WNTUN ) - $ MAXWRK = MAX( MAXWRK, 2*M+( M-1 )* - $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) ) + CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, IERR ) + LWORK_ZGEBRD=DUM(1) + MAXWRK = 2*M + LWORK_ZGEBRD + IF( WNTVS .OR. WNTVO ) THEN +* Compute space needed for ZUNGBR P + CALL ZUNGBR( 'P', M, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P ) + END IF + IF( WNTVA ) THEN + CALL ZUNGBR( 'P', N, N, M, A, N, DUM(1), + $ DUM(1), -1, IERR ) + LWORK_ZUNGBR_P=DUM(1) + MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P ) + END IF + IF( .NOT.WNTUN ) THEN + MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q ) MINWRK = 2*M + N + END IF END IF END IF MAXWRK = MAX( MAXWRK, MINWRK ) diff --git a/SRC/zungbr.f b/SRC/zungbr.f index dd2d9a7f..d7d3245f 100644 --- a/SRC/zungbr.f +++ b/SRC/zungbr.f @@ -218,12 +218,21 @@ * IF( INFO.EQ.0 ) THEN IF( WANTQ ) THEN - NB = ILAENV( 1, 'ZUNGQR', ' ', M, N, K, -1 ) + IF( M.GE.K ) THEN + CALL ZUNGQR( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL ZUNGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF ELSE - NB = ILAENV( 1, 'ZUNGLQ', ' ', M, N, K, -1 ) + IF( K.LT.N ) THEN + CALL ZUNGLQ( M, N, K, A, LDA, TAU, WORK, -1, IINFO ) + ELSE + CALL ZUNGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, + $ -1, IINFO ) + END IF END IF - LWKOPT = MAX( 1, MN )*NB - WORK( 1 ) = LWKOPT + LWKOPT = WORK( 1 ) END IF * IF( INFO.NE.0 ) THEN |