Mercurial > octave-nkf

diff libcruft/lapack/slazq3.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/lapack/slazq3.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,302 @@
+      SUBROUTINE SLAZQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
+     $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
+     $                   DN2, TAU )
+*
+*  -- LAPACK auxiliary routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE
+      INTEGER            I0, ITER, N0, NDIV, NFAIL, PP, TTYPE
+      REAL               DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, QMAX,
+     $                   SIGMA, TAU
+*     ..
+*     .. Array Arguments ..
+      REAL               Z( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
+*  In case of failure it changes shifts, and tries again until output
+*  is positive.
+*
+*  Arguments
+*  =========
+*
+*  I0     (input) INTEGER
+*         First index.
+*
+*  N0     (input) INTEGER
+*         Last index.
+*
+*  Z      (input) REAL array, dimension ( 4*N )
+*         Z holds the qd array.
+*
+*  PP     (input) INTEGER
+*         PP=0 for ping, PP=1 for pong.
+*
+*  DMIN   (output) REAL
+*         Minimum value of d.
+*
+*  SIGMA  (output) REAL
+*         Sum of shifts used in current segment.
+*
+*  DESIG  (input/output) REAL
+*         Lower order part of SIGMA
+*
+*  QMAX   (input) REAL
+*         Maximum value of q.
+*
+*  NFAIL  (output) INTEGER
+*         Number of times shift was too big.
+*
+*  ITER   (output) INTEGER
+*         Number of iterations.
+*
+*  NDIV   (output) INTEGER
+*         Number of divisions.
+*
+*  IEEE   (input) LOGICAL
+*         Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).
+*
+*  TTYPE  (input/output) INTEGER
+*         Shift type.  TTYPE is passed as an argument in order to save
+*         its value between calls to SLAZQ3
+*
+*  DMIN1  (input/output) REAL
+*  DMIN2  (input/output) REAL
+*  DN     (input/output) REAL
+*  DN1    (input/output) REAL
+*  DN2    (input/output) REAL
+*  TAU    (input/output) REAL
+*         These are passed as arguments in order to save their values
+*         between calls to SLAZQ3
+*
+*  This is a thread safe version of SLASQ3, which passes TTYPE, DMIN1,
+*  DMIN2, DN, DN1. DN2 and TAU through the argument list in place of
+*  declaring them in a SAVE statment.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               CBIAS
+      PARAMETER          ( CBIAS = 1.50E0 )
+      REAL               ZERO, QURTR, HALF, ONE, TWO, HUNDRD
+      PARAMETER          ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0,
+     $                     ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            IPN4, J4, N0IN, NN
+      REAL               EPS, G, S, SAFMIN, T, TEMP, TOL, TOL2
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLASQ5, SLASQ6, SLAZQ4
+*     ..
+*     .. External Function ..
+      REAL               SLAMCH
+      EXTERNAL           SLAMCH
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MIN, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+      N0IN   = N0
+      EPS    = SLAMCH( 'Precision' )
+      SAFMIN = SLAMCH( 'Safe minimum' )
+      TOL    = EPS*HUNDRD
+      TOL2   = TOL**2
+      G      = ZERO
+*
+*     Check for deflation.
+*
+   10 CONTINUE
+*
+      IF( N0.LT.I0 )
+     $   RETURN
+      IF( N0.EQ.I0 )
+     $   GO TO 20
+      NN = 4*N0 + PP
+      IF( N0.EQ.( I0+1 ) )
+     $   GO TO 40
+*
+*     Check whether E(N0-1) is negligible, 1 eigenvalue.
+*
+      IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
+     $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
+     $   GO TO 30
+*
+   20 CONTINUE
+*
+      Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
+      N0 = N0 - 1
+      GO TO 10
+*
+*     Check  whether E(N0-2) is negligible, 2 eigenvalues.
+*
+   30 CONTINUE
+*
+      IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
+     $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
+     $   GO TO 50
+*
+   40 CONTINUE
+*
+      IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
+         S = Z( NN-3 )
+         Z( NN-3 ) = Z( NN-7 )
+         Z( NN-7 ) = S
+      END IF
+      IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
+         T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
+         S = Z( NN-3 )*( Z( NN-5 ) / T )
+         IF( S.LE.T ) THEN
+            S = Z( NN-3 )*( Z( NN-5 ) /
+     $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
+         ELSE
+            S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
+         END IF
+         T = Z( NN-7 ) + ( S+Z( NN-5 ) )
+         Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
+         Z( NN-7 ) = T
+      END IF
+      Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
+      Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
+      N0 = N0 - 2
+      GO TO 10
+*
+   50 CONTINUE
+*
+*     Reverse the qd-array, if warranted.
+*
+      IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
+         IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
+            IPN4 = 4*( I0+N0 )
+            DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
+               TEMP = Z( J4-3 )
+               Z( J4-3 ) = Z( IPN4-J4-3 )
+               Z( IPN4-J4-3 ) = TEMP
+               TEMP = Z( J4-2 )
+               Z( J4-2 ) = Z( IPN4-J4-2 )
+               Z( IPN4-J4-2 ) = TEMP
+               TEMP = Z( J4-1 )
+               Z( J4-1 ) = Z( IPN4-J4-5 )
+               Z( IPN4-J4-5 ) = TEMP
+               TEMP = Z( J4 )
+               Z( J4 ) = Z( IPN4-J4-4 )
+               Z( IPN4-J4-4 ) = TEMP
+   60       CONTINUE
+            IF( N0-I0.LE.4 ) THEN
+               Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
+               Z( 4*N0-PP ) = Z( 4*I0-PP )
+            END IF
+            DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
+            Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
+     $                            Z( 4*I0+PP+3 ) )
+            Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
+     $                          Z( 4*I0-PP+4 ) )
+            QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
+            DMIN = -ZERO
+         END IF
+      END IF
+*
+      IF( DMIN.LT.ZERO .OR. SAFMIN*QMAX.LT.MIN( Z( 4*N0+PP-1 ),
+     $    Z( 4*N0+PP-9 ), DMIN2+Z( 4*N0-PP ) ) ) THEN
+*
+*        Choose a shift.
+*
+         CALL SLAZQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
+     $                DN2, TAU, TTYPE, G )
+*
+*        Call dqds until DMIN > 0.
+*
+   80    CONTINUE
+*
+         CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
+     $                DN1, DN2, IEEE )
+*
+         NDIV = NDIV + ( N0-I0+2 )
+         ITER = ITER + 1
+*
+*        Check status.
+*
+         IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
+*
+*           Success.
+*
+            GO TO 100
+*
+         ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
+     $            Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
+     $            ABS( DN ).LT.TOL*SIGMA ) THEN
+*
+*           Convergence hidden by negative DN.
+*
+            Z( 4*( N0-1 )-PP+2 ) = ZERO
+            DMIN = ZERO
+            GO TO 100
+         ELSE IF( DMIN.LT.ZERO ) THEN
+*
+*           TAU too big. Select new TAU and try again.
+*
+            NFAIL = NFAIL + 1
+            IF( TTYPE.LT.-22 ) THEN
+*
+*              Failed twice. Play it safe.
+*
+               TAU = ZERO
+            ELSE IF( DMIN1.GT.ZERO ) THEN
+*
+*              Late failure. Gives excellent shift.
+*
+               TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
+               TTYPE = TTYPE - 11
+            ELSE
+*
+*              Early failure. Divide by 4.
+*
+               TAU = QURTR*TAU
+               TTYPE = TTYPE - 12
+            END IF
+            GO TO 80
+         ELSE IF( DMIN.NE.DMIN ) THEN
+*
+*           NaN.
+*
+            TAU = ZERO
+            GO TO 80
+         ELSE
+*
+*           Possible underflow. Play it safe.
+*
+            GO TO 90
+         END IF
+      END IF
+*
+*     Risk of underflow.
+*
+   90 CONTINUE
+      CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
+      NDIV = NDIV + ( N0-I0+2 )
+      ITER = ITER + 1
+      TAU = ZERO
+*
+  100 CONTINUE
+      IF( TAU.LT.SIGMA ) THEN
+         DESIG = DESIG + TAU
+         T = SIGMA + DESIG
+         DESIG = DESIG - ( T-SIGMA )
+      ELSE
+         T = SIGMA + TAU
+         DESIG = SIGMA - ( T-TAU ) + DESIG
+      END IF
+      SIGMA = T
+*
+      RETURN
+*
+*     End of SLAZQ3
+*
+      END