[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.2 and 1.9

version 1.2, 2016/08/27 15:27:13 version 1.9, 2023/08/07 08:39:21
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
 *>          MV is INTEGER  *>          MV is INTEGER
 *>          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a  *>          If JOBV = 'A', then MV rows of V are post-multipled by a
 *>                           sequence of Jacobi rotations.  *>                           sequence of Jacobi rotations.
 *>          If JOBV = 'N',   then MV is not referenced.  *>          If JOBV = 'N',   then MV is not referenced.
 *> \endverbatim  *> \endverbatim
Line 124 Line 125
 *> \param[in,out] V  *> \param[in,out] V
 *> \verbatim  *> \verbatim
 *>          V is COMPLEX*16 array, dimension (LDV,N)  *>          V is COMPLEX*16 array, dimension (LDV,N)
 *>          If JOBV .EQ. 'V' then N rows of V are post-multipled by a  *>          If JOBV = 'V' then N rows of V are post-multipled by a
 *>                           sequence of Jacobi rotations.  *>                           sequence of Jacobi rotations.
 *>          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a  *>          If JOBV = 'A' then MV rows of V are post-multipled by a
 *>                           sequence of Jacobi rotations.  *>                           sequence of Jacobi rotations.
 *>          If JOBV = 'N',   then V is not referenced.  *>          If JOBV = 'N',   then V is not referenced.
 *> \endverbatim  *> \endverbatim
Line 135 Line 136
 *> \verbatim  *> \verbatim
 *>          LDV is INTEGER  *>          LDV is INTEGER
 *>          The leading dimension of the array V,  LDV >= 1.  *>          The leading dimension of the array V,  LDV >= 1.
 *>          If JOBV = 'V', LDV .GE. N.  *>          If JOBV = 'V', LDV >= N.
 *>          If JOBV = 'A', LDV .GE. MV.  *>          If JOBV = 'A', LDV >= MV.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] EPS  *> \param[in] EPS
Line 156 Line 157
 *>          TOL is DOUBLE PRECISION  *>          TOL is DOUBLE PRECISION
 *>          TOL is the threshold for Jacobi rotations. For a pair  *>          TOL is the threshold for Jacobi rotations. For a pair
 *>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is  *>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
 *>          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.  *>          applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] NSWEEP  *> \param[in] NSWEEP
Line 168 Line 169
 *>  *>
 *> \param[out] WORK  *> \param[out] WORK
 *> \verbatim  *> \verbatim
 *>          WORK is COMPLEX*16 array, dimension LWORK.  *>          WORK is COMPLEX*16 array, dimension (LWORK)
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LWORK  *> \param[in] LWORK
 *> \verbatim  *> \verbatim
 *>          LWORK is INTEGER  *>          LWORK is INTEGER
 *>          LWORK is the dimension of WORK. LWORK .GE. M.  *>          LWORK is the dimension of WORK. LWORK >= M.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] INFO  *> \param[out] INFO
 *> \verbatim  *> \verbatim
 *>          INFO is INTEGER  *>          INFO is INTEGER
 *>          = 0 : successful exit.  *>          = 0:  successful exit.
 *>          < 0 : if INFO = -i, then the i-th argument had an illegal value  *>          < 0:  if INFO = -i, then the i-th argument had an illegal value
 *> \endverbatim  *> \endverbatim
 *  *
 *  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  
 *  *
 *> \ingroup complex16OTHERcomputational  *> \ingroup complex16OTHERcomputational
 *>  *>
Line 202 Line 201
 *> 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 216
       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 --
 *  -- 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  
 *  *
       IMPLICIT NONE        IMPLICIT NONE
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
Line 230 Line 228
 *     ..  *     ..
 *     .. 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 252
 *     ..  *     ..
 *     ..  *     ..
 *     .. 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 267 Line 265
 *     .. External Subroutines ..  *     .. External Subroutines ..
 *     ..  *     ..
 *     from BLAS  *     from BLAS
       EXTERNAL           ZCOPY, ZROT, ZSWAP        EXTERNAL           ZCOPY, ZROT, ZSWAP, ZAXPY
 *     from LAPACK  *     from LAPACK
       EXTERNAL           ZLASCL, ZLASSQ, XERBLA        EXTERNAL           ZLASCL, ZLASSQ, XERBLA
 *     ..  *     ..
Line 287 Line 285
          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 311
       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 335
 *     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 347
       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 381
 *  *
             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 416
 *        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 434
 *  *
                      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 444
                            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 457
                            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 495
                                  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 550
      $                                       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 569
                                     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 581
                                     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 616
                   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 636
 *        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 660
                                  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 678
                                  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 692
                               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 714
 *  *
                                  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 767
                                     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 779
                                     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 803
                                     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 815
                                     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 854
                   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 865
 *     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 880
             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 888
          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 908
 *  *
       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.2  
changed lines
  Added in v.1.9

CVSweb interface <joel.bertrand@systella.fr>