[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.21

version 1.11, 2012/12/14 12:30:20 version 1.21, 2023/08/07 08:38:50
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  
 *  *
 *> \ingroup doubleGEcomputational  *> \ingroup doubleGEcomputational
 *  *
Line 279 Line 279
 *>  spectral condition number. The best performance of this Jacobi SVD  *>  spectral condition number. The best performance of this Jacobi SVD
 *>  procedure is achieved if used in an  accelerated version of Drmac and  *>  procedure is achieved if used in an  accelerated version of Drmac and
 *>  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].  *>  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
 *>  Some tunning parameters (marked with [TP]) are available for the  *>  Some tuning parameters (marked with [TP]) are available for the
 *>  implementer.  *>  implementer.
 *>  The computational range for the nonzero singular values is the  machine  *>  The computational range for the nonzero singular values is the  machine
 *>  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even  *>  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
Line 335 Line 335
       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 --
 *  -- 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  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N        INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
Line 374 Line 373
       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 427
          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 593
       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 635
       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 688
 *     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 700
       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 711
 *     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 809
 *  *
             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 862
 *  *
                      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 901
                               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 934
      $                                              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 950
                                     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 1067
      $                                       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 1135
                   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 1155
 *        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 1212
                               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 1240
      $                                              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 1255
      $                                  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 1375
                                     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 1393
                                     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 1465
                   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 1476
 *     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.21

CVSweb interface <joel.bertrand@systella.fr>