Moving lapack to right place

1 files changed, 0 insertions, 295 deletions

diff --git a/src/lib/lapack/zlaic1.f b/src/lib/lapack/zlaic1.f
deleted file mode 100644
index 589f0889..00000000
--- a/src/lib/lapack/zlaic1.f
+++ /dev/null
@@ -1,295 +0,0 @@
-      SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
-*
-*  -- LAPACK auxiliary routine (version 3.1) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      INTEGER            J, JOB
-      DOUBLE PRECISION   SEST, SESTPR
-      COMPLEX*16         C, GAMMA, S
-*     ..
-*     .. Array Arguments ..
-      COMPLEX*16         W( J ), X( J )
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  ZLAIC1 applies one step of incremental condition estimation in
-*  its simplest version:
-*
-*  Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
-*  lower triangular matrix L, such that
-*           twonorm(L*x) = sest
-*  Then ZLAIC1 computes sestpr, s, c such that
-*  the vector
-*                  [ s*x ]
-*           xhat = [  c  ]
-*  is an approximate singular vector of
-*                  [ L     0  ]
-*           Lhat = [ w' gamma ]
-*  in the sense that
-*           twonorm(Lhat*xhat) = sestpr.
-*
-*  Depending on JOB, an estimate for the largest or smallest singular
-*  value is computed.
-*
-*  Note that [s c]' and sestpr**2 is an eigenpair of the system
-*
-*      diag(sest*sest, 0) + [alpha  gamma] * [ conjg(alpha) ]
-*                                            [ conjg(gamma) ]
-*
-*  where  alpha =  conjg(x)'*w.
-*
-*  Arguments
-*  =========
-*
-*  JOB     (input) INTEGER
-*          = 1: an estimate for the largest singular value is computed.
-*          = 2: an estimate for the smallest singular value is computed.
-*
-*  J       (input) INTEGER
-*          Length of X and W
-*
-*  X       (input) COMPLEX*16 array, dimension (J)
-*          The j-vector x.
-*
-*  SEST    (input) DOUBLE PRECISION
-*          Estimated singular value of j by j matrix L
-*
-*  W       (input) COMPLEX*16 array, dimension (J)
-*          The j-vector w.
-*
-*  GAMMA   (input) COMPLEX*16
-*          The diagonal element gamma.
-*
-*  SESTPR  (output) DOUBLE PRECISION
-*          Estimated singular value of (j+1) by (j+1) matrix Lhat.
-*
-*  S       (output) COMPLEX*16
-*          Sine needed in forming xhat.
-*
-*  C       (output) COMPLEX*16
-*          Cosine needed in forming xhat.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   ZERO, ONE, TWO
-      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
-      DOUBLE PRECISION   HALF, FOUR
-      PARAMETER          ( HALF = 0.5D0, FOUR = 4.0D0 )
-*     ..
-*     .. Local Scalars ..
-      DOUBLE PRECISION   ABSALP, ABSEST, ABSGAM, B, EPS, NORMA, S1, S2,
-     $                   SCL, T, TEST, TMP, ZETA1, ZETA2
-      COMPLEX*16         ALPHA, COSINE, SINE
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          ABS, DCONJG, MAX, SQRT
-*     ..
-*     .. External Functions ..
-      DOUBLE PRECISION   DLAMCH
-      COMPLEX*16         ZDOTC
-      EXTERNAL           DLAMCH, ZDOTC
-*     ..
-*     .. Executable Statements ..
-*
-      EPS = DLAMCH( 'Epsilon' )
-      ALPHA = ZDOTC( J, X, 1, W, 1 )
-*
-      ABSALP = ABS( ALPHA )
-      ABSGAM = ABS( GAMMA )
-      ABSEST = ABS( SEST )
-*
-      IF( JOB.EQ.1 ) THEN
-*
-*        Estimating largest singular value
-*
-*        special cases
-*
-         IF( SEST.EQ.ZERO ) THEN
-            S1 = MAX( ABSGAM, ABSALP )
-            IF( S1.EQ.ZERO ) THEN
-               S = ZERO
-               C = ONE
-               SESTPR = ZERO
-            ELSE
-               S = ALPHA / S1
-               C = GAMMA / S1
-               TMP = SQRT( S*DCONJG( S )+C*DCONJG( C ) )
-               S = S / TMP
-               C = C / TMP
-               SESTPR = S1*TMP
-            END IF
-            RETURN
-         ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
-            S = ONE
-            C = ZERO
-            TMP = MAX( ABSEST, ABSALP )
-            S1 = ABSEST / TMP
-            S2 = ABSALP / TMP
-            SESTPR = TMP*SQRT( S1*S1+S2*S2 )
-            RETURN
-         ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
-            S1 = ABSGAM
-            S2 = ABSEST
-            IF( S1.LE.S2 ) THEN
-               S = ONE
-               C = ZERO
-               SESTPR = S2
-            ELSE
-               S = ZERO
-               C = ONE
-               SESTPR = S1
-            END IF
-            RETURN
-         ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
-            S1 = ABSGAM
-            S2 = ABSALP
-            IF( S1.LE.S2 ) THEN
-               TMP = S1 / S2
-               SCL = SQRT( ONE+TMP*TMP )
-               SESTPR = S2*SCL
-               S = ( ALPHA / S2 ) / SCL
-               C = ( GAMMA / S2 ) / SCL
-            ELSE
-               TMP = S2 / S1
-               SCL = SQRT( ONE+TMP*TMP )
-               SESTPR = S1*SCL
-               S = ( ALPHA / S1 ) / SCL
-               C = ( GAMMA / S1 ) / SCL
-            END IF
-            RETURN
-         ELSE
-*
-*           normal case
-*
-            ZETA1 = ABSALP / ABSEST
-            ZETA2 = ABSGAM / ABSEST
-*
-            B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
-            C = ZETA1*ZETA1
-            IF( B.GT.ZERO ) THEN
-               T = C / ( B+SQRT( B*B+C ) )
-            ELSE
-               T = SQRT( B*B+C ) - B
-            END IF
-*
-            SINE = -( ALPHA / ABSEST ) / T
-            COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
-            TMP = SQRT( SINE*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
-            S = SINE / TMP
-            C = COSINE / TMP
-            SESTPR = SQRT( T+ONE )*ABSEST
-            RETURN
-         END IF
-*
-      ELSE IF( JOB.EQ.2 ) THEN
-*
-*        Estimating smallest singular value
-*
-*        special cases
-*
-         IF( SEST.EQ.ZERO ) THEN
-            SESTPR = ZERO
-            IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
-               SINE = ONE
-               COSINE = ZERO
-            ELSE
-               SINE = -DCONJG( GAMMA )
-               COSINE = DCONJG( ALPHA )
-            END IF
-            S1 = MAX( ABS( SINE ), ABS( COSINE ) )
-            S = SINE / S1
-            C = COSINE / S1
-            TMP = SQRT( S*DCONJG( S )+C*DCONJG( C ) )
-            S = S / TMP
-            C = C / TMP
-            RETURN
-         ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
-            S = ZERO
-            C = ONE
-            SESTPR = ABSGAM
-            RETURN
-         ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
-            S1 = ABSGAM
-            S2 = ABSEST
-            IF( S1.LE.S2 ) THEN
-               S = ZERO
-               C = ONE
-               SESTPR = S1
-            ELSE
-               S = ONE
-               C = ZERO
-               SESTPR = S2
-            END IF
-            RETURN
-         ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
-            S1 = ABSGAM
-            S2 = ABSALP
-            IF( S1.LE.S2 ) THEN
-               TMP = S1 / S2
-               SCL = SQRT( ONE+TMP*TMP )
-               SESTPR = ABSEST*( TMP / SCL )
-               S = -( DCONJG( GAMMA ) / S2 ) / SCL
-               C = ( DCONJG( ALPHA ) / S2 ) / SCL
-            ELSE
-               TMP = S2 / S1
-               SCL = SQRT( ONE+TMP*TMP )
-               SESTPR = ABSEST / SCL
-               S = -( DCONJG( GAMMA ) / S1 ) / SCL
-               C = ( DCONJG( ALPHA ) / S1 ) / SCL
-            END IF
-            RETURN
-         ELSE
-*
-*           normal case
-*
-            ZETA1 = ABSALP / ABSEST
-            ZETA2 = ABSGAM / ABSEST
-*
-            NORMA = MAX( ONE+ZETA1*ZETA1+ZETA1*ZETA2,
-     $              ZETA1*ZETA2+ZETA2*ZETA2 )
-*
-*           See if root is closer to zero or to ONE
-*
-            TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
-            IF( TEST.GE.ZERO ) THEN
-*
-*              root is close to zero, compute directly
-*
-               B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
-               C = ZETA2*ZETA2
-               T = C / ( B+SQRT( ABS( B*B-C ) ) )
-               SINE = ( ALPHA / ABSEST ) / ( ONE-T )
-               COSINE = -( GAMMA / ABSEST ) / T
-               SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
-            ELSE
-*
-*              root is closer to ONE, shift by that amount
-*
-               B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
-               C = ZETA1*ZETA1
-               IF( B.GE.ZERO ) THEN
-                  T = -C / ( B+SQRT( B*B+C ) )
-               ELSE
-                  T = B - SQRT( B*B+C )
-               END IF
-               SINE = -( ALPHA / ABSEST ) / T
-               COSINE = -( GAMMA / ABSEST ) / ( ONE+T )
-               SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
-            END IF
-            TMP = SQRT( SINE*DCONJG( SINE )+COSINE*DCONJG( COSINE ) )
-            S = SINE / TMP
-            C = COSINE / TMP
-            RETURN
-*
-         END IF
-      END IF
-      RETURN
-*
-*     End of ZLAIC1
-*
-      END