--- rpl/lapack/lapack/dlasq5.f	2011/11/21 22:19:35	1.9
+++ rpl/lapack/lapack/dlasq5.f	2018/05/29 06:55:20	1.18
@@ -1,35 +1,35 @@
-*> \brief \b DLASQ5
+*> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
 *
 *  =========== DOCUMENTATION ===========
 *
-* Online html documentation available at 
-*            http://www.netlib.org/lapack/explore-html/ 
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
 *
 *> \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"> 
+*> 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 
+*> \endhtmlonly
 *
 *  Definition:
 *  ===========
 *
-*       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-*                          DNM1, DNM2, IEEE )
-* 
+*       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
+*       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
 *       ..
 *       .. Array Arguments ..
 *       DOUBLE PRECISION   Z( * )
 *       ..
-*  
+*
 *
 *> \par Purpose:
 *  =============
@@ -74,6 +74,12 @@
 *>        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
@@ -115,32 +121,39 @@
 *>          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. 
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
 *
-*> \date November 2011
+*> \date June 2017
 *
 *> \ingroup auxOTHERcomputational
 *
 *  =====================================================================
-      SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
-     $                   DNM1, DNM2, IEEE )
+      SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
+     $                   DN, DNM1, DNM2, IEEE, EPS )
 *
-*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK computational routine (version 3.7.1) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2011
+*     June 2017
 *
 *     .. 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( * )
@@ -149,12 +162,12 @@
 *  =====================================================================
 *
 *     .. 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
@@ -164,8 +177,11 @@
       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 ) 
+      EMIN = Z( J4+4 )
       D = Z( J4 ) - TAU
       DMIN = D
       DMIN1 = -Z( J4 )
@@ -176,7 +192,7 @@
 *
          IF( PP.EQ.0 ) THEN
             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 )
                D = D*TEMP - TAU
                DMIN = MIN( DMIN, D )
@@ -185,7 +201,7 @@
    10       CONTINUE
          ELSE
             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 )
                D = D*TEMP - TAU
                DMIN = MIN( DMIN, D )
@@ -194,7 +210,7 @@
    20       CONTINUE
          END IF
 *
-*        Unroll last two steps. 
+*        Unroll last two steps.
 *
          DNM2 = D
          DMIN2 = DMIN
@@ -219,10 +235,10 @@
 *
          IF( PP.EQ.0 ) THEN
             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
                   RETURN
-               ELSE 
+               ELSE
                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
                END IF
@@ -231,10 +247,10 @@
    30       CONTINUE
          ELSE
             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
                   RETURN
-               ELSE 
+               ELSE
                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
                END IF
@@ -243,7 +259,7 @@
    40       CONTINUE
          END IF
 *
-*        Unroll last two steps. 
+*        Unroll last two steps.
 *
          DNM2 = D
          DMIN2 = DMIN
@@ -271,6 +287,119 @@
          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