Include Zlatko changes following Paul Roberts from NAG report

Message from Zlatko:

The following changes are made in the
current version of the code:

1. in dgejsv, sgejsv:

- a typo in checking the parameters LSVEC.OR.LSVEC
  has been changed to LSVEC.OR.RSVEC
- the length of WORK, LWORK, its description
  and the minimal length for different JOBs
  have been revised and corrected
- in a call to xGESVJ with WORK(N+1), the length
  of the workspace is set to the correct value of
  LWORK-N, instead of the incorrect LWORK
- a missing RETURN after a call to XERBLA
  has been inserted
- In the case of zero matrix on input, IWORK(3)
  is set to ZERO, to correspond to the description
  of IWORK in other nontrivial cases
2. in all routines
- simple editing so that the single and
  the corresponding double routines have matching
  lines in the source codes.

diff --git a/SRC/dgesvj.f b/SRC/dgesvj.f
index d48975e..a8e154f 100644 (file)
--- a/SRC/dgesvj.f
+++ b/SRC/dgesvj.f
@@ -1,11 +1,11 @@
       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
      +                   LDV, WORK, LWORK, INFO )
 *
-*  -- LAPACK routine (version 3.3.0)                                    --
+*  -- LAPACK routine (version 3.3.1)                                  --
 *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
-*     November 2010
+*     January 2011
 *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
@@ -133,7 +133,7 @@
 *                  referenced
 *
 *  M       (input) INTEGER
-*          The number of rows of the input matrix A.  M >= 0.
+*          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.  
 *
 *  N       (input) INTEGER
 *          The number of columns of the input matrix A.
@@ -383,7 +383,7 @@
       ROOTTOL = DSQRT( TOL )
 *
       IF( DBLE( M )*EPSLN.GE.ONE ) THEN
-         INFO = -5
+         INFO = -4
          CALL XERBLA( 'DGESVJ', -INFO )
          RETURN
       END IF
@@ -824,8 +824,7 @@
 *
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
-                                 THETA = -HALF*DABS( AQOAP-APOAQ ) /
-     +                                   AAPQ
+                                 THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
 *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
 *
@@ -1131,8 +1130,7 @@
 *
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
-                                 THETA = -HALF*DABS( AQOAP-APOAQ ) /
-     +                                   AAPQ
+                                 THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
 *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN

diff --git a/SRC/dgsvj0.f b/SRC/dgsvj0.f
index e03aeb6..c26f929 100644 (file)
--- a/SRC/dgsvj0.f
+++ b/SRC/dgsvj0.f
@@ -1,11 +1,11 @@
       SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
      +                   SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 *
-*  -- LAPACK routine (version 3.3.0)                                    --
+*  -- LAPACK routine (version 3.3.1)                                  --
 *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
-*     November 2010
+*     January 2011
 *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
@@ -16,6 +16,7 @@
 * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
 *
       IMPLICIT NONE
+*     ..
 *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
       DOUBLE PRECISION   EPS, SFMIN, TOL
@@ -173,6 +174,8 @@
 *     ..
 *     .. Executable Statements ..
 *
+*     Test the input parameters.
+*
       APPLV = LSAME( JOBV, 'A' )
       RSVEC = LSAME( JOBV, 'V' )
       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
@@ -219,7 +222,6 @@
       BIGTHETA = ONE / ROOTEPS
       ROOTTOL = DSQRT( TOL )
 *
-*
 *     -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#-
 *
       EMPTSW = ( N*( N-1 ) ) / 2
@@ -376,8 +378,7 @@
 *
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
-                                 THETA = -HALF*DABS( AQOAP-APOAQ ) /
-     +                                   AAPQ
+                                 THETA = -HALF*DABS( AQOAP-APOAQ )/AAPQ
 *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
 *
@@ -676,8 +677,7 @@
 *
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
-                                 THETA = -HALF*DABS( AQOAP-APOAQ ) /
-     +                                   AAPQ
+                                 THETA = -HALF*DABS( AQOAP-APOAQ )/AAPQ
                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
 *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN

diff --git a/SRC/dgsvj1.f b/SRC/dgsvj1.f
index 1a51f15..89704e7 100644 (file)
--- a/SRC/dgsvj1.f
+++ b/SRC/dgsvj1.f
@@ -1,11 +1,11 @@
       SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
      +                   EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 *
-*  -- LAPACK routine (version 3.3.0)                                    --
+*  -- LAPACK routine (version 3.3.1)                                  --
 *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
-*     November 2010
+*     January 2011
 *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
@@ -375,8 +375,7 @@
 *
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
-                                 THETA = -HALF*DABS( AQOAP-APOAQ ) /
-     +                                   AAPQ
+                                 THETA = -HALF*DABS(AQOAP-APOAQ) / AAPQ
                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
 
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN

diff --git a/SRC/sgejsv.f b/SRC/sgejsv.f
index e2138cf..aaef93f 100644 (file)
--- a/SRC/sgejsv.f
+++ b/SRC/sgejsv.f
@@ -2,11 +2,11 @@
      &                   M, N, A, LDA, SVA, U, LDU, V, LDV,
      &                   WORK, LWORK, IWORK, INFO )
 *
-*  -- LAPACK routine (version 3.3.0)                                    --
+*  -- LAPACK routine (version 3.3.1)                                  --
 *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
-*     November 2010
+*     January 2011
 *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
@@ -31,6 +31,7 @@
 *
 *  Purpose
 *  =======
+*
 *  SGEJSV computes the singular value decomposition (SVD) of a real M-by-N
 *  matrix [A], where M >= N. The SVD of [A] is written as
 *
@@ -253,33 +254,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,
@@ -439,12 +457,16 @@
          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.
+     & (.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*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)))
+     & (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
@@ -455,6 +477,7 @@
       IF ( INFO .NE. 0 ) THEN
 *       #:(
          CALL XERBLA( 'SGEJSV', - INFO )
+         RETURN
       END IF
 *
 *     Quick return for void matrix (Y3K safe)
@@ -476,6 +499,7 @@
       SFMIN = SLAMCH('SafeMinimum')
       SMALL = SFMIN / EPSLN
       BIG   = SLAMCH('O')
+*     BIG   = ONE / SFMIN
 *
 *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
 *
@@ -535,6 +559,7 @@
          END IF
          IWORK(1) = 0
          IWORK(2) = 0
+         IWORK(3) = 0
          RETURN
       END IF
 *
@@ -1028,7 +1053,7 @@
             CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
 *
             CALL SGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
-     &                  LDU, WORK(N+1), LWORK, INFO )
+     &                  LDU, WORK(N+1), LWORK-N, INFO )
             SCALEM  = WORK(N+1)
             NUMRANK = NINT(WORK(N+2))
             IF ( NR .LT. N ) THEN

diff --git a/SRC/sgesvj.f b/SRC/sgesvj.f
index 01a3d24..f5fa397 100644 (file)
--- a/SRC/sgesvj.f
+++ b/SRC/sgesvj.f
@@ -1,11 +1,11 @@
       SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
      +                   LDV, WORK, LWORK, INFO )
 *
-*  -- LAPACK routine (version 3.3.0)                                    --
+*  -- LAPACK routine (version 3.3.1)                                  --
 *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
-*     November 2010
+*     January 2011
 *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
@@ -133,7 +133,7 @@
 *                  referenced
 *
 *  M       (input) INTEGER
-*          The number of rows of the input matrix A.  M >= 0.
+*          The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.
 *
 *  N       (input) INTEGER
 *          The number of columns of the input matrix A.
@@ -241,7 +241,8 @@
 *                    Jacobi rotation angles in the last sweep. It can be
 *                    useful for a post festum analysis.
 *
-*  LWORK   length of WORK, WORK >= MAX(6,M+N)
+*  LWORK  (input) INTEGER 
+*         length of WORK, WORK >= MAX(6,M+N)
 *
 *  INFO    (output) INTEGER
 *          = 0 : successful exit.
@@ -249,6 +250,7 @@
 *          > 0 : SGESVJ did not converge in the maximal allowed number (30)
 *                of sweeps. The output may still be useful. See the
 *                description of WORK.
+*
 *  =====================================================================
 *
 *     .. Local Parameters ..
@@ -278,6 +280,7 @@
       INTRINSIC          ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT
 *     ..
 *     .. External Functions ..
+*     ..
 *     from BLAS
       REAL               SDOT, SNRM2
       EXTERNAL           SDOT, SNRM2
@@ -290,6 +293,7 @@
       EXTERNAL           LSAME
 *     ..
 *     .. External Subroutines ..
+*     ..
 *     from BLAS
       EXTERNAL           SAXPY, SCOPY, SROTM, SSCAL, SSWAP
 *     from LAPACK
@@ -370,6 +374,7 @@
       ROOTSFMIN = SQRT( SFMIN )
       SMALL = SFMIN / EPSLN
       BIG = SLAMCH( 'Overflow' )
+*     BIG         = ONE    / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
       LARGE = BIG / SQRT( FLOAT( M*N ) )
       BIGTHETA = ONE / ROOTEPS
@@ -378,7 +383,7 @@
       ROOTTOL = SQRT( TOL )
 *
       IF( FLOAT( M )*EPSLN.GE.ONE ) THEN
-         INFO = -5
+         INFO = -4
          CALL XERBLA( 'SGESVJ', -INFO )
          RETURN
       END IF
@@ -689,6 +694,7 @@
 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 *
       DO 1993 i = 1, NSWEEP
+*
 *     .. go go go ...
 *
          MXAAPQ = ZERO

diff --git a/SRC/sgsvj0.f b/SRC/sgsvj0.f
index 1279162..07c4d03 100644 (file)
--- a/SRC/sgsvj0.f
+++ b/SRC/sgsvj0.f
@@ -1,11 +1,11 @@
       SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
      +                   SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 *
-*  -- LAPACK routine (version 3.3.0)                                    --
+*  -- LAPACK routine (version 3.3.1)                                  --
 *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
-*     November 2010
+*     January 2011
 *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
@@ -222,7 +222,6 @@
       BIGTHETA = ONE / ROOTEPS
       ROOTTOL = SQRT( TOL )
 *
-*
 *     .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
 *
       EMPTSW = ( N*( N-1 ) ) / 2