--- rpl/lapack/lapack/dlasq5.f 2010/12/21 13:53:33 1.7
+++ rpl/lapack/lapack/dlasq5.f 2014/01/27 09:28:23 1.14
@@ -1,80 +1,173 @@
- SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
- $ DNM1, DNM2, IEEE )
+*> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
*
-* -- LAPACK routine (version 3.2) --
+* =========== DOCUMENTATION ===========
*
-* -- Contributed by Osni Marques of the Lawrence Berkeley National --
-* -- Laboratory and Beresford Parlett of the Univ. of California at --
-* -- Berkeley --
-* -- November 2008 --
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
+*> \htmlonly
+*> Download DLASQ5 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \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 September 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.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* September 2012
*
* .. Scalar Arguments ..
LOGICAL IEEE
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 ..
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 ..
- DOUBLE PRECISION ZERO
- PARAMETER ( ZERO = 0.0D0 )
+ DOUBLE PRECISION ZERO, HALF
+ PARAMETER ( ZERO = 0.0D0, HALF = 0.5 )
* ..
* .. Local Scalars ..
INTEGER J4, J4P2
- DOUBLE PRECISION D, EMIN, TEMP
+ DOUBLE PRECISION D, EMIN, TEMP, DTHRESH
* ..
* .. Intrinsic Functions ..
INTRINSIC MIN
@@ -84,6 +177,9 @@
IF( ( N0-I0-1 ).LE.0 )
$ RETURN
*
+ DTHRESH = EPS*(SIGMA+TAU)
+ IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
+ IF( TAU.NE.ZERO ) THEN
J4 = 4*I0 + PP - 3
EMIN = Z( J4+4 )
D = Z( J4 ) - TAU
@@ -191,7 +287,120 @@
DMIN = MIN( DMIN, DN )
*
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( 4*N0-PP ) = EMIN
RETURN