--- rpl/lapack/lapack/dlasd4.f	2012/12/14 12:30:25	1.13
+++ rpl/lapack/lapack/dlasd4.f	2017/06/17 10:53:57	1.18
@@ -140,9 +140,9 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \date September 2012
+*> \date December 2016
 *
-*> \ingroup auxOTHERauxiliary
+*> \ingroup OTHERauxiliary
 *
 *> \par Contributors:
 *  ==================
@@ -153,10 +153,10 @@
 *  =====================================================================
       SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
 *
-*  -- LAPACK auxiliary routine (version 3.4.2) --
+*  -- LAPACK auxiliary routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     September 2012
+*     December 2016
 *
 *     .. Scalar Arguments ..
       INTEGER            I, INFO, N
@@ -223,6 +223,7 @@
 *
       EPS = DLAMCH( 'Epsilon' )
       RHOINV = ONE / RHO
+      TAU2= ZERO
 *
 *     The case I = N
 *
@@ -275,6 +276,7 @@
                ELSE
                   TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
                END IF
+               TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
             END IF
 *
 *           It can be proved that
@@ -293,6 +295,8 @@
             ELSE
                TAU2 = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
             END IF
+            TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
+
 *
 *           It can be proved that
 *           D(N)^2 < D(N)^2+TAU2 < SIGMA(N)^2 < D(N)^2+RHO/2
@@ -301,7 +305,7 @@
 *
 *        The following TAU is to approximate SIGMA_n - D( N )
 *
-         TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
+*         TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) )
 *
          SIGMA = D( N ) + TAU
          DO 30 J = 1, N
@@ -327,7 +331,7 @@
          TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )
          PHI = Z( N )*TEMP
          DPHI = TEMP*TEMP
-         ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV  
+         ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
 *    $          + ABS( TAU2 )*( DPSI+DPHI )
 *
          W = RHOINV + PHI + PSI
@@ -396,7 +400,7 @@
          TEMP = Z( N ) / TAU2
          PHI = Z( N )*TEMP
          DPHI = TEMP*TEMP
-         ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV  
+         ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
 *    $          + ABS( TAU2 )*( DPSI+DPHI )
 *
          W = RHOINV + PHI + PSI
@@ -466,7 +470,7 @@
             TEMP = Z( N ) / TAU2
             PHI = Z( N )*TEMP
             DPHI = TEMP*TEMP
-            ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV  
+            ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV
 *    $             + ABS( TAU2 )*( DPSI+DPHI )
 *
             W = RHOINV + PHI + PSI
@@ -618,8 +622,8 @@
          DW = DPSI + DPHI + TEMP*TEMP
          TEMP = Z( II )*TEMP
          W = W + TEMP
-         ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV 
-     $          + THREE*ABS( TEMP ) 
+         ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
+     $          + THREE*ABS( TEMP )
 *    $          + ABS( TAU2 )*DW
 *
 *        Test for convergence
@@ -699,7 +703,7 @@
 *
             IF( INFO.NE.0 ) THEN
 *
-*              If INFO is not 0, i.e., DLAED6 failed, switch back 
+*              If INFO is not 0, i.e., DLAED6 failed, switch back
 *              to 2 pole interpolation.
 *
                SWTCH3 = .FALSE.
@@ -799,8 +803,8 @@
          DW = DPSI + DPHI + TEMP*TEMP
          TEMP = Z( II )*TEMP
          W = RHOINV + PHI + PSI + TEMP
-         ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV 
-     $          + THREE*ABS( TEMP ) 
+         ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
+     $          + THREE*ABS( TEMP )
 *    $          + ABS( TAU2 )*DW
 *
          SWTCH = .FALSE.
@@ -918,7 +922,7 @@
 *
                IF( INFO.NE.0 ) THEN
 *
-*                 If INFO is not 0, i.e., DLAED6 failed, switch 
+*                 If INFO is not 0, i.e., DLAED6 failed, switch
 *                 back to two pole interpolation
 *
                   SWTCH3 = .FALSE.
@@ -1034,8 +1038,8 @@
             DW = DPSI + DPHI + TEMP*TEMP
             TEMP = Z( II )*TEMP
             W = RHOINV + PHI + PSI + TEMP
-            ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV 
-     $             + THREE*ABS( TEMP ) 
+            ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV
+     $             + THREE*ABS( TEMP )
 *    $             + ABS( TAU2 )*DW
 *
             IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )