bugfixes form lapack svn for bugs #142 - #155
authorWerner Saar <wernsaar@googlemail.com>
Mon, 7 Mar 2016 09:10:00 +0000 (10:10 +0100)
committerWerner Saar <wernsaar@googlemail.com>
Mon, 7 Mar 2016 09:10:00 +0000 (10:10 +0100)
25 files changed:
lapack-netlib/LAPACKE/src/lapacke_clantr.c
lapack-netlib/LAPACKE/src/lapacke_dlantr.c
lapack-netlib/LAPACKE/src/lapacke_dlantr_work.c
lapack-netlib/LAPACKE/src/lapacke_dormbr_work.c
lapack-netlib/LAPACKE/src/lapacke_dormlq_work.c
lapack-netlib/LAPACKE/src/lapacke_slantr.c
lapack-netlib/LAPACKE/src/lapacke_slantr_work.c
lapack-netlib/LAPACKE/src/lapacke_sormbr_work.c
lapack-netlib/LAPACKE/src/lapacke_sormlq_work.c
lapack-netlib/LAPACKE/src/lapacke_zlantr.c
lapack-netlib/LAPACKE/src/lapacke_zlantr_work.c
lapack-netlib/SRC/cgeev.f
lapack-netlib/SRC/cgesvdx.f
lapack-netlib/SRC/cgetc2.f
lapack-netlib/SRC/cggev3.f
lapack-netlib/SRC/dgeev.f
lapack-netlib/SRC/dgesvdx.f
lapack-netlib/SRC/dgetc2.f
lapack-netlib/SRC/sgeev.f
lapack-netlib/SRC/sgesvdx.f
lapack-netlib/SRC/sgetc2.f
lapack-netlib/SRC/zgeev.f
lapack-netlib/SRC/zgesvdx.f
lapack-netlib/SRC/zgetc2.f
lapack-netlib/SRC/zggev3.f

index 4ba4367..33e6e57 100644 (file)
@@ -51,8 +51,7 @@ float LAPACKE_clantr( int matrix_layout, char norm, char uplo, char diag,
     }
 #endif
     /* 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)) );
         if( work == NULL ) {
             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,
                                 work );
     /* 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 );
     }
 exit_level_0:
index 65802f3..8fd1120 100644 (file)
@@ -51,8 +51,7 @@ double LAPACKE_dlantr( int matrix_layout, char norm, char uplo, char diag,
     }
 #endif
     /* 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)) );
         if( work == NULL ) {
             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,
                                 work );
     /* 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 );
     }
 exit_level_0:
index 59eef38..ee5e665 100644 (file)
@@ -38,10 +38,10 @@ double LAPACKE_dlantr_work( int matrix_layout, char norm, char uplo,
                                 const double* a, lapack_int lda, double* work )
 {
     lapack_int info = 0;
-       double res = 0.;
+    double res = 0.;
     if( matrix_layout == LAPACK_COL_MAJOR ) {
         /* 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 ) {
             info = info - 1;
         }
index 9db92ce..dcd8842 100644 (file)
@@ -74,11 +74,10 @@ lapack_int LAPACKE_dormbr_work( int matrix_layout, char vect, char side,
         }
         /* Allocate memory for temporary array(s) */
         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 {
-          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 ) {
             info = LAPACK_TRANSPOSE_MEMORY_ERROR;
             goto exit_level_0;
@@ -89,11 +88,7 @@ lapack_int LAPACKE_dormbr_work( int matrix_layout, char vect, char side,
             goto exit_level_1;
         }
         /* 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 );
         /* Call LAPACK function and adjust info */
         LAPACK_dormbr( &vect, &side, &trans, &m, &n, &k, a_t, &lda_t, tau, c_t,
index 2a59cd5..f46c6d3 100644 (file)
@@ -87,12 +87,7 @@ lapack_int LAPACKE_dormlq_work( int matrix_layout, char side, char trans,
             goto exit_level_1;
         }
         /* 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 );
         /* Call LAPACK function and adjust info */
         LAPACK_dormlq( &side, &trans, &m, &n, &k, a_t, &lda_t, tau, c_t, &ldc_t,
index 99942b3..23b19b1 100644 (file)
@@ -51,8 +51,7 @@ float LAPACKE_slantr( int matrix_layout, char norm, char uplo, char diag,
     }
 #endif
     /* 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)) );
         if( work == NULL ) {
             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,
                                 work );
     /* 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 );
     }
 exit_level_0:
index 79c71a0..92a0e40 100644 (file)
@@ -41,7 +41,7 @@ float LAPACKE_slantr_work( int matrix_layout, char norm, char uplo,
     float res = 0.;
     if( matrix_layout == LAPACK_COL_MAJOR ) {
         /* 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 ) {
             info = info - 1;
         }
index f7a2ff3..34bc41f 100644 (file)
@@ -73,8 +73,11 @@ lapack_int LAPACKE_sormbr_work( int matrix_layout, char vect, char side,
             return (info < 0) ? (info - 1) : info;
         }
         /* 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 ) {
             info = LAPACK_TRANSPOSE_MEMORY_ERROR;
             goto exit_level_0;
index a277436..b02a2d1 100644 (file)
@@ -72,7 +72,11 @@ lapack_int LAPACKE_sormlq_work( int matrix_layout, char side, char trans,
             return (info < 0) ? (info - 1) : info;
         }
         /* 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 ) {
             info = LAPACK_TRANSPOSE_MEMORY_ERROR;
             goto exit_level_0;
index b2c6371..2b645e7 100644 (file)
@@ -51,8 +51,7 @@ double LAPACKE_zlantr( int matrix_layout, char norm, char uplo, char diag,
     }
 #endif
     /* 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)) );
         if( work == NULL ) {
             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,
                                 work );
     /* 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 );
     }
 exit_level_0:
index 5fd9f84..0988bf6 100644 (file)
@@ -39,7 +39,7 @@ double LAPACKE_zlantr_work( int matrix_layout, char norm, char uplo,
                                 double* work )
 {
     lapack_int info = 0;
-       double res = 0.;
+    double res = 0.;
     if( matrix_layout == LAPACK_COL_MAJOR ) {
         /* Call LAPACK function and adjust info */
         res = LAPACK_zlantr( &norm, &uplo, &diag, &m, &n, a, &lda, work );
index b79b645..0f48322 100644 (file)
      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
       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
 *
       IF( WANTVL .OR. WANTVR ) THEN
index 235426a..87ea986 100644 (file)
 *>          vectors, stored columnwise) as specified by RANGE; if 
 *>          JOBU = 'N', U is not referenced.
 *>          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.
 *> \endverbatim
 *>
       CHARACTER          JOBZ, RNGTGK
       LOGICAL            ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
       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
 *     ..
 *     .. Local Arrays ..
          IF( INFO.EQ.0 ) THEN
             IF( WANTU .AND. LDU.LT.M ) THEN
                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
 *
 *                 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
 *
 *                 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
             ELSE
                MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *
 *                 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
 *
 *                 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
             CALL CGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), 
      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
      $                   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.
 *           (Workspace: need 2*N*N+14*N)          
 *                   
             CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
-     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMP ),
+     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)
 *
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + N
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
             CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), 
      $                   WORK( ITAUQ ), WORK( ITAUP ), 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.
 *           (Workspace: need 2*N*N+14*N)          
 *                     
             CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, 
-     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMP ), 
+     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), 
      $                    IWORK, INFO)
 *
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + N
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
             CALL CGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ),
      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), 
      $                   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.
 *           (Workspace: need 2*M*M+14*M)          
 *
             CALL SBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, 
-     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMP ), 
+     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), 
      $                    IWORK, INFO)
 *
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + M
                END DO
-               CALL CLASET( 'A', M, N-M, CZERO, CZERO, 
+               CALL CLASET( 'A', NS, N-M, CZERO, CZERO,
      $                      VT( 1,M+1 ), LDVT )
 *
 *              Call CUNMBR to compute (VB**T)*(PB**T)
             CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), 
      $                   WORK( ITAUQ ), WORK( ITAUP ), 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.
 *           (Workspace: need 2*M*M+14*M)          
 *          
             CALL SBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), 
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, 
-     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMP ), 
+     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), 
      $                    IWORK, INFO)
 * 
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + M
                END DO
-               CALL CLASET( 'A', M, N-M, CZERO, CZERO, 
+               CALL CLASET( 'A', NS, N-M, CZERO, CZERO,
      $                      VT( 1,M+1 ), LDVT )
 *
 *              Call CUNMBR to compute VB**T * PB**T
index fac6b56..99eb69d 100644 (file)
 *     ..
 *     .. Executable Statements ..
 *
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
 *     Set constants to control overflow
 *
-      INFO = 0
       EPS = SLAMCH( 'P' )
       SMLNUM = SLAMCH( 'S' ) / EPS
       BIGNUM = ONE / SMLNUM
       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.
 *     Set pivots less than SMIN to SMIN
 *
index 4a000fe..decdae5 100644 (file)
      $                   LDVL, VR, LDVR, WORK, -1, IERR )
             LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
             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 ) ) )
          ELSE
             CALL CGGHD3( 'N', 'N', N, 1, N, A, LDA, B, LDB, VL, LDVL,
      $                   VR, LDVR, WORK, -1, IERR )
             LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
             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 ) ) )
          END IF
          WORK( 1 ) = CMPLX( LWKOPT )
index dd60db6..328eaa3 100644 (file)
      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
       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
 *
       IF( WANTVL .OR. WANTVR ) THEN
index cfa2ff0..4588083 100644 (file)
 *>          vectors, stored columnwise) as specified by RANGE; if 
 *>          JOBU = 'N', U is not referenced.
 *>          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.
 *> \endverbatim
 *>
          IF( INFO.EQ.0 ) THEN
             IF( WANTU .AND. LDU.LT.M ) THEN
                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
 *
 *                 Path 1 (M much larger than N)
 *
-                  MAXWRK = N*(N*2+16) + 
+                  MAXWRK = N + 
      $                     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 ) )
-                  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
 *
 *                 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 )
-                  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
             ELSE
                MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *
 *                 Path 1t (N much larger than M)
 *
-                  MAXWRK = M*(M*2+16) + 
+                  MAXWRK = M + 
      $                     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 ) )
-                  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
 *
-*                 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 )
-                  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
                   CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
                   J = J + N*2
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
                   CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
                   J = J + N*2
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
                   CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
                   J = J + M*2
                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)
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
                   CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
                   J = J + M*2
                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
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
index 7e43a02..3cd7eeb 100644 (file)
 *     ..
 *     .. Executable Statements ..
 *
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
 *     Set constants to control overflow
 *
-      INFO = 0
       EPS = DLAMCH( 'P' )
       SMLNUM = DLAMCH( 'S' ) / EPS
       BIGNUM = ONE / SMLNUM
       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.
 *     Set pivots less than SMIN to SMIN.
 *
index 89dbe08..667de0a 100644 (file)
      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
       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
 *
       IF( WANTVL .OR. WANTVR ) THEN
index aae8b07..9128a7c 100644 (file)
 *>          vectors, stored columnwise) as specified by RANGE; if 
 *>          JOBU = 'N', U is not referenced.
 *>          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.
 *> \endverbatim
 *>
          IF( INFO.EQ.0 ) THEN
             IF( WANTU .AND. LDU.LT.M ) THEN
                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
 *
 *                 Path 1 (M much larger than N)
 *
-                  MAXWRK = N*(N*2+16) + 
+                  MAXWRK = N + 
      $                     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 ) )
-                  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
 *
 *                 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 )
-                  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
             ELSE
                MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *
 *                 Path 1t (N much larger than M)
 *
-                  MAXWRK = M*(M*2+16) + 
+                  MAXWRK = M + 
      $                     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 ) )
-                  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
 *
-*                 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 )
-                  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
                   CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
                   J = J + N*2
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
                   CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
                   J = J + N*2
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
                   CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
                   J = J + M*2
                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)
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
                   CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
                   J = J + M*2
                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
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
index 3c3880d..5984465 100644 (file)
 *     ..
 *     .. Executable Statements ..
 *
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
 *     Set constants to control overflow
 *
-      INFO = 0
       EPS = SLAMCH( 'P' )
       SMLNUM = SLAMCH( 'S' ) / EPS
       BIGNUM = ONE / SMLNUM
       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.
 *     Set pivots less than SMIN to SMIN.
 *
index d452080..a518b4c 100644 (file)
      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
       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
 *
       IF( WANTVL .OR. WANTVR ) THEN
index 6f7d5ba..c9509e4 100644 (file)
 *     ..
 *
 *
-*  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:
 *  ==========
 *>
 *> \param[in,out] A
 *> \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 exit, the contents of A are destroyed.
 *> \endverbatim
 *>          vectors, stored columnwise) as specified by RANGE; if 
 *>          JOBU = 'N', U is not referenced.
 *>          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.
 *> \endverbatim
 *>
       CHARACTER          JOBZ, RNGTGK
       LOGICAL            ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
       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
 *     ..
 *     .. Local Arrays ..
          IF( INFO.EQ.0 ) THEN
             IF( WANTU .AND. LDU.LT.M ) THEN
                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
 *
 *                 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
 *
 *                 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
             ELSE
                MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
 *
 *                 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
 *
 *                 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
             CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), 
      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
      $                   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.
 *           (Workspace: need 2*N*N+14*N)          
 *                   
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
-     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMP ),
+     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)
 *
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + N
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
             CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), 
      $                   WORK( ITAUQ ), WORK( ITAUP ), 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.
 *           (Workspace: need 2*N*N+14*N)          
 *                     
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, 
-     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMP ), 
+     $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ), 
      $                    IWORK, INFO)
 *
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + N
                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.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
             CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ),
      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ), 
      $                   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.
 *           (Workspace: need 2*M*M+14*M)          
 *
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, 
-     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMP ), 
+     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), 
      $                    IWORK, INFO)
 *
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + M
                END DO
-               CALL ZLASET( 'A', M, N-M, CZERO, CZERO, 
+               CALL ZLASET( 'A', NS, N-M, CZERO, CZERO,
      $                      VT( 1,M+1 ), LDVT )
 *
 *              Call ZUNMBR to compute (VB**T)*(PB**T)
             CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), 
      $                   WORK( ITAUQ ), WORK( ITAUP ), 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.
 *           (Workspace: need 2*M*M+14*M)          
 *          
             CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), 
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S, 
-     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMP ), 
+     $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ), 
      $                    IWORK, INFO)
 * 
 *           If needed, compute left singular vectors.
                   END DO
                   K = K + M
                END DO
-               CALL ZLASET( 'A', M, N-M, CZERO, CZERO, 
+               CALL ZLASET( 'A', NS, N-M, CZERO, CZERO,
      $                      VT( 1,M+1 ), LDVT )
 *
 *              Call ZUNMBR to compute VB**T * PB**T
index 3179612..bf59415 100644 (file)
 *     ..
 *     .. Executable Statements ..
 *
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
 *     Set constants to control overflow
 *
-      INFO = 0
       EPS = DLAMCH( 'P' )
       SMLNUM = DLAMCH( 'S' ) / EPS
       BIGNUM = ONE / SMLNUM
       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.
 *     Set pivots less than SMIN to SMIN
 *
index 1c4e832..78337fd 100644 (file)
             LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
             CALL ZHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
      $                   ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1,
-     $                   WORK, IERR )
+     $                   RWORK, IERR )
             LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
          ELSE
             CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
             LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
             CALL ZHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
      $                   ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1,
-     $                   WORK, IERR )
+     $                   RWORK, IERR )
             LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
          END IF
          WORK( 1 ) = DCMPLX( LWKOPT )