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

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

version 1.3, 2016/08/27 15:34:48 version 1.4, 2017/06/17 10:54:13
Line 1 Line 1
 *> \brief \b ZGSVJ0 pre-processor for the routine zgesvj.  *> \brief <b> ZGSVJ0 pre-processor for the routine zgesvj. </b>
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download ZGSVJ0 + dependencies   *> Download ZGSVJ0 + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgsvj0.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgsvj0.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgsvj0.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgsvj0.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgsvj0.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgsvj0.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,  *       SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
 *                          SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )  *                          SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP  *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
 *       DOUBLE PRECISION   EPS, SFMIN, TOL  *       DOUBLE PRECISION   EPS, SFMIN, TOL
Line 30 Line 30
 *       COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )  *       COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
 *       DOUBLE PRECISION   SVA( N )  *       DOUBLE PRECISION   SVA( N )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 112 Line 112
 *>          the matrix A*diag(D).  *>          the matrix A*diag(D).
 *>          On exit, SVA contains the Euclidean norms of the columns of  *>          On exit, SVA contains the Euclidean norms of the columns of
 *>          the matrix A_onexit*diag(D_onexit).  *>          the matrix A_onexit*diag(D_onexit).
   *> \endverbatim
 *>  *>
 *> \param[in] MV  *> \param[in] MV
 *> \verbatim  *> \verbatim
Line 187 Line 188
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  *
 *> \date June 2016  *> \date June 2016
 *  *
Line 202 Line 203
 *> ZGSVJ0 is used just to enable ZGESVJ to call a simplified version of  *> ZGSVJ0 is used just to enable ZGESVJ to call a simplified version of
 *> itself to work on a submatrix of the original matrix.  *> itself to work on a submatrix of the original matrix.
 *>  *>
 *> Contributors:  *> Contributor:
 * =============  * =============
 *>  *>
 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)  *> Zlatko Drmac (Zagreb, Croatia)
 *>  *>
 *> Bugs, Examples and Comments:  *> \par Bugs, Examples and Comments:
 * ============================  * ============================
 *>  *>
 *> Please report all bugs and send interesting test examples and comments to  *> Please report all bugs and send interesting test examples and comments to
Line 217 Line 218
       SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,        SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
      $                   SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )       $                   SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.6.1) --  *  -- LAPACK computational routine (version 3.7.0) --
 *  -- 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..--
 *     June 2016  *     June 2016
Line 230 Line 231
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )        COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
       DOUBLE PRECISION   SVA( N )         DOUBLE PRECISION   SVA( N )
 *     ..  *     ..
 *  *
 *  =====================================================================  *  =====================================================================
Line 254 Line 255
 *     ..  *     ..
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC ABS, DMAX1, DCONJG, DBLE, MIN0, DSIGN, DSQRT        INTRINSIC ABS, MAX, CONJG, DBLE, MIN, SIGN, SQRT
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       DOUBLE PRECISION   DZNRM2        DOUBLE PRECISION   DZNRM2
Line 287 Line 288
          INFO = -5           INFO = -5
       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN        ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
          INFO = -8           INFO = -8
       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 = -10           INFO = -10
       ELSE IF( TOL.LE.EPS ) THEN        ELSE IF( TOL.LE.EPS ) THEN
Line 313 Line 314
       END IF        END IF
       RSVEC = RSVEC .OR. APPLV        RSVEC = RSVEC .OR. APPLV
   
       ROOTEPS = DSQRT( EPS )        ROOTEPS = SQRT( EPS )
       ROOTSFMIN = DSQRT( SFMIN )        ROOTSFMIN = SQRT( SFMIN )
       SMALL = SFMIN / EPS        SMALL = SFMIN / EPS
       BIG = ONE / SFMIN        BIG = ONE / SFMIN
       ROOTBIG = ONE / ROOTSFMIN        ROOTBIG = ONE / ROOTSFMIN
       BIGTHETA = ONE / ROOTEPS        BIGTHETA = ONE / ROOTEPS
       ROOTTOL = DSQRT( TOL )        ROOTTOL = SQRT( TOL )
 *  *
 *     .. Row-cyclic Jacobi SVD algorithm with column pivoting ..  *     .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
 *  *
Line 337 Line 338
 *     The boundaries are determined dynamically, based on the number of  *     The boundaries are determined dynamically, based on the number of
 *     pivots above a threshold.  *     pivots above a threshold.
 *  *
       KBL = MIN0( 8, N )        KBL = MIN( 8, N )
 *[TP] KBL is a tuning parameter that defines the tile size in the  *[TP] KBL is a tuning parameter that defines the tile size in the
 *     tiling of the p-q loops of pivot pairs. In general, an optimal  *     tiling of the p-q loops of pivot pairs. In general, an optimal
 *     value of KBL depends on the matrix dimensions and on the  *     value of KBL depends on the matrix dimensions and on the
Line 349 Line 350
       BLSKIP = KBL**2        BLSKIP = KBL**2
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.  *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 *  *
       ROWSKIP = MIN0( 5, KBL )        ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.  *[TP] ROWSKIP is a tuning parameter.
 *  *
       LKAHEAD = 1        LKAHEAD = 1
Line 383 Line 384
 *  *
             igl = ( ibr-1 )*KBL + 1              igl = ( ibr-1 )*KBL + 1
 *  *
             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )              DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
 *  *
                igl = igl + ir1*KBL                 igl = igl + ir1*KBL
 *  *
                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )                 DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
 *  *
 *     .. de Rijk's pivoting  *     .. de Rijk's pivoting
 *  *
                   q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1                    q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
                   IF( p.NE.q ) THEN                    IF( p.NE.q ) THEN
                      CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )                       CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
                      IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1,                         IF( RSVEC )CALL ZSWAP( 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 )
Line 418 Line 419
 *        If properly implemented DZNRM2 is available, the IF-THEN-ELSE-END IF  *        If properly implemented DZNRM2 is available, the IF-THEN-ELSE-END IF
 *        below should be replaced with "AAPP = DZNRM2( M, A(1,p), 1 )".  *        below should be replaced with "AAPP = DZNRM2( M, A(1,p), 1 )".
 *  *
                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.                            IF( ( SVA( p ).LT.ROOTBIG ) .AND.
      $                    ( SVA( p ).GT.ROOTSFMIN ) ) THEN       $                    ( SVA( p ).GT.ROOTSFMIN ) ) THEN
                         SVA( p ) = DZNRM2( M, A( 1, p ), 1 )                          SVA( p ) = DZNRM2( M, A( 1, p ), 1 )
                      ELSE                       ELSE
                         TEMP1 = ZERO                          TEMP1 = ZERO
                         AAPP = ONE                          AAPP = ONE
                         CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )                          CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*DSQRT( AAPP )                          SVA( p ) = TEMP1*SQRT( AAPP )
                      END IF                       END IF
                      AAPP = SVA( p )                       AAPP = SVA( p )
                   ELSE                    ELSE
Line 436 Line 437
 *  *
                      PSKIPPED = 0                       PSKIPPED = 0
 *  *
                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )                       DO 2002 q = p + 1, MIN( igl+KBL-1, N )
 *  *
                         AAQQ = SVA( q )                          AAQQ = SVA( q )
 *  *
Line 446 Line 447
                            IF( AAQQ.GE.ONE ) THEN                             IF( AAQQ.GE.ONE ) THEN
                               ROTOK = ( SMALL*AAPP ).LE.AAQQ                                ROTOK = ( SMALL*AAPP ).LE.AAQQ
                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN                                IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,                                    AAPQ = ( ZDOTC( M, A( 1, p ), 1,
      $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP       $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL ZCOPY( M, A( 1, p ), 1,                                      CALL ZCOPY( M, A( 1, p ), 1,
      $                                        WORK, 1 )       $                                        WORK, 1 )
                                  CALL ZLASCL( 'G', 0, 0, AAPP, ONE,                                    CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
      $                                M, 1, WORK, LDA, IERR )       $                                M, 1, WORK, LDA, IERR )
                                  AAPQ = ZDOTC( M, WORK, 1,                                   AAPQ = ZDOTC( M, WORK, 1,
      $                                   A( 1, q ), 1 ) / AAQQ       $                                   A( 1, q ), 1 ) / AAQQ
Line 459 Line 460
                            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 = ( ZDOTC( M, A( 1, p ), 1,                                    AAPQ = ( ZDOTC( M, A( 1, p ), 1,
      $                                    A( 1, q ), 1 ) / AAQQ ) / AAPP       $                                    A( 1, q ), 1 ) / AAPP ) / AAQQ
                               ELSE                                ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,                                      CALL ZCOPY( M, A( 1, q ), 1,
      $                                        WORK, 1 )       $                                        WORK, 1 )
                                  CALL ZLASCL( 'G', 0, 0, AAQQ,                                   CALL ZLASCL( 'G', 0, 0, AAQQ,
      $                                         ONE, M, 1,       $                                         ONE, M, 1,
      $                                         WORK, LDA, IERR )       $                                         WORK, LDA, IERR )
                                  AAPQ = ZDOTC( M, A( 1, p ), 1,                                      AAPQ = ZDOTC( M, A( 1, p ), 1,
      $                                   WORK, 1 ) / AAPP       $                                   WORK, 1 ) / AAPP
                               END IF                                END IF
                            END IF                             END IF
 *  *
                            OMPQ = AAPQ / ABS(AAPQ)   *                           AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
 *                           AAPQ = AAPQ * DCONJG( CWORK(p) ) * CWORK(q)                              AAPQ1  = -ABS(AAPQ)
                            AAPQ1  = -ABS(AAPQ)                              MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
                            MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )  
 *  *
 *        TO rotate or NOT to rotate, THAT is the question ...  *        TO rotate or NOT to rotate, THAT is the question ...
 *  *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN                             IF( ABS( AAPQ1 ).GT.TOL ) THEN
                                 OMPQ = AAPQ / ABS(AAPQ)
 *  *
 *           .. rotate  *           .. rotate
 *[RTD]      ROTATED = ROTATED + ONE  *[RTD]      ROTATED = ROTATED + ONE
Line 497 Line 498
                                  THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1                                   THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1
 *  *
                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN                                   IF( ABS( THETA ).GT.BIGTHETA ) THEN
 *   *
                                     T  = HALF / THETA                                      T  = HALF / THETA
                                     CS = ONE                                      CS = ONE
   
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,                                      CALL ZROT( M, A(1,p), 1, A(1,q), 1,
      $                                          CS, DCONJG(OMPQ)*T )       $                                          CS, CONJG(OMPQ)*T )
                                     IF ( RSVEC ) THEN                                      IF ( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,                                           CALL ZROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )       $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF                                      END IF
                                       
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                       SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )       $                                          ONE+T*APOAQ*AAPQ1 ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*SQRT( MAX( ZERO,
      $                                          ONE-T*AQOAP*AAPQ1 ) )       $                                          ONE-T*AQOAP*AAPQ1 ) )
                                     MXSINJ = DMAX1( MXSINJ, ABS( T ) )                                      MXSINJ = MAX( MXSINJ, ABS( T ) )
 *  *
                                  ELSE                                   ELSE
 *  *
 *                 .. choose correct signum for THETA and rotate  *                 .. choose correct signum for THETA and rotate
 *  *
                                     THSIGN = -DSIGN( ONE, AAPQ1 )                                      THSIGN = -SIGN( ONE, AAPQ1 )
                                     T = ONE / ( THETA+THSIGN*                                             T = ONE / ( THETA+THSIGN*
      $                                   DSQRT( ONE+THETA*THETA ) )       $                                   SQRT( ONE+THETA*THETA ) )
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
 *  *
                                     MXSINJ = DMAX1( MXSINJ, ABS( SN ) )                                      MXSINJ = MAX( MXSINJ, ABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )       $                                          ONE+T*APOAQ*AAPQ1 ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                        AAPP = AAPP*SQRT( MAX( ZERO,
      $                                      ONE-T*AQOAP*AAPQ1 ) )       $                                      ONE-T*AQOAP*AAPQ1 ) )
 *  *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,                                      CALL ZROT( M, A(1,p), 1, A(1,q), 1,
      $                                          CS, DCONJG(OMPQ)*SN )       $                                          CS, CONJG(OMPQ)*SN )
                                     IF ( RSVEC ) THEN                                      IF ( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,                                           CALL ZROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )       $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF                                         END IF
                                  END IF                                    END IF
                                  D(p) = -D(q) * OMPQ                                    D(p) = -D(q) * OMPQ
 *  *
                                  ELSE                                   ELSE
 *              .. have to use modified Gram-Schmidt like transformation  *              .. have to use modified Gram-Schmidt like transformation
Line 552 Line 553
      $                                       A( 1, q ), 1 )       $                                       A( 1, q ), 1 )
                                  CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,                                   CALL ZLASCL( '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*SQRT( MAX( ZERO,
      $                                      ONE-AAPQ1*AAPQ1 ) )       $                                      ONE-AAPQ1*AAPQ1 ) )
                                  MXSINJ = DMAX1( MXSINJ, SFMIN )                                   MXSINJ = MAX( MXSINJ, SFMIN )
                               END IF                                END IF
 *           END IF ROTOK THEN ... ELSE  *           END IF ROTOK THEN ... ELSE
 *  *
Line 571 Line 572
                                     AAQQ = ONE                                      AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,                                      CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )       $                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )                                      SVA( q ) = T*SQRT( AAQQ )
                                  END IF                                   END IF
                               END IF                                END IF
                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN                                IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
Line 583 Line 584
                                     AAPP = ONE                                      AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,                                      CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )       $                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )                                      AAPP = T*SQRT( AAPP )
                                  END IF                                   END IF
                                  SVA( p ) = AAPP                                   SVA( p ) = AAPP
                               END IF                                END IF
Line 618 Line 619
                   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 + MIN( igl+KBL-1, N ) - p
                   END IF                    END IF
 *  *
  2001          CONTINUE   2001          CONTINUE
Line 638 Line 639
 *        doing the block at ( ibr, jbc )  *        doing the block at ( ibr, jbc )
 *  *
                IJBLSK = 0                 IJBLSK = 0
                DO 2100 p = igl, MIN0( igl+KBL-1, N )                 DO 2100 p = igl, MIN( igl+KBL-1, N )
 *  *
                   AAPP = SVA( p )                    AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN                    IF( AAPP.GT.ZERO ) THEN
 *  *
                      PSKIPPED = 0                       PSKIPPED = 0
 *  *
                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )                       DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *  *
                         AAQQ = SVA( q )                          AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN                          IF( AAQQ.GT.ZERO ) THEN
Line 662 Line 663
                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP                                   ROTOK = ( SMALL*AAQQ ).LE.AAPP
                               END IF                                END IF
                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN                                IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,                                    AAPQ = ( ZDOTC( M, A( 1, p ), 1,
      $                                  A( 1, q ), 1 ) / AAQQ ) / AAPP       $                                  A( 1, q ), 1 ) / AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL ZCOPY( M, A( 1, p ), 1,                                   CALL ZCOPY( M, A( 1, p ), 1,
Line 680 Line 681
                                  ROTOK = AAQQ.LE.( AAPP / SMALL )                                   ROTOK = AAQQ.LE.( AAPP / SMALL )
                               END IF                                END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN                                IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,                                    AAPQ = ( ZDOTC( M, A( 1, p ), 1,
      $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP       $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
        $                                               / MIN(AAQQ,AAPP)
                               ELSE                                ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,                                   CALL ZCOPY( M, A( 1, q ), 1,
      $                                       WORK, 1 )       $                                       WORK, 1 )
Line 693 Line 695
                               END IF                                END IF
                            END IF                             END IF
 *  *
                            OMPQ = AAPQ / ABS(AAPQ)   *                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
 *                           AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)     
                            AAPQ1  = -ABS(AAPQ)                             AAPQ1  = -ABS(AAPQ)
                            MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )                             MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *  *
 *        TO rotate or NOT to rotate, THAT is the question ...  *        TO rotate or NOT to rotate, THAT is the question ...
 *  *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN                             IF( ABS( AAPQ1 ).GT.TOL ) THEN
                                 OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0                                NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1  *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0                                PSKIPPED = 0
Line 715 Line 717
 *  *
                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN                                   IF( ABS( THETA ).GT.BIGTHETA ) THEN
                                     T  = HALF / THETA                                      T  = HALF / THETA
                                     CS = ONE                                       CS = ONE
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,                                      CALL ZROT( M, A(1,p), 1, A(1,q), 1,
      $                                          CS, DCONJG(OMPQ)*T )       $                                          CS, CONJG(OMPQ)*T )
                                     IF( RSVEC ) THEN                                      IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,                                           CALL ZROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )       $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF                                      END IF
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )       $                                         ONE+T*APOAQ*AAPQ1 ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )       $                                     ONE-T*AQOAP*AAPQ1 ) )
                                     MXSINJ = DMAX1( MXSINJ, ABS( T ) )                                      MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE                                   ELSE
 *  *
 *                 .. choose correct signum for THETA and rotate  *                 .. choose correct signum for THETA and rotate
 *  *
                                     THSIGN = -DSIGN( ONE, AAPQ1 )                                      THSIGN = -SIGN( ONE, AAPQ1 )
                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN                                      IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                                     T = ONE / ( THETA+THSIGN*                                      T = ONE / ( THETA+THSIGN*
      $                                  DSQRT( ONE+THETA*THETA ) )       $                                  SQRT( ONE+THETA*THETA ) )
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
                                     MXSINJ = DMAX1( MXSINJ, ABS( SN ) )                                      MXSINJ = MAX( MXSINJ, ABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )       $                                         ONE+T*APOAQ*AAPQ1 ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                        AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )       $                                         ONE-T*AQOAP*AAPQ1 ) )
 *  *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,                                      CALL ZROT( M, A(1,p), 1, A(1,q), 1,
      $                                          CS, DCONJG(OMPQ)*SN )        $                                          CS, CONJG(OMPQ)*SN )
                                     IF( RSVEC ) THEN                                      IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,                                           CALL ZROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )       $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF                                      END IF
                                  END IF                                   END IF
                                  D(p) = -D(q) * OMPQ                                   D(p) = -D(q) * OMPQ
Line 768 Line 770
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,                                      CALL ZLASCL( '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*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )       $                                         ONE-AAPQ1*AAPQ1 ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE                                 ELSE
                                    CALL ZCOPY( M, A( 1, q ), 1,                                     CALL ZCOPY( M, A( 1, q ), 1,
      $                                          WORK, 1 )       $                                          WORK, 1 )
Line 780 Line 782
                                     CALL ZLASCL( 'G', 0, 0, AAPP, ONE,                                      CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
      $                                           M, 1, A( 1, p ), LDA,       $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )       $                                           IERR )
                                     CALL ZAXPY( M, -DCONJG(AAPQ),                                       CALL ZAXPY( M, -CONJG(AAPQ),
      $                                   WORK, 1, A( 1, p ), 1 )       $                                   WORK, 1, A( 1, p ), 1 )
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAPP,                                      CALL ZLASCL( '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*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )       $                                         ONE-AAPQ1*AAPQ1 ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF                                 END IF
                               END IF                                END IF
 *           END IF ROTOK THEN ... ELSE  *           END IF ROTOK THEN ... ELSE
Line 804 Line 806
                                     AAQQ = ONE                                      AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,                                      CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )       $                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )                                      SVA( q ) = T*SQRT( AAQQ )
                                  END IF                                   END IF
                               END IF                                END IF
                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN                                IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
Line 816 Line 818
                                     AAPP = ONE                                      AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,                                      CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )       $                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )                                      AAPP = T*SQRT( AAPP )
                                  END IF                                   END IF
                                  SVA( p ) = AAPP                                   SVA( p ) = AAPP
                               END IF                                END IF
Line 855 Line 857
                   ELSE                    ELSE
 *  *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +                       IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
      $                   MIN0( jgl+KBL-1, N ) - jgl + 1       $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0                       IF( AAPP.LT.ZERO )NOTROT = 0
 *  *
                   END IF                    END IF
Line 866 Line 868
 *     end of the jbc-loop  *     end of the jbc-loop
  2011       CONTINUE   2011       CONTINUE
 *2011 bailed out of the jbc-loop  *2011 bailed out of the jbc-loop
             DO 2012 p = igl, MIN0( igl+KBL-1, N )              DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )                 SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE   2012       CONTINUE
 ***  ***
Line 881 Line 883
             T = ZERO              T = ZERO
             AAPP = ONE              AAPP = ONE
             CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )              CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*DSQRT( AAPP )              SVA( N ) = T*SQRT( AAPP )
          END IF           END IF
 *  *
 *     Additional steering devices  *     Additional steering devices
Line 889 Line 891
          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.SQRT( 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 909 Line 911
 *  *
       INFO = 0        INFO = 0
 * #:) INFO = 0 confirms successful iterations.  * #:) INFO = 0 confirms successful iterations.
  1995  CONTINUE   1995 CONTINUE
 *  *
 *     Sort the vector SVA() of column norms.  *     Sort the vector SVA() of column norms.
       DO 5991 p = 1, N - 1        DO 5991 p = 1, N - 1
Removed from v.1.3  
changed lines
  Added in v.1.4

CVSweb interface <joel.bertrand@systella.fr>