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

version 1.7, 2010/12/21 13:53:33 version 1.11, 2012/08/22 09:48:20
Line 1 Line 1
       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,  *> \brief \b DLASQ5
      $                   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 April 2012
   *
   *> \ingroup auxOTHERcomputational
   *
   *  =====================================================================
         SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
        $                   DN, DNM1, DNM2, IEEE, EPS )
   *
   *  -- LAPACK computational routine (version 3.4.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..--
   *     April 2012
 *  *
 *     .. 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
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
       RETURN        RETURN
Removed from v.1.7  
changed lines
  Added in v.1.11

CVSweb interface <joel.bertrand@systella.fr>