[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.11 and 1.20

version 1.11, 2012/12/14 12:30:20 version 1.20, 2020/05/21 21:45:57
Line 2 Line 2
 *  *
 *  =========== 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 DGESVJ + dependencies   *> Download DGESVJ + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       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 )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N  *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
 *       CHARACTER*1        JOBA, JOBU, JOBV  *       CHARACTER*1        JOBA, JOBU, JOBV
Line 29 Line 29
 *       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),  *       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
 *      $                   WORK( LWORK )  *      $                   WORK( LWORK )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 45 Line 45
 *> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements  *> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
 *> of SIGMA are the singular values of A. The columns of U and V are the  *> of SIGMA are the singular values of A. The columns of U and V are the
 *> left and the right singular vectors of A, respectively.  *> left and the right singular vectors of A, respectively.
   *> DGESVJ can sometimes compute tiny singular values and their singular vectors much
   *> more accurately than other SVD routines, see below under Further Details.
 *> \endverbatim  *> \endverbatim
 *  *
 *  Arguments:  *  Arguments:
Line 52 Line 54
 *  *
 *> \param[in] JOBA  *> \param[in] JOBA
 *> \verbatim  *> \verbatim
 *>          JOBA is CHARACTER* 1  *>          JOBA is CHARACTER*1
 *>          Specifies the structure of A.  *>          Specifies the structure of A.
 *>          = 'L': The input matrix A is lower triangular;  *>          = 'L': The input matrix A is lower triangular;
 *>          = 'U': The input matrix A is upper triangular;  *>          = 'U': The input matrix A is upper triangular;
Line 88 Line 90
 *>          JOBV is CHARACTER*1  *>          JOBV is CHARACTER*1
 *>          Specifies whether to compute the right singular vectors, that  *>          Specifies whether to compute the right singular vectors, that
 *>          is, the matrix V:  *>          is, the matrix V:
 *>          = 'V' : the matrix V is computed and returned in the array V  *>          = 'V':  the matrix V is computed and returned in the array V
 *>          = 'A' : the Jacobi rotations are applied to the MV-by-N  *>          = 'A':  the Jacobi rotations are applied to the MV-by-N
 *>                  array V. In other words, the right singular vector  *>                  array V. In other words, the right singular vector
 *>                  matrix V is not computed explicitly, instead it is  *>                  matrix V is not computed explicitly, instead it is
 *>                  applied to an MV-by-N matrix initially stored in the  *>                  applied to an MV-by-N matrix initially stored in the
 *>                  first MV rows of V.  *>                  first MV rows of V.
 *>          = 'N' : the matrix V is not computed and the array V is not  *>          = 'N':  the matrix V is not computed and the array V is not
 *>                  referenced  *>                  referenced
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] M  *> \param[in] M
 *> \verbatim  *> \verbatim
 *>          M is INTEGER  *>          M is INTEGER
 *>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.    *>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] N  *> \param[in] N
Line 116 Line 118
 *>          A is DOUBLE PRECISION array, dimension (LDA,N)  *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 *>          On entry, the M-by-N matrix A.  *>          On entry, the M-by-N matrix A.
 *>          On exit :  *>          On exit :
 *>          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :  *>          If JOBU = 'U' .OR. JOBU = 'C' :
 *>                 If INFO .EQ. 0 :  *>                 If INFO = 0 :
 *>                 RANKA orthonormal columns of U are returned in the  *>                 RANKA orthonormal columns of U are returned in the
 *>                 leading RANKA columns of the array A. Here RANKA <= N  *>                 leading RANKA columns of the array A. Here RANKA <= N
 *>                 is the number of computed singular values of A that are  *>                 is the number of computed singular values of A that are
Line 127 Line 129
 *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the  *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
 *>                 descriptions of SVA and WORK. The computed columns of U  *>                 descriptions of SVA and WORK. The computed columns of U
 *>                 are mutually numerically orthogonal up to approximately  *>                 are mutually numerically orthogonal up to approximately
 *>                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),  *>                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
 *>                 see the description of JOBU.  *>                 see the description of JOBU.
 *>                 If INFO .GT. 0 :  *>                 If INFO > 0 :
 *>                 the procedure DGESVJ did not converge in the given number  *>                 the procedure DGESVJ did not converge in the given number
 *>                 of iterations (sweeps). In that case, the computed  *>                 of iterations (sweeps). In that case, the computed
 *>                 columns of U may not be orthogonal up to TOL. The output  *>                 columns of U may not be orthogonal up to TOL. The output
Line 138 Line 140
 *>                 input matrix A in the sense that the residual  *>                 input matrix A in the sense that the residual
 *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.  *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
 *>  *>
 *>          If JOBU .EQ. 'N' :  *>          If JOBU = 'N' :
 *>                 If INFO .EQ. 0 :  *>                 If INFO = 0 :
 *>                 Note that the left singular vectors are 'for free' in the  *>                 Note that the left singular vectors are 'for free' in the
 *>                 one-sided Jacobi SVD algorithm. However, if only the  *>                 one-sided Jacobi SVD algorithm. However, if only the
 *>                 singular values are needed, the level of numerical  *>                 singular values are needed, the level of numerical
Line 148 Line 150
 *>                 numerically orthogonal up to approximately M*EPS. Thus,  *>                 numerically orthogonal up to approximately M*EPS. Thus,
 *>                 on exit, A contains the columns of U scaled with the  *>                 on exit, A contains the columns of U scaled with the
 *>                 corresponding singular values.  *>                 corresponding singular values.
 *>                 If INFO .GT. 0 :  *>                 If INFO > 0 :
 *>                 the procedure DGESVJ did not converge in the given number  *>                 the procedure DGESVJ did not converge in the given number
 *>                 of iterations (sweeps).  *>                 of iterations (sweeps).
 *> \endverbatim  *> \endverbatim
Line 163 Line 165
 *> \verbatim  *> \verbatim
 *>          SVA is DOUBLE PRECISION array, dimension (N)  *>          SVA is DOUBLE PRECISION array, dimension (N)
 *>          On exit :  *>          On exit :
 *>          If INFO .EQ. 0 :  *>          If INFO = 0 :
 *>          depending on the value SCALE = WORK(1), we have:  *>          depending on the value SCALE = WORK(1), we have:
 *>                 If SCALE .EQ. ONE :  *>                 If SCALE = ONE :
 *>                 SVA(1:N) contains the computed singular values of A.  *>                 SVA(1:N) contains the computed singular values of A.
 *>                 During the computation SVA contains the Euclidean column  *>                 During the computation SVA contains the Euclidean column
 *>                 norms of the iterated matrices in the array A.  *>                 norms of the iterated matrices in the array A.
Line 173 Line 175
 *>                 The singular values of A are SCALE*SVA(1:N), and this  *>                 The singular values of A are SCALE*SVA(1:N), and this
 *>                 factored representation is due to the fact that some of the  *>                 factored representation is due to the fact that some of the
 *>                 singular values of A might underflow or overflow.  *>                 singular values of A might underflow or overflow.
 *>          If INFO .GT. 0 :  *>          If INFO > 0 :
 *>          the procedure DGESVJ did not converge in the given number of  *>          the procedure DGESVJ did not converge in the given number of
 *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.  *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
 *> \endverbatim  *> \endverbatim
Line 181 Line 183
 *> \param[in] MV  *> \param[in] MV
 *> \verbatim  *> \verbatim
 *>          MV is INTEGER  *>          MV is INTEGER
 *>          If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ  *>          If JOBV = 'A', then the product of Jacobi rotations in DGESVJ
 *>          is applied to the first MV rows of V. See the description of JOBV.  *>          is applied to the first MV rows of V. See the description of JOBV.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 199 Line 201
 *> \param[in] LDV  *> \param[in] LDV
 *> \verbatim  *> \verbatim
 *>          LDV is INTEGER  *>          LDV is INTEGER
 *>          The leading dimension of the array V, LDV .GE. 1.  *>          The leading dimension of the array V, LDV >= 1.
 *>          If JOBV .EQ. 'V', then LDV .GE. max(1,N).  *>          If JOBV = 'V', then LDV >= max(1,N).
 *>          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .  *>          If JOBV = 'A', then LDV >= max(1,MV) .
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in,out] WORK  *> \param[in,out] WORK
 *> \verbatim  *> \verbatim
 *>          WORK is DOUBLE PRECISION array, dimension max(4,M+N).  *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 *>          On entry :  *>          On entry :
 *>          If JOBU .EQ. 'C' :  *>          If JOBU = 'C' :
 *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.  *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
 *>                    The process stops if all columns of A are mutually  *>                    The process stops if all columns of A are mutually
 *>                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').  *>                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
Line 228 Line 230
 *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.  *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
 *>                    This is useful information in cases when DGESVJ did  *>                    This is useful information in cases when DGESVJ did
 *>                    not converge, as it can be used to estimate whether  *>                    not converge, as it can be used to estimate whether
 *>                    the output is stil useful and for post festum analysis.  *>                    the output is still useful and for post festum analysis.
 *>          WORK(6) = the largest absolute value over all sines of the  *>          WORK(6) = the largest absolute value over all sines of the
 *>                    Jacobi rotation angles in the last sweep. It can be  *>                    Jacobi rotation angles in the last sweep. It can be
 *>                    useful for a post festum analysis.  *>                    useful for a post festum analysis.
Line 243 Line 245
 *> \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
 *>          > 0 : DGESVJ did not converge in the maximal allowed number (30)  *>          > 0:  DGESVJ did not converge in the maximal allowed number (30)
 *>                of sweeps. The output may still be useful. See the  *>                of sweeps. The output may still be useful. See the
 *>                description of WORK.  *>                description of WORK.
 *> \endverbatim  *> \endverbatim
Line 253 Line 255
 *  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 September 2012  *> \date June 2017
 *  *
 *> \ingroup doubleGEcomputational  *> \ingroup doubleGEcomputational
 *  *
Line 335 Line 337
       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 computational routine (version 3.4.2) --  *  -- LAPACK computational routine (version 3.7.1) --
 *  -- 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..--
 *     September 2012  *     June 2017
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N        INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
Line 374 Line 376
       DOUBLE PRECISION   FASTR( 5 )        DOUBLE PRECISION   FASTR( 5 )
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT        INTRINSIC          DABS, MAX, MIN, DBLE, DSIGN, DSQRT
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
 *     ..  *     ..
Line 428 Line 430
          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
       ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN        ELSE IF( LWORK.LT.MAX( M+N, 6 ) ) THEN
          INFO = -13           INFO = -13
       ELSE        ELSE
          INFO = 0           INFO = 0
Line 594 Line 596
       AAPP = ZERO        AAPP = ZERO
       AAQQ = BIG        AAQQ = BIG
       DO 4781 p = 1, N        DO 4781 p = 1, N
          IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) )           IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) )
          AAPP = DMAX1( AAPP, SVA( p ) )           AAPP = MAX( AAPP, SVA( p ) )
  4781 CONTINUE   4781 CONTINUE
 *  *
 * #:) Quick return for zero matrix  * #:) Quick return for zero matrix
Line 636 Line 638
       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 = MIN( BIG, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1  *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1  *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN        ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
          TEMP1 = DMIN1( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) )           TEMP1 = MIN( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) )
 *         AAQQ  = AAQQ*TEMP1  *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1  *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN        ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
          TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP )           TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1  *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1  *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN        ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
          TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) )           TEMP1 = MIN( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) )
 *         AAQQ  = AAQQ*TEMP1  *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1  *         AAPP  = AAPP*TEMP1
       ELSE        ELSE
Line 689 Line 691
 *     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 701 Line 703
       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 712 Line 714
 *     invokes cubic convergence. Big part of this cycle is done inside  *     invokes cubic convergence. Big part of this cycle is done inside
 *     canonical subspaces of dimensions less than M.  *     canonical subspaces of dimensions less than M.
 *  *
       IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN        IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN
 *[TP] The number of partition levels and the actual partition are  *[TP] The number of partition levels and the actual partition are
 *     tuning parameters.  *     tuning parameters.
          N4 = N / 4           N4 = N / 4
Line 810 Line 812
 *  *
             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
 *  *
Line 863 Line 865
 *  *
                      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 902 Line 904
                               END IF                                END IF
                            END IF                             END IF
 *  *
                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )                             MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) )
 *  *
 *        TO rotate or NOT to rotate, THAT is the question ...  *        TO rotate or NOT to rotate, THAT is the question ...
 *  *
Line 935 Line 937
      $                                              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( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = MAX( MXSINJ, DABS( T ) )
 *  *
                                  ELSE                                   ELSE
 *  *
Line 951 Line 953
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = DSQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
 *  *
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )                                      MXSINJ = MAX( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
 *  *
                                     APOAQ = WORK( p ) / WORK( q )                                      APOAQ = WORK( p ) / WORK( q )
Line 1068 Line 1070
      $                                       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( MAX( ZERO,
      $                                      ONE-AAPQ*AAPQ ) )       $                                      ONE-AAPQ*AAPQ ) )
                                  MXSINJ = DMAX1( MXSINJ, SFMIN )                                   MXSINJ = MAX( MXSINJ, SFMIN )
                               END IF                                END IF
 *           END IF ROTOK THEN ... ELSE  *           END IF ROTOK THEN ... ELSE
 *  *
Line 1136 Line 1138
                   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 1156 Line 1158
 *        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 1213 Line 1215
                               END IF                                END IF
                            END IF                             END IF
 *  *
                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )                             MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) )
 *  *
 *        TO rotate or NOT to rotate, THAT is the question ...  *        TO rotate or NOT to rotate, THAT is the question ...
 *  *
Line 1241 Line 1243
      $                                              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( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = MAX( MXSINJ, DABS( T ) )
                                  ELSE                                   ELSE
 *  *
 *                 .. choose correct signum for THETA and rotate  *                 .. choose correct signum for THETA and rotate
Line 1256 Line 1258
      $                                  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 = MAX( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                       AAPP = AAPP*DSQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
 *  *
                                     APOAQ = WORK( p ) / WORK( q )                                      APOAQ = WORK( p ) / WORK( q )
Line 1376 Line 1378
                                     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( MAX( ZERO,
      $                                         ONE-AAPQ*AAPQ ) )       $                                         ONE-AAPQ*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = MAX( 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 )
Line 1394 Line 1396
                                     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( MAX( ZERO,
      $                                         ONE-AAPQ*AAPQ ) )       $                                         ONE-AAPQ*AAPQ ) )
                                     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 1466 Line 1468
                   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 1477 Line 1479
 *     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 ) = DABS( SVA( p ) )                 SVA( p ) = DABS( SVA( p ) )
  2012       CONTINUE   2012       CONTINUE
 ***  ***
Removed from v.1.11  
changed lines
  Added in v.1.20

CVSweb interface <joel.bertrand@systella.fr>