[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.5 and 1.6

version 1.5, 2010/12/21 13:53:26 version 1.6, 2011/07/22 07:38: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.3.0)                                    --  *  -- LAPACK routine (version 3.3.1)                                  --
 *  *
 *  -- 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                  --
 *     November 2010  *  -- April 2011                                                      --
 *  *
 *  -- 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 23 Line 23
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),        DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
      +                   WORK( LWORK )       $                   WORK( LWORK )
 *     ..  *     ..
 *  *
 *  Purpose  *  Purpose
Line 133 Line 133
 *                  referenced  *                  referenced
 *  *
 *  M       (input) INTEGER  *  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  *  N       (input) INTEGER
 *          The number of columns of the input matrix A.  *          The number of columns of the input matrix A.
Line 256 Line 256
 *     .. Local Parameters ..  *     .. Local Parameters ..
       DOUBLE PRECISION   ZERO, HALF, ONE, TWO        DOUBLE PRECISION   ZERO, HALF, ONE, TWO
       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,        PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
      +                   TWO = 2.0D0 )       $                   TWO = 2.0D0 )
       INTEGER            NSWEEP        INTEGER            NSWEEP
       PARAMETER          ( NSWEEP = 30 )        PARAMETER          ( NSWEEP = 30 )
 *     ..  *     ..
 *     .. 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, EPSLN, LARGE, MXAAPQ,       $                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,       $                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
      +                   SKL, 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,
      +                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,       $                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
      +                   SWBAND       $                   SWBAND
       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,        LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
      +                   RSVEC, UCTOL, UPPER       $                   RSVEC, UCTOL, UPPER
 *     ..  *     ..
 *     .. Local Arrays ..  *     .. Local Arrays ..
       DOUBLE PRECISION   FASTR( 5 )        DOUBLE PRECISION   FASTR( 5 )
Line 327 Line 327
       ELSE IF( MV.LT.0 ) THEN        ELSE IF( MV.LT.0 ) THEN
          INFO = -9           INFO = -9
       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.        ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
      +         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN       $         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
          INFO = -11           INFO = -11
       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN        ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
          INFO = -12           INFO = -12
Line 383 Line 383
       ROOTTOL = DSQRT( TOL )        ROOTTOL = DSQRT( TOL )
 *  *
       IF( DBLE( M )*EPSLN.GE.ONE ) THEN        IF( DBLE( M )*EPSLN.GE.ONE ) THEN
          INFO = -5           INFO = -4
          CALL XERBLA( 'DGESVJ', -INFO )           CALL XERBLA( 'DGESVJ', -INFO )
          RETURN           RETURN
       END IF        END IF
Line 518 Line 518
 *  *
       IF( N.EQ.1 ) THEN        IF( N.EQ.1 ) THEN
          IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, 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 / SKL           WORK( 1 ) = ONE / SKL
          IF( SVA( 1 ).GE.SFMIN ) THEN           IF( SVA( 1 ).GE.SFMIN ) THEN
             WORK( 2 ) = ONE              WORK( 2 ) = ONE
Line 538 Line 538
       SN = DSQRT( SFMIN / EPSLN )        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
          TEMP1 = DMIN1( BIG, TEMP1 / AAPP )           TEMP1 = DMIN1( BIG, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1  *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1  *         AAPP  = AAPP*TEMP1
Line 638 Line 638
 *     [+ + x x]                    [x x].             [x x]  *     [+ + x x]                    [x x].             [x x]
 *  *
             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, EPSLN, 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, EPSLN, 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, EPSLN, 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, EPSLN, 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,
      +                   EPSLN, 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, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),       $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )       $                   LWORK-N, IERR )
 *  *
 *  *
          ELSE IF( UPPER ) THEN           ELSE IF( UPPER ) THEN
 *  *
 *  *
             CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,              CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
      +                   EPSLN, 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,
      +                   EPSLN, 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, EPSLN, 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, EPSLN, 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 725 Line 725
                   IF( p.NE.q ) THEN                    IF( p.NE.q ) THEN
                      CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )                       CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
                      IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,                       IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
      +                                      V( 1, q ), 1 )       $                                      V( 1, q ), 1 )
                      TEMP1 = SVA( p )                       TEMP1 = SVA( p )
                      SVA( p ) = SVA( q )                       SVA( p ) = SVA( q )
                      SVA( q ) = TEMP1                       SVA( q ) = TEMP1
Line 749 Line 749
 *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".  *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".
 *  *
                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.                       IF( ( SVA( p ).LT.ROOTBIG ) .AND.
      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN       $                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
                         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
Line 777 Line 777
                               ROTOK = ( SMALL*AAPP ).LE.AAQQ                                ROTOK = ( SMALL*AAPP ).LE.AAQQ
                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN                                IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, p ), 1,                                   CALL DCOPY( M, A( 1, p ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAPP,                                   CALL DLASCL( 'G', 0, 0, AAPP,
      +                                        WORK( p ), M, 1,       $                                        WORK( p ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ       $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
                               END IF                                END IF
                            ELSE                             ELSE
                               ROTOK = AAPP.LE.( AAQQ / SMALL )                                ROTOK = AAPP.LE.( AAQQ / SMALL )
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN                                IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, q ), 1,                                   CALL DCOPY( M, A( 1, q ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAQQ,                                   CALL DLASCL( 'G', 0, 0, AAQQ,
      +                                        WORK( q ), M, 1,       $                                        WORK( q ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP       $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
                               END IF                                END IF
                            END IF                             END IF
 *  *
Line 824 Line 824
 *  *
                                  AQOAP = AAQQ / AAPP                                   AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ                                   APOAQ = AAPP / AAQQ
                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /                                   THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
      +                                   AAPQ  
 *  *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN                                   IF( DABS( THETA ).GT.BIGTHETA ) THEN
 *  *
                                     T = HALF / THETA                                      T = HALF / THETA
                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )                                      FASTR( 3 ) = T*WORK( p ) / WORK( q )
                                     FASTR( 4 ) = -T*WORK( q ) /                                      FASTR( 4 ) = -T*WORK( q ) /
      +                                           WORK( p )       $                                           WORK( p )
                                     CALL DROTM( M, A( 1, p ), 1,                                      CALL DROTM( M, A( 1, p ), 1,
      +                                          A( 1, q ), 1, FASTR )       $                                          A( 1, q ), 1, FASTR )
                                     IF( RSVEC )CALL DROTM( MVL,                                      IF( RSVEC )CALL DROTM( MVL,
      +                                              V( 1, p ), 1,       $                                              V( 1, p ), 1,
      +                                              V( 1, q ), 1,       $                                              V( 1, q ), 1,
      +                                              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( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( T ) )
 *  *
                                  ELSE                                   ELSE
Line 851 Line 850
 *  *
                                     THSIGN = -DSIGN( ONE, AAPQ )                                      THSIGN = -DSIGN( ONE, AAPQ )
                                     T = ONE / ( THETA+THSIGN*                                      T = ONE / ( THETA+THSIGN*
      +                                  DSQRT( ONE+THETA*THETA ) )       $                                  DSQRT( ONE+THETA*THETA ) )
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = DSQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
 *  *
                                     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( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     ONE-T*AQOAP*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 870 Line 869
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                           CALL DROTM( M, A( 1, p ), 1,                                            CALL DROTM( M, A( 1, p ), 1,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                FASTR )       $                                                FASTR )
                                           IF( RSVEC )CALL DROTM( MVL,                                            IF( RSVEC )CALL DROTM( MVL,
      +                                        V( 1, p ), 1, V( 1, q ),       $                                        V( 1, p ), 1, V( 1, q ),
      +                                        1, FASTR )       $                                        1, FASTR )
                                        ELSE                                         ELSE
                                           CALL DAXPY( M, -T*AQOAP,                                            CALL DAXPY( M, -T*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           CALL DAXPY( M, CS*SN*APOAQ,                                            CALL DAXPY( M, CS*SN*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q ) / CS                                            WORK( q ) = WORK( q ) / CS
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, -T*AQOAP,                                               CALL DAXPY( MVL, -T*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   CS*SN*APOAQ,       $                                                   CS*SN*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                           END IF                                            END IF
                                        END IF                                         END IF
                                     ELSE                                      ELSE
                                        IF( WORK( q ).GE.ONE ) THEN                                         IF( WORK( q ).GE.ONE ) THEN
                                           CALL DAXPY( M, T*APOAQ,                                            CALL DAXPY( M, T*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           CALL DAXPY( M, -CS*SN*AQOAP,                                            CALL DAXPY( M, -CS*SN*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           WORK( p ) = WORK( p ) / CS                                            WORK( p ) = WORK( p ) / CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, T*APOAQ,                                               CALL DAXPY( MVL, T*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                           END IF                                            END IF
                                        ELSE                                         ELSE
                                           IF( WORK( p ).GE.WORK( q ) )                                            IF( WORK( p ).GE.WORK( q ) )
      +                                        THEN       $                                        THEN
                                              CALL DAXPY( M, -T*AQOAP,                                               CALL DAXPY( M, -T*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              CALL DAXPY( M, CS*SN*APOAQ,                                               CALL DAXPY( M, CS*SN*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              WORK( p ) = WORK( p )*CS                                               WORK( p ) = WORK( p )*CS
                                              WORK( q ) = WORK( q ) / CS                                               WORK( q ) = WORK( q ) / CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -T*AQOAP,       $                                               -T*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               CS*SN*APOAQ,       $                                               CS*SN*APOAQ,
      +                                               V( 1, p ), 1,       $                                               V( 1, p ), 1,
      +                                               V( 1, q ), 1 )       $                                               V( 1, q ), 1 )
                                              END IF                                               END IF
                                           ELSE                                            ELSE
                                              CALL DAXPY( M, T*APOAQ,                                               CALL DAXPY( M, T*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              CALL DAXPY( M,                                               CALL DAXPY( M,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              WORK( p ) = WORK( p ) / CS                                               WORK( p ) = WORK( p ) / CS
                                              WORK( q ) = WORK( q )*CS                                               WORK( q ) = WORK( q )*CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               T*APOAQ, V( 1, p ),       $                                               T*APOAQ, V( 1, p ),
      +                                               1, V( 1, q ), 1 )       $                                               1, V( 1, q ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -CS*SN*AQOAP,       $                                               -CS*SN*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                              END IF                                               END IF
                                           END IF                                            END IF
                                        END IF                                         END IF
Line 961 Line 960
                               ELSE                                ELSE
 *              .. have to use modified Gram-Schmidt like transformation  *              .. have to use modified Gram-Schmidt like transformation
                                  CALL DCOPY( M, A( 1, p ), 1,                                   CALL DCOPY( M, A( 1, p ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,                                   CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
      +                                        1, WORK( N+1 ), LDA,       $                                        1, WORK( N+1 ), LDA,
      +                                        IERR )       $                                        IERR )
                                  CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,                                   CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
      +                                        1, A( 1, q ), LDA, IERR )       $                                        1, A( 1, q ), LDA, IERR )
                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )                                   TEMP1 = -AAPQ*WORK( p ) / WORK( q )
                                  CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,                                   CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,
      +                                       A( 1, q ), 1 )       $                                       A( 1, q ), 1 )
                                  CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,                                   CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
      +                                        1, A( 1, q ), LDA, IERR )       $                                        1, A( 1, q ), LDA, IERR )
                                  SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                   SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                      ONE-AAPQ*AAPQ ) )       $                                      ONE-AAPQ*AAPQ ) )
                                  MXSINJ = DMAX1( MXSINJ, SFMIN )                                   MXSINJ = DMAX1( MXSINJ, SFMIN )
                               END IF                                END IF
 *           END IF ROTOK THEN ... ELSE  *           END IF ROTOK THEN ... ELSE
Line 982 Line 981
 *           recompute SVA(q), SVA(p).  *           recompute SVA(q), SVA(p).
 *  *
                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )                                IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
      +                            THEN       $                            THEN
                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.                                   IF( ( AAQQ.LT.ROOTBIG ) .AND.
      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN       $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*                                      SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
      +                                         WORK( q )       $                                         WORK( q )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAQQ = ONE                                      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 )
                                  END IF                                   END IF
                               END IF                                END IF
                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN                                IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
                                  IF( ( AAPP.LT.ROOTBIG ) .AND.                                   IF( ( AAPP.LT.ROOTBIG ) .AND.
      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN       $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
                                     AAPP = DNRM2( M, A( 1, p ), 1 )*                                      AAPP = DNRM2( M, A( 1, p ), 1 )*
      +                                     WORK( p )       $                                     WORK( p )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAPP = ONE                                      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 )
                                  END IF                                   END IF
                                  SVA( p ) = AAPP                                   SVA( p ) = AAPP
Line 1023 Line 1022
                         END IF                          END IF
 *  *
                         IF( ( i.LE.SWBAND ) .AND.                          IF( ( i.LE.SWBAND ) .AND.
      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN       $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
                            IF( ir1.EQ.0 )AAPP = -AAPP                             IF( ir1.EQ.0 )AAPP = -AAPP
                            NOTROT = 0                             NOTROT = 0
                            GO TO 2103                             GO TO 2103
Line 1040 Line 1039
                   ELSE                    ELSE
                      SVA( p ) = AAPP                       SVA( p ) = AAPP
                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )                       IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p       $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
                   END IF                    END IF
 *  *
  2001          CONTINUE   2001          CONTINUE
Line 1085 Line 1084
                               END IF                                END IF
                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN                                IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, p ), 1,                                   CALL DCOPY( M, A( 1, p ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAPP,                                   CALL DLASCL( 'G', 0, 0, AAPP,
      +                                        WORK( p ), M, 1,       $                                        WORK( p ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ       $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
                               END IF                                END IF
                            ELSE                             ELSE
                               IF( AAPP.GE.AAQQ ) THEN                                IF( AAPP.GE.AAQQ ) THEN
Line 1104 Line 1103
                               END IF                                END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN                                IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, q ), 1,                                   CALL DCOPY( M, A( 1, q ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAQQ,                                   CALL DLASCL( 'G', 0, 0, AAQQ,
      +                                        WORK( q ), M, 1,       $                                        WORK( q ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP       $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
                               END IF                                END IF
                            END IF                             END IF
 *  *
Line 1131 Line 1130
 *  *
                                  AQOAP = AAQQ / AAPP                                   AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ                                   APOAQ = AAPP / AAQQ
                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /                                   THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
      +                                   AAPQ  
                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA                                   IF( AAQQ.GT.AAPP0 )THETA = -THETA
 *  *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN                                   IF( DABS( THETA ).GT.BIGTHETA ) THEN
                                     T = HALF / THETA                                      T = HALF / THETA
                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )                                      FASTR( 3 ) = T*WORK( p ) / WORK( q )
                                     FASTR( 4 ) = -T*WORK( q ) /                                      FASTR( 4 ) = -T*WORK( q ) /
      +                                           WORK( p )       $                                           WORK( p )
                                     CALL DROTM( M, A( 1, p ), 1,                                      CALL DROTM( M, A( 1, p ), 1,
      +                                          A( 1, q ), 1, FASTR )       $                                          A( 1, q ), 1, FASTR )
                                     IF( RSVEC )CALL DROTM( MVL,                                      IF( RSVEC )CALL DROTM( MVL,
      +                                              V( 1, p ), 1,       $                                              V( 1, p ), 1,
      +                                              V( 1, q ), 1,       $                                              V( 1, q ), 1,
      +                                              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( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( T ) )
                                  ELSE                                   ELSE
 *  *
Line 1158 Line 1156
                                     THSIGN = -DSIGN( ONE, AAPQ )                                      THSIGN = -DSIGN( ONE, AAPQ )
                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN                                      IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                                     T = ONE / ( THETA+THSIGN*                                      T = ONE / ( THETA+THSIGN*
      +                                  DSQRT( ONE+THETA*THETA ) )       $                                  DSQRT( ONE+THETA*THETA ) )
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = DSQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
                                     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( DMAX1( ZERO,                                       AAPP = AAPP*DSQRT( DMAX1( ZERO, 
      +                                     ONE-T*AQOAP*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 1177 Line 1175
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                           CALL DROTM( M, A( 1, p ), 1,                                            CALL DROTM( M, A( 1, p ), 1,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                FASTR )       $                                                FASTR )
                                           IF( RSVEC )CALL DROTM( MVL,                                            IF( RSVEC )CALL DROTM( MVL,
      +                                        V( 1, p ), 1, V( 1, q ),       $                                        V( 1, p ), 1, V( 1, q ),
      +                                        1, FASTR )       $                                        1, FASTR )
                                        ELSE                                         ELSE
                                           CALL DAXPY( M, -T*AQOAP,                                            CALL DAXPY( M, -T*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           CALL DAXPY( M, CS*SN*APOAQ,                                            CALL DAXPY( M, CS*SN*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, -T*AQOAP,                                               CALL DAXPY( MVL, -T*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   CS*SN*APOAQ,       $                                                   CS*SN*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                           END IF                                            END IF
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q ) / CS                                            WORK( q ) = WORK( q ) / CS
Line 1204 Line 1202
                                     ELSE                                      ELSE
                                        IF( WORK( q ).GE.ONE ) THEN                                         IF( WORK( q ).GE.ONE ) THEN
                                           CALL DAXPY( M, T*APOAQ,                                            CALL DAXPY( M, T*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           CALL DAXPY( M, -CS*SN*AQOAP,                                            CALL DAXPY( M, -CS*SN*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, T*APOAQ,                                               CALL DAXPY( MVL, T*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                           END IF                                            END IF
                                           WORK( p ) = WORK( p ) / CS                                            WORK( p ) = WORK( p ) / CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                        ELSE                                         ELSE
                                           IF( WORK( p ).GE.WORK( q ) )                                            IF( WORK( p ).GE.WORK( q ) )
      +                                        THEN       $                                        THEN
                                              CALL DAXPY( M, -T*AQOAP,                                               CALL DAXPY( M, -T*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              CALL DAXPY( M, CS*SN*APOAQ,                                               CALL DAXPY( M, CS*SN*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              WORK( p ) = WORK( p )*CS                                               WORK( p ) = WORK( p )*CS
                                              WORK( q ) = WORK( q ) / CS                                               WORK( q ) = WORK( q ) / CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -T*AQOAP,       $                                               -T*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               CS*SN*APOAQ,       $                                               CS*SN*APOAQ,
      +                                               V( 1, p ), 1,       $                                               V( 1, p ), 1,
      +                                               V( 1, q ), 1 )       $                                               V( 1, q ), 1 )
                                              END IF                                               END IF
                                           ELSE                                            ELSE
                                              CALL DAXPY( M, T*APOAQ,                                               CALL DAXPY( M, T*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              CALL DAXPY( M,                                               CALL DAXPY( M,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              WORK( p ) = WORK( p ) / CS                                               WORK( p ) = WORK( p ) / CS
                                              WORK( q ) = WORK( q )*CS                                               WORK( q ) = WORK( q )*CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               T*APOAQ, V( 1, p ),       $                                               T*APOAQ, V( 1, p ),
      +                                               1, V( 1, q ), 1 )       $                                               1, V( 1, q ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -CS*SN*AQOAP,       $                                               -CS*SN*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                              END IF                                               END IF
                                           END IF                                            END IF
                                        END IF                                         END IF
Line 1268 Line 1266
                               ELSE                                ELSE
                                  IF( AAPP.GT.AAQQ ) THEN                                   IF( AAPP.GT.AAQQ ) THEN
                                     CALL DCOPY( M, A( 1, p ), 1,                                      CALL DCOPY( M, A( 1, p ), 1,
      +                                          WORK( N+1 ), 1 )       $                                          WORK( N+1 ), 1 )
                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,                                      CALL DLASCL( 'G', 0, 0, AAPP, ONE,
      +                                           M, 1, WORK( N+1 ), LDA,       $                                           M, 1, WORK( N+1 ), LDA,
      +                                           IERR )       $                                           IERR )
                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,                                      CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
      +                                           M, 1, A( 1, q ), LDA,       $                                           M, 1, A( 1, q ), LDA,
      +                                           IERR )       $                                           IERR )
                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )                                      TEMP1 = -AAPQ*WORK( p ) / WORK( q )
                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),                                      CALL DAXPY( M, TEMP1, WORK( N+1 ),
      +                                          1, A( 1, q ), 1 )       $                                          1, A( 1, q ), 1 )
                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,                                      CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
      +                                           M, 1, A( 1, q ), LDA,       $                                           M, 1, A( 1, q ), LDA,
      +                                           IERR )       $                                           IERR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE-AAPQ*AAPQ ) )       $                                         ONE-AAPQ*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = DMAX1( MXSINJ, SFMIN )
                                  ELSE                                   ELSE
                                     CALL DCOPY( M, A( 1, q ), 1,                                      CALL DCOPY( M, A( 1, q ), 1,
      +                                          WORK( N+1 ), 1 )       $                                          WORK( N+1 ), 1 )
                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,                                      CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
      +                                           M, 1, WORK( N+1 ), LDA,       $                                           M, 1, WORK( N+1 ), LDA,
      +                                           IERR )       $                                           IERR )
                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,                                      CALL DLASCL( 'G', 0, 0, AAPP, ONE,
      +                                           M, 1, A( 1, p ), LDA,       $                                           M, 1, A( 1, p ), LDA,
      +                                           IERR )       $                                           IERR )
                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )                                      TEMP1 = -AAPQ*WORK( q ) / WORK( p )
                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),                                      CALL DAXPY( M, TEMP1, WORK( N+1 ),
      +                                          1, A( 1, p ), 1 )       $                                          1, A( 1, p ), 1 )
                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,                                      CALL DLASCL( 'G', 0, 0, ONE, AAPP,
      +                                           M, 1, A( 1, p ), LDA,       $                                           M, 1, A( 1, p ), LDA,
      +                                           IERR )       $                                           IERR )
                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,                                      SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
      +                                         ONE-AAPQ*AAPQ ) )       $                                         ONE-AAPQ*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = DMAX1( MXSINJ, SFMIN )
                                  END IF                                   END IF
                               END IF                                END IF
Line 1309 Line 1307
 *           In the case of cancellation in updating SVA(q)  *           In the case of cancellation in updating SVA(q)
 *           .. recompute SVA(q)  *           .. recompute SVA(q)
                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )                                IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
      +                            THEN       $                            THEN
                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.                                   IF( ( AAQQ.LT.ROOTBIG ) .AND.
      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN       $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*                                      SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
      +                                         WORK( q )       $                                         WORK( q )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAQQ = ONE                                      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 )
                                  END IF                                   END IF
                               END IF                                END IF
                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN                                IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
                                  IF( ( AAPP.LT.ROOTBIG ) .AND.                                   IF( ( AAPP.LT.ROOTBIG ) .AND.
      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN       $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
                                     AAPP = DNRM2( M, A( 1, p ), 1 )*                                      AAPP = DNRM2( M, A( 1, p ), 1 )*
      +                                     WORK( p )       $                                     WORK( p )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAPP = ONE                                      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 )
                                  END IF                                   END IF
                                  SVA( p ) = AAPP                                   SVA( p ) = AAPP
Line 1350 Line 1348
                         END IF                          END IF
 *  *
                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )                          IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
      +                      THEN       $                      THEN
                            SVA( p ) = AAPP                             SVA( p ) = AAPP
                            NOTROT = 0                             NOTROT = 0
                            GO TO 2011                             GO TO 2011
                         END IF                          END IF
                         IF( ( i.LE.SWBAND ) .AND.                          IF( ( i.LE.SWBAND ) .AND.
      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN       $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
                            AAPP = -AAPP                             AAPP = -AAPP
                            NOTROT = 0                             NOTROT = 0
                            GO TO 2203                             GO TO 2203
Line 1371 Line 1369
                   ELSE                    ELSE
 *  *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +                       IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
      +                   MIN0( jgl+KBL-1, N ) - jgl + 1       $                   MIN0( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0                       IF( AAPP.LT.ZERO )NOTROT = 0
 *  *
                   END IF                    END IF
Line 1391 Line 1389
 *  *
 *     .. update SVA(N)  *     .. update SVA(N)
          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )           IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
      +       THEN       $       THEN
             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
Line 1403 Line 1401
 *     Additional steering devices  *     Additional steering devices
 *  *
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.           IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      +       ( ISWROT.LE.N ) ) )SWBAND = i       $       ( ISWROT.LE.N ) ) )SWBAND = i
 *  *
          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*           IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
      +       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN       $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994              GO TO 1994
          END IF           END IF
 *  *
Line 1479 Line 1477
 *  *
 *     Undo scaling, if necessary (and possible).  *     Undo scaling, if necessary (and possible).
       IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /        IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
      +    SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.       $    SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
      +    ( SFMIN / SKL) ) ) ) THEN       $    ( SFMIN / SKL) ) ) ) THEN
          DO 2400 p = 1, N           DO 2400 p = 1, N
             SVA( p ) = SKL*SVA( p )              SVA( p ) = SKL*SVA( p )
  2400    CONTINUE   2400    CONTINUE
Removed from v.1.5  
changed lines
  Added in v.1.6

CVSweb interface <joel.bertrand@systella.fr>