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

Diff for /rpl/lapack/lapack/dlasq5.f between versions 1.7 and 1.18

version 1.7, 2010/12/21 13:53:33 version 1.18, 2018/05/29 06:55:20
Line 1 Line 1
       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,  *> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
      $                   DNM1, DNM2, IEEE )  
 *  *
 *  -- LAPACK routine (version 3.2)                                    --  *  =========== DOCUMENTATION ===========
 *  *
 *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --  * Online html documentation available at
 *  -- Laboratory and Beresford Parlett of the Univ. of California at  --  *            http://www.netlib.org/lapack/explore-html/
 *  -- Berkeley                                                        --  
 *  -- November 2008                                                   --  
 *  *
   *> \htmlonly
   *> Download DLASQ5 + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
   *                          DNM1, DNM2, IEEE, EPS )
   *
   *       .. Scalar Arguments ..
   *       LOGICAL            IEEE
   *       INTEGER            I0, N0, PP
   *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   Z( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DLASQ5 computes one dqds transform in ping-pong form, one
   *> version for IEEE machines another for non IEEE machines.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] I0
   *> \verbatim
   *>          I0 is INTEGER
   *>        First index.
   *> \endverbatim
   *>
   *> \param[in] N0
   *> \verbatim
   *>          N0 is INTEGER
   *>        Last index.
   *> \endverbatim
   *>
   *> \param[in] Z
   *> \verbatim
   *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
   *>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
   *>        an extra argument.
   *> \endverbatim
   *>
   *> \param[in] PP
   *> \verbatim
   *>          PP is INTEGER
   *>        PP=0 for ping, PP=1 for pong.
   *> \endverbatim
   *>
   *> \param[in] TAU
   *> \verbatim
   *>          TAU is DOUBLE PRECISION
   *>        This is the shift.
   *> \endverbatim
   *>
   *> \param[in] SIGMA
   *> \verbatim
   *>          SIGMA is DOUBLE PRECISION
   *>        This is the accumulated shift up to this step.
   *> \endverbatim
   *>
   *> \param[out] DMIN
   *> \verbatim
   *>          DMIN is DOUBLE PRECISION
   *>        Minimum value of d.
   *> \endverbatim
   *>
   *> \param[out] DMIN1
   *> \verbatim
   *>          DMIN1 is DOUBLE PRECISION
   *>        Minimum value of d, excluding D( N0 ).
   *> \endverbatim
   *>
   *> \param[out] DMIN2
   *> \verbatim
   *>          DMIN2 is DOUBLE PRECISION
   *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
   *> \endverbatim
   *>
   *> \param[out] DN
   *> \verbatim
   *>          DN is DOUBLE PRECISION
   *>        d(N0), the last value of d.
   *> \endverbatim
   *>
   *> \param[out] DNM1
   *> \verbatim
   *>          DNM1 is DOUBLE PRECISION
   *>        d(N0-1).
   *> \endverbatim
   *>
   *> \param[out] DNM2
   *> \verbatim
   *>          DNM2 is DOUBLE PRECISION
   *>        d(N0-2).
   *> \endverbatim
   *>
   *> \param[in] IEEE
   *> \verbatim
   *>          IEEE is LOGICAL
   *>        Flag for IEEE or non IEEE arithmetic.
   *> \endverbatim
   *>
   *> \param[in] EPS
   *> \verbatim
   *>          EPS is DOUBLE PRECISION
   *>        This is the value of epsilon used.
   *> \endverbatim
   *>
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \date June 2017
   *
   *> \ingroup auxOTHERcomputational
   *
   *  =====================================================================
         SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
        $                   DN, DNM1, DNM2, IEEE, EPS )
   *
   *  -- 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..--
   *     June 2017
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       LOGICAL            IEEE        LOGICAL            IEEE
       INTEGER            I0, N0, PP        INTEGER            I0, N0, PP
       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU        DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
        $                   SIGMA, EPS
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       DOUBLE PRECISION   Z( * )        DOUBLE PRECISION   Z( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DLASQ5 computes one dqds transform in ping-pong form, one  
 *  version for IEEE machines another for non IEEE machines.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  I0    (input) INTEGER  
 *        First index.  
 *  
 *  N0    (input) INTEGER  
 *        Last index.  
 *  
 *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )  
 *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid  
 *        an extra argument.  
 *  
 *  PP    (input) INTEGER  
 *        PP=0 for ping, PP=1 for pong.  
 *  
 *  TAU   (input) DOUBLE PRECISION  
 *        This is the shift.  
 *  
 *  DMIN  (output) DOUBLE PRECISION  
 *        Minimum value of d.  
 *  
 *  DMIN1 (output) DOUBLE PRECISION  
 *        Minimum value of d, excluding D( N0 ).  
 *  
 *  DMIN2 (output) DOUBLE PRECISION  
 *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).  
 *  
 *  DN    (output) DOUBLE PRECISION  
 *        d(N0), the last value of d.  
 *  
 *  DNM1  (output) DOUBLE PRECISION  
 *        d(N0-1).  
 *  
 *  DNM2  (output) DOUBLE PRECISION  
 *        d(N0-2).  
 *  
 *  IEEE  (input) LOGICAL  
 *        Flag for IEEE or non IEEE arithmetic.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameter ..  *     .. Parameter ..
       DOUBLE PRECISION   ZERO        DOUBLE PRECISION   ZERO, HALF
       PARAMETER          ( ZERO = 0.0D0 )        PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       INTEGER            J4, J4P2        INTEGER            J4, J4P2
       DOUBLE PRECISION   D, EMIN, TEMP        DOUBLE PRECISION   D, EMIN, TEMP, DTHRESH
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          MIN        INTRINSIC          MIN
Line 84 Line 177
       IF( ( N0-I0-1 ).LE.0 )        IF( ( N0-I0-1 ).LE.0 )
      $   RETURN       $   RETURN
 *  *
         DTHRESH = EPS*(SIGMA+TAU)
         IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
         IF( TAU.NE.ZERO ) THEN
       J4 = 4*I0 + PP - 3        J4 = 4*I0 + PP - 3
       EMIN = Z( J4+4 )         EMIN = Z( J4+4 )
       D = Z( J4 ) - TAU        D = Z( J4 ) - TAU
       DMIN = D        DMIN = D
       DMIN1 = -Z( J4 )        DMIN1 = -Z( J4 )
Line 96 Line 192
 *  *
          IF( PP.EQ.0 ) THEN           IF( PP.EQ.0 ) THEN
             DO 10 J4 = 4*I0, 4*( N0-3 ), 4              DO 10 J4 = 4*I0, 4*( N0-3 ), 4
                Z( J4-2 ) = D + Z( J4-1 )                  Z( J4-2 ) = D + Z( J4-1 )
                TEMP = Z( J4+1 ) / Z( J4-2 )                 TEMP = Z( J4+1 ) / Z( J4-2 )
                D = D*TEMP - TAU                 D = D*TEMP - TAU
                DMIN = MIN( DMIN, D )                 DMIN = MIN( DMIN, D )
Line 105 Line 201
    10       CONTINUE     10       CONTINUE
          ELSE           ELSE
             DO 20 J4 = 4*I0, 4*( N0-3 ), 4              DO 20 J4 = 4*I0, 4*( N0-3 ), 4
                Z( J4-3 ) = D + Z( J4 )                  Z( J4-3 ) = D + Z( J4 )
                TEMP = Z( J4+2 ) / Z( J4-3 )                 TEMP = Z( J4+2 ) / Z( J4-3 )
                D = D*TEMP - TAU                 D = D*TEMP - TAU
                DMIN = MIN( DMIN, D )                 DMIN = MIN( DMIN, D )
Line 114 Line 210
    20       CONTINUE     20       CONTINUE
          END IF           END IF
 *  *
 *        Unroll last two steps.   *        Unroll last two steps.
 *  *
          DNM2 = D           DNM2 = D
          DMIN2 = DMIN           DMIN2 = DMIN
Line 139 Line 235
 *  *
          IF( PP.EQ.0 ) THEN           IF( PP.EQ.0 ) THEN
             DO 30 J4 = 4*I0, 4*( N0-3 ), 4              DO 30 J4 = 4*I0, 4*( N0-3 ), 4
                Z( J4-2 ) = D + Z( J4-1 )                  Z( J4-2 ) = D + Z( J4-1 )
                IF( D.LT.ZERO ) THEN                 IF( D.LT.ZERO ) THEN
                   RETURN                    RETURN
                ELSE                  ELSE
                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )                    Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU                    D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
                END IF                 END IF
Line 151 Line 247
    30       CONTINUE     30       CONTINUE
          ELSE           ELSE
             DO 40 J4 = 4*I0, 4*( N0-3 ), 4              DO 40 J4 = 4*I0, 4*( N0-3 ), 4
                Z( J4-3 ) = D + Z( J4 )                  Z( J4-3 ) = D + Z( J4 )
                IF( D.LT.ZERO ) THEN                 IF( D.LT.ZERO ) THEN
                   RETURN                    RETURN
                ELSE                  ELSE
                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )                    Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU                    D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
                END IF                 END IF
Line 163 Line 259
    40       CONTINUE     40       CONTINUE
          END IF           END IF
 *  *
 *        Unroll last two steps.   *        Unroll last two steps.
 *  *
          DNM2 = D           DNM2 = D
          DMIN2 = DMIN           DMIN2 = DMIN
Line 191 Line 287
          DMIN = MIN( DMIN, DN )           DMIN = MIN( DMIN, DN )
 *  *
       END IF        END IF
         ELSE
   *     This is the version that sets d's to zero if they are small enough
            J4 = 4*I0 + PP - 3
            EMIN = Z( J4+4 )
            D = Z( J4 ) - TAU
            DMIN = D
            DMIN1 = -Z( J4 )
            IF( IEEE ) THEN
   *
   *     Code for IEEE arithmetic.
   *
               IF( PP.EQ.0 ) THEN
                  DO 50 J4 = 4*I0, 4*( N0-3 ), 4
                     Z( J4-2 ) = D + Z( J4-1 )
                     TEMP = Z( J4+1 ) / Z( J4-2 )
                     D = D*TEMP - TAU
                     IF( D.LT.DTHRESH ) D = ZERO
                     DMIN = MIN( DMIN, D )
                     Z( J4 ) = Z( J4-1 )*TEMP
                     EMIN = MIN( Z( J4 ), EMIN )
    50            CONTINUE
               ELSE
                  DO 60 J4 = 4*I0, 4*( N0-3 ), 4
                     Z( J4-3 ) = D + Z( J4 )
                     TEMP = Z( J4+2 ) / Z( J4-3 )
                     D = D*TEMP - TAU
                     IF( D.LT.DTHRESH ) D = ZERO
                     DMIN = MIN( DMIN, D )
                     Z( J4-1 ) = Z( J4 )*TEMP
                     EMIN = MIN( Z( J4-1 ), EMIN )
    60            CONTINUE
               END IF
   *
   *     Unroll last two steps.
   *
               DNM2 = D
               DMIN2 = DMIN
               J4 = 4*( N0-2 ) - PP
               J4P2 = J4 + 2*PP - 1
               Z( J4-2 ) = DNM2 + Z( J4P2 )
               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
               DMIN = MIN( DMIN, DNM1 )
   *
               DMIN1 = DMIN
               J4 = J4 + 4
               J4P2 = J4 + 2*PP - 1
               Z( J4-2 ) = DNM1 + Z( J4P2 )
               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
               DMIN = MIN( DMIN, DN )
   *
            ELSE
   *
   *     Code for non IEEE arithmetic.
   *
               IF( PP.EQ.0 ) THEN
                  DO 70 J4 = 4*I0, 4*( N0-3 ), 4
                     Z( J4-2 ) = D + Z( J4-1 )
                     IF( D.LT.ZERO ) THEN
                        RETURN
                     ELSE
                        Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
                        D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
                     END IF
                     IF( D.LT.DTHRESH) D = ZERO
                     DMIN = MIN( DMIN, D )
                     EMIN = MIN( EMIN, Z( J4 ) )
    70            CONTINUE
               ELSE
                  DO 80 J4 = 4*I0, 4*( N0-3 ), 4
                     Z( J4-3 ) = D + Z( J4 )
                     IF( D.LT.ZERO ) THEN
                        RETURN
                     ELSE
                        Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
                        D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
                     END IF
                     IF( D.LT.DTHRESH) D = ZERO
                     DMIN = MIN( DMIN, D )
                     EMIN = MIN( EMIN, Z( J4-1 ) )
    80            CONTINUE
               END IF
   *
   *     Unroll last two steps.
   *
               DNM2 = D
               DMIN2 = DMIN
               J4 = 4*( N0-2 ) - PP
               J4P2 = J4 + 2*PP - 1
               Z( J4-2 ) = DNM2 + Z( J4P2 )
               IF( DNM2.LT.ZERO ) THEN
                  RETURN
               ELSE
                  Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
                  DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
               END IF
               DMIN = MIN( DMIN, DNM1 )
   *
               DMIN1 = DMIN
               J4 = J4 + 4
               J4P2 = J4 + 2*PP - 1
               Z( J4-2 ) = DNM1 + Z( J4P2 )
               IF( DNM1.LT.ZERO ) THEN
                  RETURN
               ELSE
                  Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
                  DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
               END IF
               DMIN = MIN( DMIN, DN )
   *
            END IF
         END IF
 *  *
       Z( J4+2 ) = DN        Z( J4+2 ) = DN
       Z( 4*N0-PP ) = EMIN        Z( 4*N0-PP ) = EMIN
Removed from v.1.7  
changed lines
  Added in v.1.18

CVSweb interface <joel.bertrand@systella.fr>