While looking if bug0022 was corrected, found that some comments were not updated

bug0022 was indeed corrected by Zlatko

3 files changed, 80 insertions, 59 deletions

diff --git a/SRC/dgejsv.f b/SRC/dgejsv.f
index e61b27a9..fee34ddc 100644
--- a/SRC/dgejsv.f
+++ b/SRC/dgejsv.f
@@ -213,7 +213,7 @@
 *          If JOBV = 'V' or 'J' or 'W', then LDV >= N.
 *
 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension at least LWORK.
-*          On exit,
+*          On exit, if N.GT.0 .AND. M.GT.0 (else not referenced), 
 *          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
 *                    that SCALE*SVA(1:N) are the computed singular values
 *                    of A. (See the description of SVA().)
@@ -252,33 +252,50 @@
 *          LWORK depends on the job:
 *
 *          If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
-*            -> .. no scaled condition estimate required ( JOBE.EQ.'N'):
+*            -> .. no scaled condition estimate required (JOBE.EQ.'N'):
 *               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
-*               For optimal performance (blocked code) the optimal value
+*               ->> For optimal performance (blocked code) the optimal value
 *               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
-*               block size for xGEQP3/xGEQRF.
+*               block size for DGEQP3 and DGEQRF.
+*               In general, optimal LWORK is computed as 
+*               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).        
 *            -> .. an estimate of the scaled condition number of A is
 *               required (JOBA='E', 'G'). In this case, LWORK is the maximum
-*               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7).
+*               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7).
+*               ->> For optimal performance (blocked code) the optimal value 
+*               is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7).
+*               In general, the optimal length LWORK is computed as
+*               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 
+*                                                     N+N*N+LWORK(DPOCON),7).
 *
 *          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
-*            -> the minimal requirement is LWORK >= max(2*N+M,7).
-*            -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7),
-*               where NB is the optimal block size.
+*            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
+*            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
+*               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ,
+*               DORMLQ. In general, the optimal length LWORK is computed as
+*               LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), 
+*                       N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
 *
 *          If SIGMA and the left singular vectors are needed
-*            -> the minimal requirement is LWORK >= max(2*N+M,7).
-*            -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7),
-*               where NB is the optimal block size.
-*
-*          If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and
-*            -> .. the singular vectors are computed without explicit
-*               accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N
-*            -> .. in the iterative part, the Jacobi rotations are
-*               explicitly accumulated (option, see the description of JOBV),
-*               then the minimal requirement is LWORK >= max(M+3*N+N*N,7).
-*               For better performance, if NB is the optimal block size,
-*               LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7).
+*            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
+*            -> For optimal performance:
+*               if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
+*               if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
+*               where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
+*               In general, the optimal length LWORK is computed as
+*               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
+*                        2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). 
+*               Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or 
+*               M*NB (for JOBU.EQ.'F').
+*               
+*          If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and 
+*            -> if JOBV.EQ.'V'  
+*               the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). 
+*            -> if JOBV.EQ.'J' the minimal requirement is 
+*               LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
+*            -> For optimal performance, LWORK should be additionally
+*               larger than N+M*NB, where NB is the optimal block size
+*               for DORMQR.
 *
 *  IWORK   (workspace/output) INTEGER array, dimension M+3*N.
 *          On exit,
@@ -300,8 +317,8 @@
 *  Further Details
 *  ===============
 *
-*  DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3,
-*  SGEQRF, and SGELQF as preprocessors and preconditioners. Optionally, an
+*  DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses DGEQP3,
+*  DGEQRF, and DGELQF as preprocessors and preconditioners. Optionally, an
 *  additional row pivoting can be used as a preprocessor, which in some
 *  cases results in much higher accuracy. An example is matrix A with the
 *  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
@@ -319,14 +336,14 @@
 *  the singular values in the range of normalized IEEE numbers is that the
 *  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
 *  overflow. This code (DGEJSV) is best used in this restricted range,
-*  meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are
+*  meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
 *  returned as zeros. See JOBR for details on this.
 *     Further, this implementation is somewhat slower than the one described
 *  in [1,2] due to replacement of some non-LAPACK components, and because
 *  the choice of some tuning parameters in the iterative part (DGESVJ) is
 *  left to the implementer on a particular machine.
-*     The rank revealing QR factorization (in this code: SGEQP3) should be
-*  implemented as in [3]. We have a new version of SGEQP3 under development
+*     The rank revealing QR factorization (in this code: DGEQP3) should be
+*  implemented as in [3]. We have a new version of DGEQP3 under development
 *  that is more robust than the current one in LAPACK, with a cleaner cut in
 *  rank defficient cases. It will be available in the SIGMA library [4].
 *  If M is much larger than N, it is obvious that the inital QRF with
@@ -437,14 +454,18 @@
       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
          INFO = - 14
       ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
-     $                           (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR.
-     $ (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND.
-     $                         (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR.
-     $ (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR.
-     $ (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR.
-     $ (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N))
-     $ .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N)))
-     $   THEN
+     &                           (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR.
+     & (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND.
+     &                         (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR.
+     & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1)))
+     & .OR.
+     & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1)))
+     & .OR.
+     & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. 
+     &                          (LWORK.LT.MAX0(2*M+N,6*N+2*N*N)))
+     & .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND.
+     &                          LWORK.LT.MAX0(2*M+N,4*N+N*N,2*N+N*N+6)))
+     &   THEN
          INFO = - 17
       ELSE
 *        #:)
@@ -454,6 +475,7 @@
       IF ( INFO .NE. 0 ) THEN
 *       #:(
          CALL XERBLA( 'DGEJSV', - INFO )
+         RETURN
       END IF
 *
 *     Quick return for void matrix (Y3K safe)
@@ -471,7 +493,6 @@
 *
 *!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
 *
-
       EPSLN = DLAMCH('Epsilon')
       SFMIN = DLAMCH('SafeMinimum')
       SMALL = SFMIN / EPSLN
@@ -536,6 +557,7 @@
          END IF
          IWORK(1) = 0
          IWORK(2) = 0
+         IWORK(3) = 0
          RETURN
       END IF
 *
@@ -916,7 +938,6 @@
 *
 *     Phase 3:
 *
-
       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
 *
 *         Singular Values only
diff --git a/SRC/dgsvj1.f b/SRC/dgsvj1.f
index e9a3494e..ae4f8bfe 100644
--- a/SRC/dgsvj1.f
+++ b/SRC/dgsvj1.f
@@ -44,12 +44,12 @@
 *  block-entries (tiles) of the (1,2) off-diagonal block are marked by the
 *  [x]'s in the following scheme:
 *
-*     | *   C   * [x] [x] [x]|
-*     | *   C   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
-*     | *   C   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
+*     | *   *   * [x] [x] [x]|
+*     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
+*     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
 *
 *  In terms of the columns of A, the first N1 columns are rotated 'against'
 *  the remaining N-N1 columns, trying to increase the angle between the
@@ -271,12 +271,12 @@
 *     Jacobi SVD algorithm SGESVJ.
 *
 *
-*     | *   C   * [x] [x] [x]|
-*     | *   C   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
-*     | *   C   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
+*     | *   *   * [x] [x] [x]|
+*     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
+*     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
 *
 *
       DO 1993 i = 1, NSWEEP
diff --git a/SRC/sgsvj1.f b/SRC/sgsvj1.f
index bb36ea15..c1249b79 100644
--- a/SRC/sgsvj1.f
+++ b/SRC/sgsvj1.f
@@ -44,12 +44,12 @@
 *  block-entries (tiles) of the (1,2) off-diagonal block are marked by the
 *  [x]'s in the following scheme:
 *
-*     | *   C   * [x] [x] [x]|
-*     | *   C   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
-*     | *   C   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
+*     | *   *   * [x] [x] [x]|
+*     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
+*     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
 *
 *  In terms of the columns of A, the first N1 columns are rotated 'against'
 *  the remaining N-N1 columns, trying to increase the angle between the
@@ -271,12 +271,12 @@
 *     Jacobi SVD algorithm SGESVJ.
 *
 *
-*     | *   C   * [x] [x] [x]|
-*     | *   C   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
-*     | *   C   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
-*     |[x] [x] [x] *   C   * |
+*     | *   *   * [x] [x] [x]|
+*     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
+*     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
+*     |[x] [x] [x] *   *   * |
 *
 *
       DO 1993 i = 1, NSWEEP