Original Version

1 files changed, 295 insertions, 0 deletions

diff --git a/src/fortran/lapack/zlaic1.f b/src/fortran/lapack/zlaic1.f
new file mode 100644
index 0000000..589f088
--- /dev/null
+++ b/src/fortran/lapack/zlaic1.f
@@ -0,0 +1,295 @@
+      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