| @@ -51,8 +51,7 @@ float LAPACKE_clantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| } | } | ||||
| #endif | #endif | ||||
| /* Allocate memory for working array(s) */ | /* Allocate memory for working array(s) */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| work = (float*)LAPACKE_malloc( sizeof(float) * MAX(1,MAX(m,n)) ); | work = (float*)LAPACKE_malloc( sizeof(float) * MAX(1,MAX(m,n)) ); | ||||
| if( work == NULL ) { | if( work == NULL ) { | ||||
| info = LAPACK_WORK_MEMORY_ERROR; | info = LAPACK_WORK_MEMORY_ERROR; | ||||
| @@ -63,8 +62,7 @@ float LAPACKE_clantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| res = LAPACKE_clantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | res = LAPACKE_clantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | ||||
| work ); | work ); | ||||
| /* Release memory and exit */ | /* Release memory and exit */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| LAPACKE_free( work ); | LAPACKE_free( work ); | ||||
| } | } | ||||
| exit_level_0: | exit_level_0: | ||||
| @@ -51,8 +51,7 @@ double LAPACKE_dlantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| } | } | ||||
| #endif | #endif | ||||
| /* Allocate memory for working array(s) */ | /* Allocate memory for working array(s) */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| work = (double*)LAPACKE_malloc( sizeof(double) * MAX(1,MAX(m,n)) ); | work = (double*)LAPACKE_malloc( sizeof(double) * MAX(1,MAX(m,n)) ); | ||||
| if( work == NULL ) { | if( work == NULL ) { | ||||
| info = LAPACK_WORK_MEMORY_ERROR; | info = LAPACK_WORK_MEMORY_ERROR; | ||||
| @@ -63,8 +62,7 @@ double LAPACKE_dlantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| res = LAPACKE_dlantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | res = LAPACKE_dlantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | ||||
| work ); | work ); | ||||
| /* Release memory and exit */ | /* Release memory and exit */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| LAPACKE_free( work ); | LAPACKE_free( work ); | ||||
| } | } | ||||
| exit_level_0: | exit_level_0: | ||||
| @@ -38,10 +38,10 @@ double LAPACKE_dlantr_work( int matrix_layout, char norm, char uplo, | |||||
| const double* a, lapack_int lda, double* work ) | const double* a, lapack_int lda, double* work ) | ||||
| { | { | ||||
| lapack_int info = 0; | lapack_int info = 0; | ||||
| double res = 0.; | |||||
| double res = 0.; | |||||
| if( matrix_layout == LAPACK_COL_MAJOR ) { | if( matrix_layout == LAPACK_COL_MAJOR ) { | ||||
| /* Call LAPACK function and adjust info */ | /* Call LAPACK function and adjust info */ | ||||
| LAPACK_dlantr( &norm, &uplo, &diag, &m, &n, a, &lda, work ); | |||||
| res = LAPACK_dlantr( &norm, &uplo, &diag, &m, &n, a, &lda, work ); | |||||
| if( info < 0 ) { | if( info < 0 ) { | ||||
| info = info - 1; | info = info - 1; | ||||
| } | } | ||||
| @@ -74,11 +74,10 @@ lapack_int LAPACKE_dormbr_work( int matrix_layout, char vect, char side, | |||||
| } | } | ||||
| /* Allocate memory for temporary array(s) */ | /* Allocate memory for temporary array(s) */ | ||||
| if( LAPACKE_lsame( vect, 'q' ) ) { | if( LAPACKE_lsame( vect, 'q' ) ) { | ||||
| a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * k ); | |||||
| a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * MAX(1,k) ); | |||||
| } else { | } else { | ||||
| a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * nq ); | |||||
| a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * MAX(1,nq) ); | |||||
| } | } | ||||
| if( a_t == NULL ) { | if( a_t == NULL ) { | ||||
| info = LAPACK_TRANSPOSE_MEMORY_ERROR; | info = LAPACK_TRANSPOSE_MEMORY_ERROR; | ||||
| goto exit_level_0; | goto exit_level_0; | ||||
| @@ -89,11 +88,7 @@ lapack_int LAPACKE_dormbr_work( int matrix_layout, char vect, char side, | |||||
| goto exit_level_1; | goto exit_level_1; | ||||
| } | } | ||||
| /* Transpose input matrices */ | /* Transpose input matrices */ | ||||
| if( LAPACKE_lsame( vect, 'q' ) ) { | |||||
| LAPACKE_dge_trans( matrix_layout, nq, k, a, lda, a_t, lda_t ); | |||||
| } else { | |||||
| LAPACKE_dge_trans( matrix_layout, k, nq, a, lda, a_t, lda_t ); | |||||
| } | |||||
| LAPACKE_dge_trans( matrix_layout, r, MIN(nq,k), a, lda, a_t, lda_t ); | |||||
| LAPACKE_dge_trans( matrix_layout, m, n, c, ldc, c_t, ldc_t ); | LAPACKE_dge_trans( matrix_layout, m, n, c, ldc, c_t, ldc_t ); | ||||
| /* Call LAPACK function and adjust info */ | /* Call LAPACK function and adjust info */ | ||||
| LAPACK_dormbr( &vect, &side, &trans, &m, &n, &k, a_t, &lda_t, tau, c_t, | LAPACK_dormbr( &vect, &side, &trans, &m, &n, &k, a_t, &lda_t, tau, c_t, | ||||
| @@ -87,12 +87,7 @@ lapack_int LAPACKE_dormlq_work( int matrix_layout, char side, char trans, | |||||
| goto exit_level_1; | goto exit_level_1; | ||||
| } | } | ||||
| /* Transpose input matrices */ | /* Transpose input matrices */ | ||||
| if( LAPACKE_lsame( side, 'l' ) ){ | |||||
| LAPACKE_dge_trans( matrix_layout, k, m, a, lda, a_t, lda_t ); | |||||
| } else { | |||||
| LAPACKE_dge_trans( matrix_layout, k, n, a, lda, a_t, lda_t ); | |||||
| } | |||||
| LAPACKE_dge_trans( matrix_layout, k, m, a, lda, a_t, lda_t ); | |||||
| LAPACKE_dge_trans( matrix_layout, m, n, c, ldc, c_t, ldc_t ); | LAPACKE_dge_trans( matrix_layout, m, n, c, ldc, c_t, ldc_t ); | ||||
| /* Call LAPACK function and adjust info */ | /* Call LAPACK function and adjust info */ | ||||
| LAPACK_dormlq( &side, &trans, &m, &n, &k, a_t, &lda_t, tau, c_t, &ldc_t, | LAPACK_dormlq( &side, &trans, &m, &n, &k, a_t, &lda_t, tau, c_t, &ldc_t, | ||||
| @@ -51,8 +51,7 @@ float LAPACKE_slantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| } | } | ||||
| #endif | #endif | ||||
| /* Allocate memory for working array(s) */ | /* Allocate memory for working array(s) */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| work = (float*)LAPACKE_malloc( sizeof(float) * MAX(1,MAX(m,n)) ); | work = (float*)LAPACKE_malloc( sizeof(float) * MAX(1,MAX(m,n)) ); | ||||
| if( work == NULL ) { | if( work == NULL ) { | ||||
| info = LAPACK_WORK_MEMORY_ERROR; | info = LAPACK_WORK_MEMORY_ERROR; | ||||
| @@ -63,8 +62,7 @@ float LAPACKE_slantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| res = LAPACKE_slantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | res = LAPACKE_slantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | ||||
| work ); | work ); | ||||
| /* Release memory and exit */ | /* Release memory and exit */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| LAPACKE_free( work ); | LAPACKE_free( work ); | ||||
| } | } | ||||
| exit_level_0: | exit_level_0: | ||||
| @@ -41,7 +41,7 @@ float LAPACKE_slantr_work( int matrix_layout, char norm, char uplo, | |||||
| float res = 0.; | float res = 0.; | ||||
| if( matrix_layout == LAPACK_COL_MAJOR ) { | if( matrix_layout == LAPACK_COL_MAJOR ) { | ||||
| /* Call LAPACK function and adjust info */ | /* Call LAPACK function and adjust info */ | ||||
| LAPACK_slantr( &norm, &uplo, &diag, &m, &n, a, &lda, work ); | |||||
| res = LAPACK_slantr( &norm, &uplo, &diag, &m, &n, a, &lda, work ); | |||||
| if( info < 0 ) { | if( info < 0 ) { | ||||
| info = info - 1; | info = info - 1; | ||||
| } | } | ||||
| @@ -73,8 +73,11 @@ lapack_int LAPACKE_sormbr_work( int matrix_layout, char vect, char side, | |||||
| return (info < 0) ? (info - 1) : info; | return (info < 0) ? (info - 1) : info; | ||||
| } | } | ||||
| /* Allocate memory for temporary array(s) */ | /* Allocate memory for temporary array(s) */ | ||||
| a_t = (float*) | |||||
| LAPACKE_malloc( sizeof(float) * lda_t * MAX(1,MIN(nq,k)) ); | |||||
| if( LAPACKE_lsame( vect, 'q' ) ) { | |||||
| a_t = (float*)LAPACKE_malloc( sizeof(float) * lda_t * MAX(1,k) ); | |||||
| } else { | |||||
| a_t = (float*)LAPACKE_malloc( sizeof(float) * lda_t * MAX(1,nq) ); | |||||
| } | |||||
| if( a_t == NULL ) { | if( a_t == NULL ) { | ||||
| info = LAPACK_TRANSPOSE_MEMORY_ERROR; | info = LAPACK_TRANSPOSE_MEMORY_ERROR; | ||||
| goto exit_level_0; | goto exit_level_0; | ||||
| @@ -72,7 +72,11 @@ lapack_int LAPACKE_sormlq_work( int matrix_layout, char side, char trans, | |||||
| return (info < 0) ? (info - 1) : info; | return (info < 0) ? (info - 1) : info; | ||||
| } | } | ||||
| /* Allocate memory for temporary array(s) */ | /* Allocate memory for temporary array(s) */ | ||||
| a_t = (float*)LAPACKE_malloc( sizeof(float) * lda_t * MAX(1,m) ); | |||||
| if( LAPACKE_lsame( side, 'l' ) ) { | |||||
| a_t = (float*)LAPACKE_malloc( sizeof(float) * lda_t * MAX(1,m) ); | |||||
| } else { | |||||
| a_t = (float*)LAPACKE_malloc( sizeof(float) * lda_t * MAX(1,n) ); | |||||
| } | |||||
| if( a_t == NULL ) { | if( a_t == NULL ) { | ||||
| info = LAPACK_TRANSPOSE_MEMORY_ERROR; | info = LAPACK_TRANSPOSE_MEMORY_ERROR; | ||||
| goto exit_level_0; | goto exit_level_0; | ||||
| @@ -51,8 +51,7 @@ double LAPACKE_zlantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| } | } | ||||
| #endif | #endif | ||||
| /* Allocate memory for working array(s) */ | /* Allocate memory for working array(s) */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| work = (double*)LAPACKE_malloc( sizeof(double) * MAX(1,MAX(m,n)) ); | work = (double*)LAPACKE_malloc( sizeof(double) * MAX(1,MAX(m,n)) ); | ||||
| if( work == NULL ) { | if( work == NULL ) { | ||||
| info = LAPACK_WORK_MEMORY_ERROR; | info = LAPACK_WORK_MEMORY_ERROR; | ||||
| @@ -63,8 +62,7 @@ double LAPACKE_zlantr( int matrix_layout, char norm, char uplo, char diag, | |||||
| res = LAPACKE_zlantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | res = LAPACKE_zlantr_work( matrix_layout, norm, uplo, diag, m, n, a, lda, | ||||
| work ); | work ); | ||||
| /* Release memory and exit */ | /* Release memory and exit */ | ||||
| if( LAPACKE_lsame( norm, 'i' ) || LAPACKE_lsame( norm, '1' ) || | |||||
| LAPACKE_lsame( norm, 'O' ) ) { | |||||
| if( LAPACKE_lsame( norm, 'i' ) ) { | |||||
| LAPACKE_free( work ); | LAPACKE_free( work ); | ||||
| } | } | ||||
| exit_level_0: | exit_level_0: | ||||
| @@ -39,7 +39,7 @@ double LAPACKE_zlantr_work( int matrix_layout, char norm, char uplo, | |||||
| double* work ) | double* work ) | ||||
| { | { | ||||
| lapack_int info = 0; | lapack_int info = 0; | ||||
| double res = 0.; | |||||
| double res = 0.; | |||||
| if( matrix_layout == LAPACK_COL_MAJOR ) { | if( matrix_layout == LAPACK_COL_MAJOR ) { | ||||
| /* Call LAPACK function and adjust info */ | /* Call LAPACK function and adjust info */ | ||||
| res = LAPACK_zlantr( &norm, &uplo, &diag, &m, &n, a, &lda, work ); | res = LAPACK_zlantr( &norm, &uplo, &diag, &m, &n, a, &lda, work ); | ||||
| @@ -405,9 +405,9 @@ | |||||
| $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | ||||
| END IF | END IF | ||||
| * | * | ||||
| * If INFO > 0 from CHSEQR, then quit | |||||
| * If INFO .NE. 0 from CHSEQR, then quit | |||||
| * | * | ||||
| IF( INFO.GT.0 ) | |||||
| IF( INFO.NE.0 ) | |||||
| $ GO TO 50 | $ GO TO 50 | ||||
| * | * | ||||
| IF( WANTVL .OR. WANTVR ) THEN | IF( WANTVL .OR. WANTVR ) THEN | ||||
| @@ -170,7 +170,7 @@ | |||||
| *> vectors, stored columnwise) as specified by RANGE; if | *> vectors, stored columnwise) as specified by RANGE; if | ||||
| *> JOBU = 'N', U is not referenced. | *> JOBU = 'N', U is not referenced. | ||||
| *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | ||||
| *> the exact value of NS is not known ILQFin advance and an upper | |||||
| *> the exact value of NS is not known in advance and an upper | |||||
| *> bound must be used. | *> bound must be used. | ||||
| *> \endverbatim | *> \endverbatim | ||||
| *> | *> | ||||
| @@ -294,8 +294,8 @@ | |||||
| CHARACTER JOBZ, RNGTGK | CHARACTER JOBZ, RNGTGK | ||||
| LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT | LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT | ||||
| INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, | INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, | ||||
| $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK, | |||||
| $ J, K, MAXWRK, MINMN, MINWRK, MNTHR | |||||
| $ ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ, | |||||
| $ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR | |||||
| REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM | REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM | ||||
| * .. | * .. | ||||
| * .. Local Arrays .. | * .. Local Arrays .. | ||||
| @@ -367,8 +367,14 @@ | |||||
| IF( INFO.EQ.0 ) THEN | IF( INFO.EQ.0 ) THEN | ||||
| IF( WANTU .AND. LDU.LT.M ) THEN | IF( WANTU .AND. LDU.LT.M ) THEN | ||||
| INFO = -15 | INFO = -15 | ||||
| ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN | |||||
| INFO = -16 | |||||
| ELSE IF( WANTVT ) THEN | |||||
| IF( INDS ) THEN | |||||
| IF( LDVT.LT.IU-IL+1 ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| ELSE IF( LDVT.LT.MINMN ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -390,18 +396,24 @@ | |||||
| * | * | ||||
| * Path 1 (M much larger than N) | * Path 1 (M much larger than N) | ||||
| * | * | ||||
| MAXWRK = N + N* | |||||
| $ ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) | |||||
| MAXWRK = MAX( MAXWRK, N*N + N + 2*N* | |||||
| $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) | |||||
| MINWRK = N*(N+4) | |||||
| MINWRK = N*(N+5) | |||||
| MAXWRK = N + N*ILAENV(1,'CGEQRF',' ',M,N,-1,-1) | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ N*N+2*N+2*N*ILAENV(1,'CGEBRD',' ',N,N,-1,-1)) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ N*N+2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1)) | |||||
| END IF | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2 (M at least N, but not much larger) | * Path 2 (M at least N, but not much larger) | ||||
| * | * | ||||
| MAXWRK = 2*N + ( M+N )* | |||||
| $ ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 ) | |||||
| MINWRK = 2*N + M | |||||
| MINWRK = 3*N + M | |||||
| MAXWRK = 2*N + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ 2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1)) | |||||
| END IF | |||||
| END IF | END IF | ||||
| ELSE | ELSE | ||||
| MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) | MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) | ||||
| @@ -409,18 +421,25 @@ | |||||
| * | * | ||||
| * Path 1t (N much larger than M) | * Path 1t (N much larger than M) | ||||
| * | * | ||||
| MAXWRK = M + M* | |||||
| $ ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) | |||||
| MAXWRK = MAX( MAXWRK, M*M + M + 2*M* | |||||
| $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) | |||||
| MINWRK = M*(M+4) | |||||
| MINWRK = M*(M+5) | |||||
| MAXWRK = M + M*ILAENV(1,'CGELQF',' ',M,N,-1,-1) | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ M*M+2*M+2*M*ILAENV(1,'CGEBRD',' ',M,M,-1,-1)) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ M*M+2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1)) | |||||
| END IF | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2t (N greater than M, but not much larger) | * Path 2t (N greater than M, but not much larger) | ||||
| * | * | ||||
| MAXWRK = M*(M*2+19) + ( M+N )* | |||||
| $ ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 ) | |||||
| MINWRK = 2*M + N | |||||
| * | |||||
| MINWRK = 3*M + N | |||||
| MAXWRK = 2*M + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ 2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1)) | |||||
| END IF | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -518,14 +537,14 @@ | |||||
| CALL CGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), | CALL CGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), | ||||
| $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | ||||
| $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + N*(N*2+1) | |||||
| ITEMPR = ITGKZ + N*(N*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*N*N+14*N) | * (Workspace: need 2*N*N+14*N) | ||||
| * | * | ||||
| CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -539,7 +558,7 @@ | |||||
| END DO | END DO | ||||
| K = K + N | K = K + N | ||||
| END DO | END DO | ||||
| CALL CLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU ) | |||||
| CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) | |||||
| * | * | ||||
| * Call CUNMBR to compute QB*UB. | * Call CUNMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -594,14 +613,14 @@ | |||||
| CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | ||||
| $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | ||||
| $ LWORK-ITEMP+1, INFO ) | $ LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + N*(N*2+1) | |||||
| ITEMPR = ITGKZ + N*(N*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*N*N+14*N) | * (Workspace: need 2*N*N+14*N) | ||||
| * | * | ||||
| CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -615,7 +634,7 @@ | |||||
| END DO | END DO | ||||
| K = K + N | K = K + N | ||||
| END DO | END DO | ||||
| CALL CLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU ) | |||||
| CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) | |||||
| * | * | ||||
| * Call CUNMBR to compute QB*UB. | * Call CUNMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -681,14 +700,14 @@ | |||||
| CALL CGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), | CALL CGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), | ||||
| $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | ||||
| $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + M*(M*2+1) | |||||
| ITEMPR = ITGKZ + M*(M*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*M*M+14*M) | * (Workspace: need 2*M*M+14*M) | ||||
| * | * | ||||
| CALL SBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), | CALL SBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -722,7 +741,7 @@ | |||||
| END DO | END DO | ||||
| K = K + M | K = K + M | ||||
| END DO | END DO | ||||
| CALL CLASET( 'A', M, N-M, CZERO, CZERO, | |||||
| CALL CLASET( 'A', NS, N-M, CZERO, CZERO, | |||||
| $ VT( 1,M+1 ), LDVT ) | $ VT( 1,M+1 ), LDVT ) | ||||
| * | * | ||||
| * Call CUNMBR to compute (VB**T)*(PB**T) | * Call CUNMBR to compute (VB**T)*(PB**T) | ||||
| @@ -758,14 +777,14 @@ | |||||
| CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | ||||
| $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | ||||
| $ LWORK-ITEMP+1, INFO ) | $ LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + M*(M*2+1) | |||||
| ITEMPR = ITGKZ + M*(M*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*M*M+14*M) | * (Workspace: need 2*M*M+14*M) | ||||
| * | * | ||||
| CALL SBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), | CALL SBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -799,7 +818,7 @@ | |||||
| END DO | END DO | ||||
| K = K + M | K = K + M | ||||
| END DO | END DO | ||||
| CALL CLASET( 'A', M, N-M, CZERO, CZERO, | |||||
| CALL CLASET( 'A', NS, N-M, CZERO, CZERO, | |||||
| $ VT( 1,M+1 ), LDVT ) | $ VT( 1,M+1 ), LDVT ) | ||||
| * | * | ||||
| * Call CUNMBR to compute VB**T * PB**T | * Call CUNMBR to compute VB**T * PB**T | ||||
| @@ -145,15 +145,33 @@ | |||||
| INTRINSIC ABS, CMPLX, MAX | INTRINSIC ABS, CMPLX, MAX | ||||
| * .. | * .. | ||||
| * .. Executable Statements .. | * .. Executable Statements .. | ||||
| * | |||||
| INFO = 0 | |||||
| * | |||||
| * Quick return if possible | |||||
| * | |||||
| IF( N.EQ.0 ) | |||||
| $ RETURN | |||||
| * | * | ||||
| * Set constants to control overflow | * Set constants to control overflow | ||||
| * | * | ||||
| INFO = 0 | |||||
| EPS = SLAMCH( 'P' ) | EPS = SLAMCH( 'P' ) | ||||
| SMLNUM = SLAMCH( 'S' ) / EPS | SMLNUM = SLAMCH( 'S' ) / EPS | ||||
| BIGNUM = ONE / SMLNUM | BIGNUM = ONE / SMLNUM | ||||
| CALL SLABAD( SMLNUM, BIGNUM ) | CALL SLABAD( SMLNUM, BIGNUM ) | ||||
| * | * | ||||
| * Handle the case N=1 by itself | |||||
| * | |||||
| IF( N.EQ.1 ) THEN | |||||
| IPIV( 1 ) = 1 | |||||
| JPIV( 1 ) = 1 | |||||
| IF( ABS( A( 1, 1 ) ).LT.SMLNUM ) THEN | |||||
| INFO = 1 | |||||
| A( 1, 1 ) = CMPLX( SMLNUM, ZERO ) | |||||
| END IF | |||||
| RETURN | |||||
| END IF | |||||
| * | |||||
| * Factorize A using complete pivoting. | * Factorize A using complete pivoting. | ||||
| * Set pivots less than SMIN to SMIN | * Set pivots less than SMIN to SMIN | ||||
| * | * | ||||
| @@ -339,16 +339,16 @@ | |||||
| $ LDVL, VR, LDVR, WORK, -1, IERR ) | $ LDVL, VR, LDVR, WORK, -1, IERR ) | ||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| CALL CHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | CALL CHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | ||||
| $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, | |||||
| $ -1, WORK, IERR ) | |||||
| $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1, | |||||
| $ RWORK, IERR ) | |||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| ELSE | ELSE | ||||
| CALL CGGHD3( 'N', 'N', N, 1, N, A, LDA, B, LDB, VL, LDVL, | CALL CGGHD3( 'N', 'N', N, 1, N, A, LDA, B, LDB, VL, LDVL, | ||||
| $ VR, LDVR, WORK, -1, IERR ) | $ VR, LDVR, WORK, -1, IERR ) | ||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| CALL CHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | CALL CHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | ||||
| $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, | |||||
| $ -1, WORK, IERR ) | |||||
| $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1, | |||||
| $ RWORK, IERR ) | |||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| END IF | END IF | ||||
| WORK( 1 ) = CMPLX( LWKOPT ) | WORK( 1 ) = CMPLX( LWKOPT ) | ||||
| @@ -418,9 +418,9 @@ | |||||
| $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | ||||
| END IF | END IF | ||||
| * | * | ||||
| * If INFO > 0 from DHSEQR, then quit | |||||
| * If INFO .NE. 0 from DHSEQR, then quit | |||||
| * | * | ||||
| IF( INFO.GT.0 ) | |||||
| IF( INFO.NE.0 ) | |||||
| $ GO TO 50 | $ GO TO 50 | ||||
| * | * | ||||
| IF( WANTVL .OR. WANTVR ) THEN | IF( WANTVL .OR. WANTVR ) THEN | ||||
| @@ -169,7 +169,7 @@ | |||||
| *> vectors, stored columnwise) as specified by RANGE; if | *> vectors, stored columnwise) as specified by RANGE; if | ||||
| *> JOBU = 'N', U is not referenced. | *> JOBU = 'N', U is not referenced. | ||||
| *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | ||||
| *> the exact value of NS is not known ILQFin advance and an upper | |||||
| *> the exact value of NS is not known in advance and an upper | |||||
| *> bound must be used. | *> bound must be used. | ||||
| *> \endverbatim | *> \endverbatim | ||||
| *> | *> | ||||
| @@ -357,8 +357,14 @@ | |||||
| IF( INFO.EQ.0 ) THEN | IF( INFO.EQ.0 ) THEN | ||||
| IF( WANTU .AND. LDU.LT.M ) THEN | IF( WANTU .AND. LDU.LT.M ) THEN | ||||
| INFO = -15 | INFO = -15 | ||||
| ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN | |||||
| INFO = -16 | |||||
| ELSE IF( WANTVT ) THEN | |||||
| IF( INDS ) THEN | |||||
| IF( LDVT.LT.IU-IL+1 ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| ELSE IF( LDVT.LT.MINMN ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -380,18 +386,34 @@ | |||||
| * | * | ||||
| * Path 1 (M much larger than N) | * Path 1 (M much larger than N) | ||||
| * | * | ||||
| MAXWRK = N*(N*2+16) + | |||||
| MAXWRK = N + | |||||
| $ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) | $ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) | ||||
| MAXWRK = MAX( MAXWRK, N*(N*2+20) + 2*N* | |||||
| MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N* | |||||
| $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) | $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) | ||||
| MINWRK = N*(N*2+21) | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*3+6)+N* | |||||
| $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*3+6)+N* | |||||
| $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = N*(N*3+20) | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2 (M at least N, but not much larger) | * Path 2 (M at least N, but not much larger) | ||||
| * | * | ||||
| MAXWRK = N*(N*2+19) + ( M+N )* | |||||
| MAXWRK = 4*N + ( M+N )* | |||||
| $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) | $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) | ||||
| MINWRK = N*(N*2+20) + M | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*2+5)+N* | |||||
| $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*2+5)+N* | |||||
| $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = MAX(N*(N*2+19),4*N+M) | |||||
| END IF | END IF | ||||
| ELSE | ELSE | ||||
| MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) | MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) | ||||
| @@ -399,18 +421,34 @@ | |||||
| * | * | ||||
| * Path 1t (N much larger than M) | * Path 1t (N much larger than M) | ||||
| * | * | ||||
| MAXWRK = M*(M*2+16) + | |||||
| MAXWRK = M + | |||||
| $ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) | $ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) | ||||
| MAXWRK = MAX( MAXWRK, M*(M*2+20) + 2*M* | |||||
| MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M* | |||||
| $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) | $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) | ||||
| MINWRK = M*(M*2+21) | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*3+6)+M* | |||||
| $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*3+6)+M* | |||||
| $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = M*(M*3+20) | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2t (N greater than M, but not much larger) | |||||
| * Path 2t (N at least M, but not much larger) | |||||
| * | * | ||||
| MAXWRK = M*(M*2+19) + ( M+N )* | |||||
| MAXWRK = 4*M + ( M+N )* | |||||
| $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) | $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) | ||||
| MINWRK = M*(M*2+20) + N | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*2+5)+M* | |||||
| $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*2+5)+M* | |||||
| $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = MAX(M*(M*2+19),4*M+N) | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -522,7 +560,7 @@ | |||||
| CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | ||||
| J = J + N*2 | J = J + N*2 | ||||
| END DO | END DO | ||||
| CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| * | * | ||||
| * Call DORMBR to compute QB*UB. | * Call DORMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -591,7 +629,7 @@ | |||||
| CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | ||||
| J = J + N*2 | J = J + N*2 | ||||
| END DO | END DO | ||||
| CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| * | * | ||||
| * Call DORMBR to compute QB*UB. | * Call DORMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -687,7 +725,7 @@ | |||||
| CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | ||||
| J = J + M*2 | J = J + M*2 | ||||
| END DO | END DO | ||||
| CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT ) | |||||
| CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) | |||||
| * | * | ||||
| * Call DORMBR to compute (VB**T)*(PB**T) | * Call DORMBR to compute (VB**T)*(PB**T) | ||||
| * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | ||||
| @@ -756,7 +794,7 @@ | |||||
| CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | ||||
| J = J + M*2 | J = J + M*2 | ||||
| END DO | END DO | ||||
| CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT ) | |||||
| CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) | |||||
| * | * | ||||
| * Call DORMBR to compute VB**T * PB**T | * Call DORMBR to compute VB**T * PB**T | ||||
| * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | ||||
| @@ -145,15 +145,33 @@ | |||||
| INTRINSIC ABS, MAX | INTRINSIC ABS, MAX | ||||
| * .. | * .. | ||||
| * .. Executable Statements .. | * .. Executable Statements .. | ||||
| * | |||||
| INFO = 0 | |||||
| * | |||||
| * Quick return if possible | |||||
| * | |||||
| IF( N.EQ.0 ) | |||||
| $ RETURN | |||||
| * | * | ||||
| * Set constants to control overflow | * Set constants to control overflow | ||||
| * | * | ||||
| INFO = 0 | |||||
| EPS = DLAMCH( 'P' ) | EPS = DLAMCH( 'P' ) | ||||
| SMLNUM = DLAMCH( 'S' ) / EPS | SMLNUM = DLAMCH( 'S' ) / EPS | ||||
| BIGNUM = ONE / SMLNUM | BIGNUM = ONE / SMLNUM | ||||
| CALL DLABAD( SMLNUM, BIGNUM ) | CALL DLABAD( SMLNUM, BIGNUM ) | ||||
| * | * | ||||
| * Handle the case N=1 by itself | |||||
| * | |||||
| IF( N.EQ.1 ) THEN | |||||
| IPIV( 1 ) = 1 | |||||
| JPIV( 1 ) = 1 | |||||
| IF( ABS( A( 1, 1 ) ).LT.SMLNUM ) THEN | |||||
| INFO = 1 | |||||
| A( 1, 1 ) = SMLNUM | |||||
| END IF | |||||
| RETURN | |||||
| END IF | |||||
| * | |||||
| * Factorize A using complete pivoting. | * Factorize A using complete pivoting. | ||||
| * Set pivots less than SMIN to SMIN. | * Set pivots less than SMIN to SMIN. | ||||
| * | * | ||||
| @@ -418,9 +418,9 @@ | |||||
| $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | ||||
| END IF | END IF | ||||
| * | * | ||||
| * If INFO > 0 from SHSEQR, then quit | |||||
| * If INFO .NE. 0 from SHSEQR, then quit | |||||
| * | * | ||||
| IF( INFO.GT.0 ) | |||||
| IF( INFO.NE.0 ) | |||||
| $ GO TO 50 | $ GO TO 50 | ||||
| * | * | ||||
| IF( WANTVL .OR. WANTVR ) THEN | IF( WANTVL .OR. WANTVR ) THEN | ||||
| @@ -169,7 +169,7 @@ | |||||
| *> vectors, stored columnwise) as specified by RANGE; if | *> vectors, stored columnwise) as specified by RANGE; if | ||||
| *> JOBU = 'N', U is not referenced. | *> JOBU = 'N', U is not referenced. | ||||
| *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | ||||
| *> the exact value of NS is not known ILQFin advance and an upper | |||||
| *> the exact value of NS is not known in advance and an upper | |||||
| *> bound must be used. | *> bound must be used. | ||||
| *> \endverbatim | *> \endverbatim | ||||
| *> | *> | ||||
| @@ -357,8 +357,14 @@ | |||||
| IF( INFO.EQ.0 ) THEN | IF( INFO.EQ.0 ) THEN | ||||
| IF( WANTU .AND. LDU.LT.M ) THEN | IF( WANTU .AND. LDU.LT.M ) THEN | ||||
| INFO = -15 | INFO = -15 | ||||
| ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN | |||||
| INFO = -16 | |||||
| ELSE IF( WANTVT ) THEN | |||||
| IF( INDS ) THEN | |||||
| IF( LDVT.LT.IU-IL+1 ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| ELSE IF( LDVT.LT.MINMN ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -380,18 +386,34 @@ | |||||
| * | * | ||||
| * Path 1 (M much larger than N) | * Path 1 (M much larger than N) | ||||
| * | * | ||||
| MAXWRK = N*(N*2+16) + | |||||
| MAXWRK = N + | |||||
| $ N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) | $ N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) | ||||
| MAXWRK = MAX( MAXWRK, N*(N*2+20) + 2*N* | |||||
| MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N* | |||||
| $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) | $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) | ||||
| MINWRK = N*(N*2+21) | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*3+6)+N* | |||||
| $ ILAENV( 1, 'SORMQR', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*3+6)+N* | |||||
| $ ILAENV( 1, 'SORMLQ', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = N*(N*3+20) | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2 (M at least N, but not much larger) | * Path 2 (M at least N, but not much larger) | ||||
| * | * | ||||
| MAXWRK = N*(N*2+19) + ( M+N )* | |||||
| MAXWRK = 4*N + ( M+N )* | |||||
| $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) | $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) | ||||
| MINWRK = N*(N*2+20) + M | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*2+5)+N* | |||||
| $ ILAENV( 1, 'SORMQR', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,N*(N*2+5)+N* | |||||
| $ ILAENV( 1, 'SORMLQ', ' ', N, N, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = MAX(N*(N*2+19),4*N+M) | |||||
| END IF | END IF | ||||
| ELSE | ELSE | ||||
| MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) | MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) | ||||
| @@ -399,18 +421,34 @@ | |||||
| * | * | ||||
| * Path 1t (N much larger than M) | * Path 1t (N much larger than M) | ||||
| * | * | ||||
| MAXWRK = M*(M*2+16) + | |||||
| MAXWRK = M + | |||||
| $ M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) | $ M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) | ||||
| MAXWRK = MAX( MAXWRK, M*(M*2+20) + 2*M* | |||||
| MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M* | |||||
| $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) | $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) | ||||
| MINWRK = M*(M*2+21) | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*3+6)+M* | |||||
| $ ILAENV( 1, 'SORMQR', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*3+6)+M* | |||||
| $ ILAENV( 1, 'SORMLQ', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = M*(M*3+20) | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2t (N greater than M, but not much larger) | |||||
| * Path 2t (N at least M, but not much larger) | |||||
| * | * | ||||
| MAXWRK = M*(M*2+19) + ( M+N )* | |||||
| MAXWRK = 4*M + ( M+N )* | |||||
| $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) | $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) | ||||
| MINWRK = M*(M*2+20) + N | |||||
| IF (WANTU) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*2+5)+M* | |||||
| $ ILAENV( 1, 'SORMQR', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| IF (WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK,M*(M*2+5)+M* | |||||
| $ ILAENV( 1, 'SORMLQ', ' ', M, M, -1, -1 ) ) | |||||
| END IF | |||||
| MINWRK = MAX(M*(M*2+19),4*M+N) | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -522,7 +560,7 @@ | |||||
| CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | ||||
| J = J + N*2 | J = J + N*2 | ||||
| END DO | END DO | ||||
| CALL SLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| CALL SLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| * | * | ||||
| * Call SORMBR to compute QB*UB. | * Call SORMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -591,7 +629,7 @@ | |||||
| CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) | ||||
| J = J + N*2 | J = J + N*2 | ||||
| END DO | END DO | ||||
| CALL SLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| CALL SLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) | |||||
| * | * | ||||
| * Call SORMBR to compute QB*UB. | * Call SORMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -687,7 +725,7 @@ | |||||
| CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | ||||
| J = J + M*2 | J = J + M*2 | ||||
| END DO | END DO | ||||
| CALL SLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT ) | |||||
| CALL SLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) | |||||
| * | * | ||||
| * Call SORMBR to compute (VB**T)*(PB**T) | * Call SORMBR to compute (VB**T)*(PB**T) | ||||
| * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | ||||
| @@ -756,7 +794,7 @@ | |||||
| CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) | ||||
| J = J + M*2 | J = J + M*2 | ||||
| END DO | END DO | ||||
| CALL SLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT ) | |||||
| CALL SLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) | |||||
| * | * | ||||
| * Call SORMBR to compute VB**T * PB**T | * Call SORMBR to compute VB**T * PB**T | ||||
| * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | * (Workspace in WORK( ITEMP ): need M, prefer M*NB) | ||||
| @@ -145,15 +145,33 @@ | |||||
| INTRINSIC ABS, MAX | INTRINSIC ABS, MAX | ||||
| * .. | * .. | ||||
| * .. Executable Statements .. | * .. Executable Statements .. | ||||
| * | |||||
| INFO = 0 | |||||
| * | |||||
| * Quick return if possible | |||||
| * | |||||
| IF( N.EQ.0 ) | |||||
| $ RETURN | |||||
| * | * | ||||
| * Set constants to control overflow | * Set constants to control overflow | ||||
| * | * | ||||
| INFO = 0 | |||||
| EPS = SLAMCH( 'P' ) | EPS = SLAMCH( 'P' ) | ||||
| SMLNUM = SLAMCH( 'S' ) / EPS | SMLNUM = SLAMCH( 'S' ) / EPS | ||||
| BIGNUM = ONE / SMLNUM | BIGNUM = ONE / SMLNUM | ||||
| CALL SLABAD( SMLNUM, BIGNUM ) | CALL SLABAD( SMLNUM, BIGNUM ) | ||||
| * | * | ||||
| * Handle the case N=1 by itself | |||||
| * | |||||
| IF( N.EQ.1 ) THEN | |||||
| IPIV( 1 ) = 1 | |||||
| JPIV( 1 ) = 1 | |||||
| IF( ABS( A( 1, 1 ) ).LT.SMLNUM ) THEN | |||||
| INFO = 1 | |||||
| A( 1, 1 ) = SMLNUM | |||||
| END IF | |||||
| RETURN | |||||
| END IF | |||||
| * | |||||
| * Factorize A using complete pivoting. | * Factorize A using complete pivoting. | ||||
| * Set pivots less than SMIN to SMIN. | * Set pivots less than SMIN to SMIN. | ||||
| * | * | ||||
| @@ -404,9 +404,9 @@ | |||||
| $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | $ WORK( IWRK ), LWORK-IWRK+1, INFO ) | ||||
| END IF | END IF | ||||
| * | * | ||||
| * If INFO > 0 from ZHSEQR, then quit | |||||
| * If INFO .NE. 0 from ZHSEQR, then quit | |||||
| * | * | ||||
| IF( INFO.GT.0 ) | |||||
| IF( INFO.NE.0 ) | |||||
| $ GO TO 50 | $ GO TO 50 | ||||
| * | * | ||||
| IF( WANTVL .OR. WANTVR ) THEN | IF( WANTVL .OR. WANTVR ) THEN | ||||
| @@ -36,27 +36,30 @@ | |||||
| * .. | * .. | ||||
| * | * | ||||
| * | * | ||||
| * Purpose | |||||
| * ======= | |||||
| * | |||||
| * ZGESVDX computes the singular value decomposition (SVD) of a complex | |||||
| * M-by-N matrix A, optionally computing the left and/or right singular | |||||
| * vectors. The SVD is written | |||||
| * | |||||
| * A = U * SIGMA * transpose(V) | |||||
| * | |||||
| * where SIGMA is an M-by-N matrix which is zero except for its | |||||
| * min(m,n) diagonal elements, U is an M-by-M unitary matrix, and | |||||
| * V is an N-by-N unitary matrix. The diagonal elements of SIGMA | |||||
| * are the singular values of A; they are real and non-negative, and | |||||
| * are returned in descending order. The first min(m,n) columns of | |||||
| * U and V are the left and right singular vectors of A. | |||||
| * | |||||
| * ZGESVDX uses an eigenvalue problem for obtaining the SVD, which | |||||
| * allows for the computation of a subset of singular values and | |||||
| * vectors. See DBDSVDX for details. | |||||
| * | |||||
| * Note that the routine returns V**T, not V. | |||||
| *> \par Purpose: | |||||
| * ============= | |||||
| *> | |||||
| *> \verbatim | |||||
| *> | |||||
| *> ZGESVDX computes the singular value decomposition (SVD) of a complex | |||||
| *> M-by-N matrix A, optionally computing the left and/or right singular | |||||
| *> vectors. The SVD is written | |||||
| *> | |||||
| *> A = U * SIGMA * transpose(V) | |||||
| *> | |||||
| *> where SIGMA is an M-by-N matrix which is zero except for its | |||||
| *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and | |||||
| *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA | |||||
| *> are the singular values of A; they are real and non-negative, and | |||||
| *> are returned in descending order. The first min(m,n) columns of | |||||
| *> U and V are the left and right singular vectors of A. | |||||
| *> | |||||
| *> ZGESVDX uses an eigenvalue problem for obtaining the SVD, which | |||||
| *> allows for the computation of a subset of singular values and | |||||
| *> vectors. See DBDSVDX for details. | |||||
| *> | |||||
| *> Note that the routine returns V**T, not V. | |||||
| *> \endverbatim | |||||
| * | * | ||||
| * Arguments: | * Arguments: | ||||
| * ========== | * ========== | ||||
| @@ -107,7 +110,7 @@ | |||||
| *> | *> | ||||
| *> \param[in,out] A | *> \param[in,out] A | ||||
| *> \verbatim | *> \verbatim | ||||
| *> A is COMPLEX array, dimension (LDA,N) | |||||
| *> A is COMPLEX*16 array, dimension (LDA,N) | |||||
| *> On entry, the M-by-N matrix A. | *> On entry, the M-by-N matrix A. | ||||
| *> On exit, the contents of A are destroyed. | *> On exit, the contents of A are destroyed. | ||||
| *> \endverbatim | *> \endverbatim | ||||
| @@ -167,7 +170,7 @@ | |||||
| *> vectors, stored columnwise) as specified by RANGE; if | *> vectors, stored columnwise) as specified by RANGE; if | ||||
| *> JOBU = 'N', U is not referenced. | *> JOBU = 'N', U is not referenced. | ||||
| *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', | ||||
| *> the exact value of NS is not known ILQFin advance and an upper | |||||
| *> the exact value of NS is not known in advance and an upper | |||||
| *> bound must be used. | *> bound must be used. | ||||
| *> \endverbatim | *> \endverbatim | ||||
| *> | *> | ||||
| @@ -291,8 +294,8 @@ | |||||
| CHARACTER JOBZ, RNGTGK | CHARACTER JOBZ, RNGTGK | ||||
| LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT | LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT | ||||
| INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, | INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL, | ||||
| $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK, | |||||
| $ J, K, MAXWRK, MINMN, MINWRK, MNTHR | |||||
| $ ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ, | |||||
| $ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR | |||||
| DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM | DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM | ||||
| * .. | * .. | ||||
| * .. Local Arrays .. | * .. Local Arrays .. | ||||
| @@ -364,8 +367,14 @@ | |||||
| IF( INFO.EQ.0 ) THEN | IF( INFO.EQ.0 ) THEN | ||||
| IF( WANTU .AND. LDU.LT.M ) THEN | IF( WANTU .AND. LDU.LT.M ) THEN | ||||
| INFO = -15 | INFO = -15 | ||||
| ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN | |||||
| INFO = -16 | |||||
| ELSE IF( WANTVT ) THEN | |||||
| IF( INDS ) THEN | |||||
| IF( LDVT.LT.IU-IL+1 ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| ELSE IF( LDVT.LT.MINMN ) THEN | |||||
| INFO = -17 | |||||
| END IF | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -387,18 +396,24 @@ | |||||
| * | * | ||||
| * Path 1 (M much larger than N) | * Path 1 (M much larger than N) | ||||
| * | * | ||||
| MAXWRK = N + N* | |||||
| $ ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) | |||||
| MAXWRK = MAX( MAXWRK, N*N + N + 2*N* | |||||
| $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) | |||||
| MINWRK = N*(N+4) | |||||
| MINWRK = N*(N+5) | |||||
| MAXWRK = N + N*ILAENV(1,'ZGEQRF',' ',M,N,-1,-1) | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ N*N+2*N+2*N*ILAENV(1,'ZGEBRD',' ',N,N,-1,-1)) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ N*N+2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1)) | |||||
| END IF | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2 (M at least N, but not much larger) | * Path 2 (M at least N, but not much larger) | ||||
| * | * | ||||
| MAXWRK = 2*N + ( M+N )* | |||||
| $ ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 ) | |||||
| MINWRK = 2*N + M | |||||
| MINWRK = 3*N + M | |||||
| MAXWRK = 2*N + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ 2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1)) | |||||
| END IF | |||||
| END IF | END IF | ||||
| ELSE | ELSE | ||||
| MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) | MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) | ||||
| @@ -406,18 +421,25 @@ | |||||
| * | * | ||||
| * Path 1t (N much larger than M) | * Path 1t (N much larger than M) | ||||
| * | * | ||||
| MAXWRK = M + M* | |||||
| $ ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) | |||||
| MAXWRK = MAX( MAXWRK, M*M + M + 2*M* | |||||
| $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) | |||||
| MINWRK = M*(M+4) | |||||
| MINWRK = M*(M+5) | |||||
| MAXWRK = M + M*ILAENV(1,'ZGELQF',' ',M,N,-1,-1) | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ M*M+2*M+2*M*ILAENV(1,'ZGEBRD',' ',M,M,-1,-1)) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ M*M+2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1)) | |||||
| END IF | |||||
| ELSE | ELSE | ||||
| * | * | ||||
| * Path 2t (N greater than M, but not much larger) | * Path 2t (N greater than M, but not much larger) | ||||
| * | * | ||||
| MAXWRK = M*(M*2+19) + ( M+N )* | |||||
| $ ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 ) | |||||
| MINWRK = 2*M + N | |||||
| * | |||||
| MINWRK = 3*M + N | |||||
| MAXWRK = 2*M + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1) | |||||
| IF (WANTU .OR. WANTVT) THEN | |||||
| MAXWRK = MAX(MAXWRK, | |||||
| $ 2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1)) | |||||
| END IF | |||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| END IF | END IF | ||||
| @@ -515,14 +537,14 @@ | |||||
| CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), | CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), | ||||
| $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | ||||
| $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + N*(N*2+1) | |||||
| ITEMPR = ITGKZ + N*(N*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*N*N+14*N) | * (Workspace: need 2*N*N+14*N) | ||||
| * | * | ||||
| CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -536,7 +558,7 @@ | |||||
| END DO | END DO | ||||
| K = K + N | K = K + N | ||||
| END DO | END DO | ||||
| CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU ) | |||||
| CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) | |||||
| * | * | ||||
| * Call ZUNMBR to compute QB*UB. | * Call ZUNMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -591,14 +613,14 @@ | |||||
| CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | ||||
| $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | ||||
| $ LWORK-ITEMP+1, INFO ) | $ LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + N*(N*2+1) | |||||
| ITEMPR = ITGKZ + N*(N*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*N*N+14*N) | * (Workspace: need 2*N*N+14*N) | ||||
| * | * | ||||
| CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -612,7 +634,7 @@ | |||||
| END DO | END DO | ||||
| K = K + N | K = K + N | ||||
| END DO | END DO | ||||
| CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU ) | |||||
| CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) | |||||
| * | * | ||||
| * Call ZUNMBR to compute QB*UB. | * Call ZUNMBR to compute QB*UB. | ||||
| * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | * (Workspace in WORK( ITEMP ): need N, prefer N*NB) | ||||
| @@ -678,14 +700,14 @@ | |||||
| CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), | CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), | ||||
| $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | $ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), | ||||
| $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | $ WORK( ITEMP ), LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + M*(M*2+1) | |||||
| ITEMPR = ITGKZ + M*(M*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*M*M+14*M) | * (Workspace: need 2*M*M+14*M) | ||||
| * | * | ||||
| CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), | CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -719,7 +741,7 @@ | |||||
| END DO | END DO | ||||
| K = K + M | K = K + M | ||||
| END DO | END DO | ||||
| CALL ZLASET( 'A', M, N-M, CZERO, CZERO, | |||||
| CALL ZLASET( 'A', NS, N-M, CZERO, CZERO, | |||||
| $ VT( 1,M+1 ), LDVT ) | $ VT( 1,M+1 ), LDVT ) | ||||
| * | * | ||||
| * Call ZUNMBR to compute (VB**T)*(PB**T) | * Call ZUNMBR to compute (VB**T)*(PB**T) | ||||
| @@ -755,14 +777,14 @@ | |||||
| CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), | ||||
| $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ), | ||||
| $ LWORK-ITEMP+1, INFO ) | $ LWORK-ITEMP+1, INFO ) | ||||
| ITEMP = ITGKZ + M*(M*2+1) | |||||
| ITEMPR = ITGKZ + M*(M*2+1) | |||||
| * | * | ||||
| * Solve eigenvalue problem TGK*Z=Z*S. | * Solve eigenvalue problem TGK*Z=Z*S. | ||||
| * (Workspace: need 2*M*M+14*M) | * (Workspace: need 2*M*M+14*M) | ||||
| * | * | ||||
| CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), | CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), | ||||
| $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | $ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, | ||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMP ), | |||||
| $ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), | |||||
| $ IWORK, INFO) | $ IWORK, INFO) | ||||
| * | * | ||||
| * If needed, compute left singular vectors. | * If needed, compute left singular vectors. | ||||
| @@ -796,7 +818,7 @@ | |||||
| END DO | END DO | ||||
| K = K + M | K = K + M | ||||
| END DO | END DO | ||||
| CALL ZLASET( 'A', M, N-M, CZERO, CZERO, | |||||
| CALL ZLASET( 'A', NS, N-M, CZERO, CZERO, | |||||
| $ VT( 1,M+1 ), LDVT ) | $ VT( 1,M+1 ), LDVT ) | ||||
| * | * | ||||
| * Call ZUNMBR to compute VB**T * PB**T | * Call ZUNMBR to compute VB**T * PB**T | ||||
| @@ -145,15 +145,33 @@ | |||||
| INTRINSIC ABS, DCMPLX, MAX | INTRINSIC ABS, DCMPLX, MAX | ||||
| * .. | * .. | ||||
| * .. Executable Statements .. | * .. Executable Statements .. | ||||
| * | |||||
| INFO = 0 | |||||
| * | |||||
| * Quick return if possible | |||||
| * | |||||
| IF( N.EQ.0 ) | |||||
| $ RETURN | |||||
| * | * | ||||
| * Set constants to control overflow | * Set constants to control overflow | ||||
| * | * | ||||
| INFO = 0 | |||||
| EPS = DLAMCH( 'P' ) | EPS = DLAMCH( 'P' ) | ||||
| SMLNUM = DLAMCH( 'S' ) / EPS | SMLNUM = DLAMCH( 'S' ) / EPS | ||||
| BIGNUM = ONE / SMLNUM | BIGNUM = ONE / SMLNUM | ||||
| CALL DLABAD( SMLNUM, BIGNUM ) | CALL DLABAD( SMLNUM, BIGNUM ) | ||||
| * | * | ||||
| * Handle the case N=1 by itself | |||||
| * | |||||
| IF( N.EQ.1 ) THEN | |||||
| IPIV( 1 ) = 1 | |||||
| JPIV( 1 ) = 1 | |||||
| IF( ABS( A( 1, 1 ) ).LT.SMLNUM ) THEN | |||||
| INFO = 1 | |||||
| A( 1, 1 ) = DCMPLX( SMLNUM, ZERO ) | |||||
| END IF | |||||
| RETURN | |||||
| END IF | |||||
| * | |||||
| * Factorize A using complete pivoting. | * Factorize A using complete pivoting. | ||||
| * Set pivots less than SMIN to SMIN | * Set pivots less than SMIN to SMIN | ||||
| * | * | ||||
| @@ -340,7 +340,7 @@ | |||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| CALL ZHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | CALL ZHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | ||||
| $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1, | $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1, | ||||
| $ WORK, IERR ) | |||||
| $ RWORK, IERR ) | |||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| ELSE | ELSE | ||||
| CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL, | CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL, | ||||
| @@ -348,7 +348,7 @@ | |||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| CALL ZHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | CALL ZHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, | ||||
| $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1, | $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1, | ||||
| $ WORK, IERR ) | |||||
| $ RWORK, IERR ) | |||||
| LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) ) | ||||
| END IF | END IF | ||||
| WORK( 1 ) = DCMPLX( LWKOPT ) | WORK( 1 ) = DCMPLX( LWKOPT ) | ||||