[BACK]Return to dgesvj.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dgesvj.f between versions 1.3 and 1.4

version 1.3, 2010/08/13 21:03:45 version 1.4, 2010/12/21 13:48:05
Line 1 Line 1
       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,        SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
      +                   LDV, WORK, LWORK, INFO )       +                   LDV, WORK, LWORK, INFO )
 *  *
 *  -- LAPACK routine (version 3.2.2)                                    --  *  -- LAPACK routine (version 3.3.0)                                    --
 *  *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --  *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --  *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
 *  -- June 2010                                                       --  *     November 2010
 *  *
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Line 222 Line 222
 *                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').  *                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
 *                    It is required that CTOL >= ONE, i.e. it is not  *                    It is required that CTOL >= ONE, i.e. it is not
 *                    allowed to force the routine to obtain orthogonality  *                    allowed to force the routine to obtain orthogonality
 *                    below EPSILON.  *                    below EPS.
 *          On exit :  *          On exit :
 *          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)  *          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
 *                    are the computed singular values of A.  *                    are the computed singular values of A.
Line 262 Line 262
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,        DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
      +                   BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,       +                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,       +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
      +                   SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,       +                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
      +                   THSIGN, TOL       +                   THSIGN, TOL
       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,        INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,       +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
Line 368 Line 368
 *     ... and the machine dependent parameters are  *     ... and the machine dependent parameters are
 *[!]  (Make sure that DLAMCH() works properly on the target machine.)  *[!]  (Make sure that DLAMCH() works properly on the target machine.)
 *  *
       EPSILON = DLAMCH( 'Epsilon' )        EPSLN = DLAMCH( 'Epsilon' )
       ROOTEPS = DSQRT( EPSILON )        ROOTEPS = DSQRT( EPSLN )
       SFMIN = DLAMCH( 'SafeMinimum' )        SFMIN = DLAMCH( 'SafeMinimum' )
       ROOTSFMIN = DSQRT( SFMIN )        ROOTSFMIN = DSQRT( SFMIN )
       SMALL = SFMIN / EPSILON        SMALL = SFMIN / EPSLN
       BIG = DLAMCH( 'Overflow' )        BIG = DLAMCH( 'Overflow' )
 *     BIG         = ONE    / SFMIN  *     BIG         = ONE    / SFMIN
       ROOTBIG = ONE / ROOTSFMIN        ROOTBIG = ONE / ROOTSFMIN
       LARGE = BIG / DSQRT( DBLE( M*N ) )        LARGE = BIG / DSQRT( DBLE( M*N ) )
       BIGTHETA = ONE / ROOTEPS        BIGTHETA = ONE / ROOTEPS
 *  *
       TOL = CTOL*EPSILON        TOL = CTOL*EPSLN
       ROOTTOL = DSQRT( TOL )        ROOTTOL = DSQRT( TOL )
 *  *
       IF( DBLE( M )*EPSILON.GE.ONE ) THEN        IF( DBLE( M )*EPSLN.GE.ONE ) THEN
          INFO = -5           INFO = -5
          CALL XERBLA( 'DGESVJ', -INFO )           CALL XERBLA( 'DGESVJ', -INFO )
          RETURN           RETURN
Line 407 Line 407
 *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries  *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 *     in A are detected, the procedure returns with INFO=-6.  *     in A are detected, the procedure returns with INFO=-6.
 *  *
       SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) )        SKL= ONE / DSQRT( DBLE( M )*DBLE( N ) )
       NOSCALE = .TRUE.        NOSCALE = .TRUE.
       GOSCALE = .TRUE.        GOSCALE = .TRUE.
 *  *
Line 415 Line 415
 *        the input matrix is M-by-N lower triangular (trapezoidal)  *        the input matrix is M-by-N lower triangular (trapezoidal)
          DO 1874 p = 1, N           DO 1874 p = 1, N
             AAPP = ZERO              AAPP = ZERO
             AAQQ = ZERO              AAQQ = ONE
             CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )              CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN              IF( AAPP.GT.BIG ) THEN
                INFO = -6                 INFO = -6
Line 427 Line 427
                SVA( p ) = AAPP*AAQQ                 SVA( p ) = AAPP*AAQQ
             ELSE              ELSE
                NOSCALE = .FALSE.                 NOSCALE = .FALSE.
                SVA( p ) = AAPP*( AAQQ*SCALE )                 SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN                 IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.                    GOSCALE = .FALSE.
                   DO 1873 q = 1, p - 1                    DO 1873 q = 1, p - 1
                      SVA( q ) = SVA( q )*SCALE                       SVA( q ) = SVA( q )*SKL
  1873             CONTINUE   1873             CONTINUE
                END IF                 END IF
             END IF              END IF
Line 440 Line 440
 *        the input matrix is M-by-N upper triangular (trapezoidal)  *        the input matrix is M-by-N upper triangular (trapezoidal)
          DO 2874 p = 1, N           DO 2874 p = 1, N
             AAPP = ZERO              AAPP = ZERO
             AAQQ = ZERO              AAQQ = ONE
             CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )              CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN              IF( AAPP.GT.BIG ) THEN
                INFO = -6                 INFO = -6
Line 452 Line 452
                SVA( p ) = AAPP*AAQQ                 SVA( p ) = AAPP*AAQQ
             ELSE              ELSE
                NOSCALE = .FALSE.                 NOSCALE = .FALSE.
                SVA( p ) = AAPP*( AAQQ*SCALE )                 SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN                 IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.                    GOSCALE = .FALSE.
                   DO 2873 q = 1, p - 1                    DO 2873 q = 1, p - 1
                      SVA( q ) = SVA( q )*SCALE                       SVA( q ) = SVA( q )*SKL
  2873             CONTINUE   2873             CONTINUE
                END IF                 END IF
             END IF              END IF
Line 465 Line 465
 *        the input matrix is M-by-N general dense  *        the input matrix is M-by-N general dense
          DO 3874 p = 1, N           DO 3874 p = 1, N
             AAPP = ZERO              AAPP = ZERO
             AAQQ = ZERO              AAQQ = ONE
             CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )              CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN              IF( AAPP.GT.BIG ) THEN
                INFO = -6                 INFO = -6
Line 477 Line 477
                SVA( p ) = AAPP*AAQQ                 SVA( p ) = AAPP*AAQQ
             ELSE              ELSE
                NOSCALE = .FALSE.                 NOSCALE = .FALSE.
                SVA( p ) = AAPP*( AAQQ*SCALE )                 SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN                 IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.                    GOSCALE = .FALSE.
                   DO 3873 q = 1, p - 1                    DO 3873 q = 1, p - 1
                      SVA( q ) = SVA( q )*SCALE                       SVA( q ) = SVA( q )*SKL
  3873             CONTINUE   3873             CONTINUE
                END IF                 END IF
             END IF              END IF
  3874    CONTINUE   3874    CONTINUE
       END IF        END IF
 *  *
       IF( NOSCALE )SCALE = ONE        IF( NOSCALE )SKL= ONE
 *  *
 *     Move the smaller part of the spectrum from the underflow threshold  *     Move the smaller part of the spectrum from the underflow threshold
 *(!)  Start by determining the position of the nonzero entries of the  *(!)  Start by determining the position of the nonzero entries of the
Line 517 Line 517
 * #:) Quick return for one-column matrix  * #:) Quick return for one-column matrix
 *  *
       IF( N.EQ.1 ) THEN        IF( N.EQ.1 ) THEN
          IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,           IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
      +                           A( 1, 1 ), LDA, IERR )       +                           A( 1, 1 ), LDA, IERR )
          WORK( 1 ) = ONE / SCALE           WORK( 1 ) = ONE / SKL
          IF( SVA( 1 ).GE.SFMIN ) THEN           IF( SVA( 1 ).GE.SFMIN ) THEN
             WORK( 2 ) = ONE              WORK( 2 ) = ONE
          ELSE           ELSE
Line 535 Line 535
 *     Protect small singular values from underflow, and try to  *     Protect small singular values from underflow, and try to
 *     avoid underflows/overflows in computing Jacobi rotations.  *     avoid underflows/overflows in computing Jacobi rotations.
 *  *
       SN = DSQRT( SFMIN / EPSILON )        SN = DSQRT( SFMIN / EPSLN )
       TEMP1 = DSQRT( BIG / DBLE( N ) )        TEMP1 = DSQRT( BIG / DBLE( N ) )
       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.        IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN       +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
Line 563 Line 563
       IF( TEMP1.NE.ONE ) THEN        IF( TEMP1.NE.ONE ) THEN
          CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )           CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
       END IF        END IF
       SCALE = TEMP1*SCALE        SKL= TEMP1*SKL
       IF( SCALE.NE.ONE ) THEN        IF( SKL.NE.ONE ) THEN
          CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )           CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
          SCALE = ONE / SCALE           SKL= ONE / SKL
       END IF        END IF
 *  *
 *     Row-cyclic Jacobi SVD algorithm with column pivoting  *     Row-cyclic Jacobi SVD algorithm with column pivoting
Line 639 Line 639
 *  *
             CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,              CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,       +                   WORK( N34+1 ), SVA( N34+1 ), MVL,
      +                   V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,       +                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
      +                   2, WORK( N+1 ), LWORK-N, IERR )       +                   2, WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,              CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,       +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,       +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
      +                   WORK( N+1 ), LWORK-N, IERR )       +                   WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,              CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,       +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,       +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )       +                   WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,              CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,       +                   WORK( N4+1 ), SVA( N4+1 ), MVL,
      +                   V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,       +                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )       +                   WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,              CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,       +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )       +                   IERR )
 *  *
             CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,              CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),       +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )       +                   LWORK-N, IERR )
 *  *
 *  *
Line 670 Line 670
 *  *
 *  *
             CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,              CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
      +                   EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,       +                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
      +                   IERR )       +                   IERR )
 *  *
             CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),              CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,       +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,       +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )       +                   IERR )
 *  *
             CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,              CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),       +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )       +                   LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,              CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,       +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,       +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )       +                   WORK( N+1 ), LWORK-N, IERR )
   
          END IF           END IF
Line 753 Line 753
                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )                          SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
                      ELSE                       ELSE
                         TEMP1 = ZERO                          TEMP1 = ZERO
                         AAPP = ZERO                          AAPP = ONE
                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )                          CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )                          SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
                      END IF                       END IF
Line 841 Line 841
      +                                              FASTR )       +                                              FASTR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )       +                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     AAPQ )       +                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( T ) )
 *  *
                                  ELSE                                   ELSE
Line 989 Line 989
      +                                         WORK( q )       +                                         WORK( q )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAQQ = ZERO                                      AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,                                      CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )       +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )                                      SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
Line 1002 Line 1002
      +                                     WORK( p )       +                                     WORK( p )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAPP = ZERO                                      AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,                                      CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )       +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*WORK( p )                                      AAPP = T*DSQRT( AAPP )*WORK( p )
Line 1164 Line 1164
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )       +                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*                                      AAPP = AAPP*DSQRT( DMAX1( ZERO, 
      +                                     AAPQ )       +                                     ONE-T*AQOAP*AAPQ ) )
 *  *
                                     APOAQ = WORK( p ) / WORK( q )                                      APOAQ = WORK( p ) / WORK( q )
                                     AQOAP = WORK( q ) / WORK( p )                                      AQOAP = WORK( q ) / WORK( p )
Line 1316 Line 1316
      +                                         WORK( q )       +                                         WORK( q )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAQQ = ZERO                                      AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,                                      CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )       +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )                                      SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
Line 1329 Line 1329
      +                                     WORK( p )       +                                     WORK( p )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAPP = ZERO                                      AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,                                      CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )       +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*WORK( p )                                      AAPP = T*DSQRT( AAPP )*WORK( p )
Line 1395 Line 1395
             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )              SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
          ELSE           ELSE
             T = ZERO              T = ZERO
             AAPP = ZERO              AAPP = ONE
             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )              CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*DSQRT( AAPP )*WORK( N )              SVA( N ) = T*DSQRT( AAPP )*WORK( N )
          END IF           END IF
Line 1446 Line 1446
          END IF           END IF
          IF( SVA( p ).NE.ZERO ) THEN           IF( SVA( p ).NE.ZERO ) THEN
             N4 = N4 + 1              N4 = N4 + 1
             IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1              IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
          END IF           END IF
  5991 CONTINUE   5991 CONTINUE
       IF( SVA( N ).NE.ZERO ) THEN        IF( SVA( N ).NE.ZERO ) THEN
          N4 = N4 + 1           N4 = N4 + 1
          IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1           IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
       END IF        END IF
 *  *
 *     Normalize the left singular vectors.  *     Normalize the left singular vectors.
Line 1478 Line 1478
       END IF        END IF
 *  *
 *     Undo scaling, if necessary (and possible).  *     Undo scaling, if necessary (and possible).
       IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /        IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
      +    SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.       +    SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
      +    ( SFMIN / SCALE ) ) ) ) THEN       +    ( SFMIN / SKL) ) ) ) THEN
          DO 2400 p = 1, N           DO 2400 p = 1, N
             SVA( p ) = SCALE*SVA( p )              SVA( p ) = SKL*SVA( p )
  2400    CONTINUE   2400    CONTINUE
          SCALE = ONE           SKL= ONE
       END IF        END IF
 *  *
       WORK( 1 ) = SCALE        WORK( 1 ) = SKL
 *     The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE  *     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
 *     then some of the singular values may overflow or underflow and  *     then some of the singular values may overflow or underflow and
 *     the spectrum is given in this factored representation.  *     the spectrum is given in this factored representation.
 *  *
Removed from v.1.3  
changed lines
  Added in v.1.4

CVSweb interface <joel.bertrand@systella.fr>