Moving lapack to right place

1 files changed, 0 insertions, 448 deletions

diff --git a/src/lib/lapack/dlasq2.f b/src/lib/lapack/dlasq2.f
deleted file mode 100644
index b6b79aeb..00000000
--- a/src/lib/lapack/dlasq2.f
+++ /dev/null
@@ -1,448 +0,0 @@
-      SUBROUTINE DLASQ2( N, Z, INFO )
-*
-*  -- LAPACK routine (version 3.1) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     Modified to call DLAZQ3 in place of DLASQ3, 13 Feb 03, SJH.
-*
-*     .. Scalar Arguments ..
-      INTEGER            INFO, N
-*     ..
-*     .. Array Arguments ..
-      DOUBLE PRECISION   Z( * )
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DLASQ2 computes all the eigenvalues of the symmetric positive 
-*  definite tridiagonal matrix associated with the qd array Z to high
-*  relative accuracy are computed to high relative accuracy, in the
-*  absence of denormalization, underflow and overflow.
-*
-*  To see the relation of Z to the tridiagonal matrix, let L be a
-*  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
-*  let U be an upper bidiagonal matrix with 1's above and diagonal
-*  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
-*  symmetric tridiagonal to which it is similar.
-*
-*  Note : DLASQ2 defines a logical variable, IEEE, which is true
-*  on machines which follow ieee-754 floating-point standard in their
-*  handling of infinities and NaNs, and false otherwise. This variable
-*  is passed to DLAZQ3.
-*
-*  Arguments
-*  =========
-*
-*  N     (input) INTEGER
-*        The number of rows and columns in the matrix. N >= 0.
-*
-*  Z     (workspace) DOUBLE PRECISION array, dimension ( 4*N )
-*        On entry Z holds the qd array. On exit, entries 1 to N hold
-*        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
-*        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
-*        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
-*        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
-*        shifts that failed.
-*
-*  INFO  (output) INTEGER
-*        = 0: successful exit
-*        < 0: if the i-th argument is a scalar and had an illegal
-*             value, then INFO = -i, if the i-th argument is an
-*             array and the j-entry had an illegal value, then
-*             INFO = -(i*100+j)
-*        > 0: the algorithm failed
-*              = 1, a split was marked by a positive value in E
-*              = 2, current block of Z not diagonalized after 30*N
-*                   iterations (in inner while loop)
-*              = 3, termination criterion of outer while loop not met 
-*                   (program created more than N unreduced blocks)
-*
-*  Further Details
-*  ===============
-*  Local Variables: I0:N0 defines a current unreduced segment of Z.
-*  The shifts are accumulated in SIGMA. Iteration count is in ITER.
-*  Ping-pong is controlled by PP (alternates between 0 and 1).
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   CBIAS
-      PARAMETER          ( CBIAS = 1.50D0 )
-      DOUBLE PRECISION   ZERO, HALF, ONE, TWO, FOUR, HUNDRD
-      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
-     $                     TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 )
-*     ..
-*     .. Local Scalars ..
-      LOGICAL            IEEE
-      INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, 
-     $                   N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
-      DOUBLE PRECISION   D, DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, E,
-     $                   EMAX, EMIN, EPS, OLDEMN, QMAX, QMIN, S, SAFMIN,
-     $                   SIGMA, T, TAU, TEMP, TOL, TOL2, TRACE, ZMAX
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL           DLAZQ3, DLASRT, XERBLA
-*     ..
-*     .. External Functions ..
-      INTEGER            ILAENV
-      DOUBLE PRECISION   DLAMCH
-      EXTERNAL           DLAMCH, ILAENV
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
-*     ..
-*     .. Executable Statements ..
-*      
-*     Test the input arguments.
-*     (in case DLASQ2 is not called by DLASQ1)
-*
-      INFO = 0
-      EPS = DLAMCH( 'Precision' )
-      SAFMIN = DLAMCH( 'Safe minimum' )
-      TOL = EPS*HUNDRD
-      TOL2 = TOL**2
-*
-      IF( N.LT.0 ) THEN
-         INFO = -1
-         CALL XERBLA( 'DLASQ2', 1 )
-         RETURN
-      ELSE IF( N.EQ.0 ) THEN
-         RETURN
-      ELSE IF( N.EQ.1 ) THEN
-*
-*        1-by-1 case.
-*
-         IF( Z( 1 ).LT.ZERO ) THEN
-            INFO = -201
-            CALL XERBLA( 'DLASQ2', 2 )
-         END IF
-         RETURN
-      ELSE IF( N.EQ.2 ) THEN
-*
-*        2-by-2 case.
-*
-         IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
-            INFO = -2
-            CALL XERBLA( 'DLASQ2', 2 )
-            RETURN
-         ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
-            D = Z( 3 )
-            Z( 3 ) = Z( 1 )
-            Z( 1 ) = D
-         END IF
-         Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
-         IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
-            T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 
-            S = Z( 3 )*( Z( 2 ) / T )
-            IF( S.LE.T ) THEN
-               S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
-            ELSE
-               S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
-            END IF
-            T = Z( 1 ) + ( S+Z( 2 ) )
-            Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
-            Z( 1 ) = T
-         END IF
-         Z( 2 ) = Z( 3 )
-         Z( 6 ) = Z( 2 ) + Z( 1 )
-         RETURN
-      END IF
-*
-*     Check for negative data and compute sums of q's and e's.
-*
-      Z( 2*N ) = ZERO
-      EMIN = Z( 2 )
-      QMAX = ZERO
-      ZMAX = ZERO
-      D = ZERO
-      E = ZERO
-*
-      DO 10 K = 1, 2*( N-1 ), 2
-         IF( Z( K ).LT.ZERO ) THEN
-            INFO = -( 200+K )
-            CALL XERBLA( 'DLASQ2', 2 )
-            RETURN
-         ELSE IF( Z( K+1 ).LT.ZERO ) THEN
-            INFO = -( 200+K+1 )
-            CALL XERBLA( 'DLASQ2', 2 )
-            RETURN
-         END IF
-         D = D + Z( K )
-         E = E + Z( K+1 )
-         QMAX = MAX( QMAX, Z( K ) )
-         EMIN = MIN( EMIN, Z( K+1 ) )
-         ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
-   10 CONTINUE
-      IF( Z( 2*N-1 ).LT.ZERO ) THEN
-         INFO = -( 200+2*N-1 )
-         CALL XERBLA( 'DLASQ2', 2 )
-         RETURN
-      END IF
-      D = D + Z( 2*N-1 )
-      QMAX = MAX( QMAX, Z( 2*N-1 ) )
-      ZMAX = MAX( QMAX, ZMAX )
-*
-*     Check for diagonality.
-*
-      IF( E.EQ.ZERO ) THEN
-         DO 20 K = 2, N
-            Z( K ) = Z( 2*K-1 )
-   20    CONTINUE
-         CALL DLASRT( 'D', N, Z, IINFO )
-         Z( 2*N-1 ) = D
-         RETURN
-      END IF
-*
-      TRACE = D + E
-*
-*     Check for zero data.
-*
-      IF( TRACE.EQ.ZERO ) THEN
-         Z( 2*N-1 ) = ZERO
-         RETURN
-      END IF
-*         
-*     Check whether the machine is IEEE conformable.
-*         
-      IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
-     $       ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1      
-*         
-*     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
-*
-      DO 30 K = 2*N, 2, -2
-         Z( 2*K ) = ZERO 
-         Z( 2*K-1 ) = Z( K ) 
-         Z( 2*K-2 ) = ZERO 
-         Z( 2*K-3 ) = Z( K-1 ) 
-   30 CONTINUE
-*
-      I0 = 1
-      N0 = N
-*
-*     Reverse the qd-array, if warranted.
-*
-      IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
-         IPN4 = 4*( I0+N0 )
-         DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
-            TEMP = Z( I4-3 )
-            Z( I4-3 ) = Z( IPN4-I4-3 )
-            Z( IPN4-I4-3 ) = TEMP
-            TEMP = Z( I4-1 )
-            Z( I4-1 ) = Z( IPN4-I4-5 )
-            Z( IPN4-I4-5 ) = TEMP
-   40    CONTINUE
-      END IF
-*
-*     Initial split checking via dqd and Li's test.
-*
-      PP = 0
-*
-      DO 80 K = 1, 2
-*
-         D = Z( 4*N0+PP-3 )
-         DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
-            IF( Z( I4-1 ).LE.TOL2*D ) THEN
-               Z( I4-1 ) = -ZERO
-               D = Z( I4-3 )
-            ELSE
-               D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
-            END IF
-   50    CONTINUE
-*
-*        dqd maps Z to ZZ plus Li's test.
-*
-         EMIN = Z( 4*I0+PP+1 )
-         D = Z( 4*I0+PP-3 )
-         DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
-            Z( I4-2*PP-2 ) = D + Z( I4-1 )
-            IF( Z( I4-1 ).LE.TOL2*D ) THEN
-               Z( I4-1 ) = -ZERO
-               Z( I4-2*PP-2 ) = D
-               Z( I4-2*PP ) = ZERO
-               D = Z( I4+1 )
-            ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
-     $               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
-               TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
-               Z( I4-2*PP ) = Z( I4-1 )*TEMP
-               D = D*TEMP
-            ELSE
-               Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
-               D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
-            END IF
-            EMIN = MIN( EMIN, Z( I4-2*PP ) )
-   60    CONTINUE 
-         Z( 4*N0-PP-2 ) = D
-*
-*        Now find qmax.
-*
-         QMAX = Z( 4*I0-PP-2 )
-         DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
-            QMAX = MAX( QMAX, Z( I4 ) )
-   70    CONTINUE
-*
-*        Prepare for the next iteration on K.
-*
-         PP = 1 - PP
-   80 CONTINUE
-*
-*     Initialise variables to pass to DLAZQ3
-*
-      TTYPE = 0
-      DMIN1 = ZERO
-      DMIN2 = ZERO
-      DN    = ZERO
-      DN1   = ZERO
-      DN2   = ZERO
-      TAU   = ZERO
-*
-      ITER = 2
-      NFAIL = 0
-      NDIV = 2*( N0-I0 )
-*
-      DO 140 IWHILA = 1, N + 1
-         IF( N0.LT.1 ) 
-     $      GO TO 150
-*
-*        While array unfinished do 
-*
-*        E(N0) holds the value of SIGMA when submatrix in I0:N0
-*        splits from the rest of the array, but is negated.
-*      
-         DESIG = ZERO
-         IF( N0.EQ.N ) THEN
-            SIGMA = ZERO
-         ELSE
-            SIGMA = -Z( 4*N0-1 )
-         END IF
-         IF( SIGMA.LT.ZERO ) THEN
-            INFO = 1
-            RETURN
-         END IF
-*
-*        Find last unreduced submatrix's top index I0, find QMAX and
-*        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
-*
-         EMAX = ZERO 
-         IF( N0.GT.I0 ) THEN
-            EMIN = ABS( Z( 4*N0-5 ) )
-         ELSE
-            EMIN = ZERO
-         END IF
-         QMIN = Z( 4*N0-3 )
-         QMAX = QMIN
-         DO 90 I4 = 4*N0, 8, -4
-            IF( Z( I4-5 ).LE.ZERO )
-     $         GO TO 100
-            IF( QMIN.GE.FOUR*EMAX ) THEN
-               QMIN = MIN( QMIN, Z( I4-3 ) )
-               EMAX = MAX( EMAX, Z( I4-5 ) )
-            END IF
-            QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
-            EMIN = MIN( EMIN, Z( I4-5 ) )
-   90    CONTINUE
-         I4 = 4 
-*
-  100    CONTINUE
-         I0 = I4 / 4
-*
-*        Store EMIN for passing to DLAZQ3.
-*
-         Z( 4*N0-1 ) = EMIN
-*
-*        Put -(initial shift) into DMIN.
-*
-         DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
-*
-*        Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong.
-*
-         PP = 0 
-*
-         NBIG = 30*( N0-I0+1 )
-         DO 120 IWHILB = 1, NBIG
-            IF( I0.GT.N0 ) 
-     $         GO TO 130
-*
-*           While submatrix unfinished take a good dqds step.
-*
-            CALL DLAZQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
-     $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
-     $                   DN2, TAU )
-*
-            PP = 1 - PP
-*
-*           When EMIN is very small check for splits.
-*
-            IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
-               IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
-     $             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
-                  SPLT = I0 - 1
-                  QMAX = Z( 4*I0-3 )
-                  EMIN = Z( 4*I0-1 )
-                  OLDEMN = Z( 4*I0 )
-                  DO 110 I4 = 4*I0, 4*( N0-3 ), 4
-                     IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
-     $                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
-                        Z( I4-1 ) = -SIGMA
-                        SPLT = I4 / 4
-                        QMAX = ZERO
-                        EMIN = Z( I4+3 )
-                        OLDEMN = Z( I4+4 )
-                     ELSE
-                        QMAX = MAX( QMAX, Z( I4+1 ) )
-                        EMIN = MIN( EMIN, Z( I4-1 ) )
-                        OLDEMN = MIN( OLDEMN, Z( I4 ) )
-                     END IF
-  110             CONTINUE
-                  Z( 4*N0-1 ) = EMIN
-                  Z( 4*N0 ) = OLDEMN
-                  I0 = SPLT + 1
-               END IF
-            END IF
-*
-  120    CONTINUE
-*
-         INFO = 2
-         RETURN
-*
-*        end IWHILB
-*
-  130    CONTINUE
-*
-  140 CONTINUE
-*
-      INFO = 3
-      RETURN
-*
-*     end IWHILA   
-*
-  150 CONTINUE
-*      
-*     Move q's to the front.
-*      
-      DO 160 K = 2, N
-         Z( K ) = Z( 4*K-3 )
-  160 CONTINUE
-*      
-*     Sort and compute sum of eigenvalues.
-*
-      CALL DLASRT( 'D', N, Z, IINFO )
-*
-      E = ZERO
-      DO 170 K = N, 1, -1
-         E = E + Z( K )
-  170 CONTINUE
-*
-*     Store trace, sum(eigenvalues) and information on performance.
-*
-      Z( 2*N+1 ) = TRACE 
-      Z( 2*N+2 ) = E
-      Z( 2*N+3 ) = DBLE( ITER )
-      Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )
-      Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER )
-      RETURN
-*
-*     End of DLASQ2
-*
-      END