LAPACK: SRC/dlasq3.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
00002      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
00003      $                   DN2, G, TAU )
00004 *
00005 *  -- LAPACK routine (version 3.2.2)                                    --
00006 *
00007 *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
00008 *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
00009 *  -- Berkeley                                                        --
00010 *  -- June 2010                                                       --
00011 *
00012 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00013 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00014 *
00015 *     .. Scalar Arguments ..
00016       LOGICAL            IEEE
00017       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
00018       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
00019      $                   QMAX, SIGMA, TAU
00020 *     ..
00021 *     .. Array Arguments ..
00022       DOUBLE PRECISION   Z( * )
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *  DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
00029 *  In case of failure it changes shifts, and tries again until output
00030 *  is positive.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  I0     (input) INTEGER
00036 *         First index.
00037 *
00038 *  N0     (input/output) INTEGER
00039 *         Last index.
00040 *
00041 *  Z      (input) DOUBLE PRECISION array, dimension ( 4*N )
00042 *         Z holds the qd array.
00043 *
00044 *  PP     (input/output) INTEGER
00045 *         PP=0 for ping, PP=1 for pong.
00046 *         PP=2 indicates that flipping was applied to the Z array   
00047 *         and that the initial tests for deflation should not be 
00048 *         performed.
00049 *
00050 *  DMIN   (output) DOUBLE PRECISION
00051 *         Minimum value of d.
00052 *
00053 *  SIGMA  (output) DOUBLE PRECISION
00054 *         Sum of shifts used in current segment.
00055 *
00056 *  DESIG  (input/output) DOUBLE PRECISION
00057 *         Lower order part of SIGMA
00058 *
00059 *  QMAX   (input) DOUBLE PRECISION
00060 *         Maximum value of q.
00061 *
00062 *  NFAIL  (output) INTEGER
00063 *         Number of times shift was too big.
00064 *
00065 *  ITER   (output) INTEGER
00066 *         Number of iterations.
00067 *
00068 *  NDIV   (output) INTEGER
00069 *         Number of divisions.
00070 *
00071 *  IEEE   (input) LOGICAL
00072 *         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
00073 *
00074 *  TTYPE  (input/output) INTEGER
00075 *         Shift type.
00076 *
00077 *  DMIN1  (input/output) DOUBLE PRECISION
00078 *
00079 *  DMIN2  (input/output) DOUBLE PRECISION
00080 *
00081 *  DN     (input/output) DOUBLE PRECISION
00082 *
00083 *  DN1    (input/output) DOUBLE PRECISION
00084 *
00085 *  DN2    (input/output) DOUBLE PRECISION
00086 *
00087 *  G      (input/output) DOUBLE PRECISION
00088 *
00089 *  TAU    (input/output) DOUBLE PRECISION
00090 *
00091 *         These are passed as arguments in order to save their values
00092 *         between calls to DLASQ3.
00093 *
00094 *  =====================================================================
00095 *
00096 *     .. Parameters ..
00097       DOUBLE PRECISION   CBIAS
00098       PARAMETER          ( CBIAS = 1.50D0 )
00099       DOUBLE PRECISION   ZERO, QURTR, HALF, ONE, TWO, HUNDRD
00100       PARAMETER          ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
00101      $                     ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
00102 *     ..
00103 *     .. Local Scalars ..
00104       INTEGER            IPN4, J4, N0IN, NN, TTYPE
00105       DOUBLE PRECISION   EPS, S, T, TEMP, TOL, TOL2
00106 *     ..
00107 *     .. External Subroutines ..
00108       EXTERNAL           DLASQ4, DLASQ5, DLASQ6
00109 *     ..
00110 *     .. External Function ..
00111       DOUBLE PRECISION   DLAMCH
00112       LOGICAL            DISNAN
00113       EXTERNAL           DISNAN, DLAMCH
00114 *     ..
00115 *     .. Intrinsic Functions ..
00116       INTRINSIC          ABS, MAX, MIN, SQRT
00117 *     ..
00118 *     .. Executable Statements ..
00119 *
00120       N0IN = N0
00121       EPS = DLAMCH( 'Precision' )
00122       TOL = EPS*HUNDRD
00123       TOL2 = TOL**2
00124 *
00125 *     Check for deflation.
00126 *
00127    10 CONTINUE
00128 *
00129       IF( N0.LT.I0 )
00130      $   RETURN
00131       IF( N0.EQ.I0 )
00132      $   GO TO 20
00133       NN = 4*N0 + PP
00134       IF( N0.EQ.( I0+1 ) )
00135      $   GO TO 40
00136 *
00137 *     Check whether E(N0-1) is negligible, 1 eigenvalue.
00138 *
00139       IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
00140      $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
00141      $   GO TO 30
00142 *
00143    20 CONTINUE
00144 *
00145       Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
00146       N0 = N0 - 1
00147       GO TO 10
00148 *
00149 *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
00150 *
00151    30 CONTINUE
00152 *
00153       IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
00154      $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
00155      $   GO TO 50
00156 *
00157    40 CONTINUE
00158 *
00159       IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
00160          S = Z( NN-3 )
00161          Z( NN-3 ) = Z( NN-7 )
00162          Z( NN-7 ) = S
00163       END IF
00164       IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
00165          T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
00166          S = Z( NN-3 )*( Z( NN-5 ) / T )
00167          IF( S.LE.T ) THEN
00168             S = Z( NN-3 )*( Z( NN-5 ) /
00169      $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
00170          ELSE
00171             S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
00172          END IF
00173          T = Z( NN-7 ) + ( S+Z( NN-5 ) )
00174          Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
00175          Z( NN-7 ) = T
00176       END IF
00177       Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
00178       Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
00179       N0 = N0 - 2
00180       GO TO 10
00181 *
00182    50 CONTINUE
00183       IF( PP.EQ.2 ) 
00184      $   PP = 0
00185 *
00186 *     Reverse the qd-array, if warranted.
00187 *
00188       IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
00189          IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
00190             IPN4 = 4*( I0+N0 )
00191             DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
00192                TEMP = Z( J4-3 )
00193                Z( J4-3 ) = Z( IPN4-J4-3 )
00194                Z( IPN4-J4-3 ) = TEMP
00195                TEMP = Z( J4-2 )
00196                Z( J4-2 ) = Z( IPN4-J4-2 )
00197                Z( IPN4-J4-2 ) = TEMP
00198                TEMP = Z( J4-1 )
00199                Z( J4-1 ) = Z( IPN4-J4-5 )
00200                Z( IPN4-J4-5 ) = TEMP
00201                TEMP = Z( J4 )
00202                Z( J4 ) = Z( IPN4-J4-4 )
00203                Z( IPN4-J4-4 ) = TEMP
00204    60       CONTINUE
00205             IF( N0-I0.LE.4 ) THEN
00206                Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
00207                Z( 4*N0-PP ) = Z( 4*I0-PP )
00208             END IF
00209             DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
00210             Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
00211      $                            Z( 4*I0+PP+3 ) )
00212             Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
00213      $                          Z( 4*I0-PP+4 ) )
00214             QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
00215             DMIN = -ZERO
00216          END IF
00217       END IF
00218 *
00219 *     Choose a shift.
00220 *
00221       CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
00222      $             DN2, TAU, TTYPE, G )
00223 *
00224 *     Call dqds until DMIN > 0.
00225 *
00226    70 CONTINUE
00227 *
00228       CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
00229      $             DN1, DN2, IEEE )
00230 *
00231       NDIV = NDIV + ( N0-I0+2 )
00232       ITER = ITER + 1
00233 *
00234 *     Check status.
00235 *
00236       IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
00237 *
00238 *        Success.
00239 *
00240          GO TO 90
00241 *
00242       ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 
00243      $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
00244      $         ABS( DN ).LT.TOL*SIGMA ) THEN
00245 *
00246 *        Convergence hidden by negative DN.
00247 *
00248          Z( 4*( N0-1 )-PP+2 ) = ZERO
00249          DMIN = ZERO
00250          GO TO 90
00251       ELSE IF( DMIN.LT.ZERO ) THEN
00252 *
00253 *        TAU too big. Select new TAU and try again.
00254 *
00255          NFAIL = NFAIL + 1
00256          IF( TTYPE.LT.-22 ) THEN
00257 *
00258 *           Failed twice. Play it safe.
00259 *
00260             TAU = ZERO
00261          ELSE IF( DMIN1.GT.ZERO ) THEN
00262 *
00263 *           Late failure. Gives excellent shift.
00264 *
00265             TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
00266             TTYPE = TTYPE - 11
00267          ELSE
00268 *
00269 *           Early failure. Divide by 4.
00270 *
00271             TAU = QURTR*TAU
00272             TTYPE = TTYPE - 12
00273          END IF
00274          GO TO 70
00275       ELSE IF( DISNAN( DMIN ) ) THEN
00276 *
00277 *        NaN.
00278 *
00279          IF( TAU.EQ.ZERO ) THEN
00280             GO TO 80
00281          ELSE
00282             TAU = ZERO
00283             GO TO 70
00284          END IF
00285       ELSE
00286 *            
00287 *        Possible underflow. Play it safe.
00288 *
00289          GO TO 80
00290       END IF
00291 *
00292 *     Risk of underflow.
00293 *
00294    80 CONTINUE
00295       CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
00296       NDIV = NDIV + ( N0-I0+2 )
00297       ITER = ITER + 1
00298       TAU = ZERO
00299 *
00300    90 CONTINUE
00301       IF( TAU.LT.SIGMA ) THEN
00302          DESIG = DESIG + TAU
00303          T = SIGMA + DESIG
00304          DESIG = DESIG - ( T-SIGMA )
00305       ELSE
00306          T = SIGMA + TAU
00307          DESIG = SIGMA - ( T-TAU ) + DESIG
00308       END IF
00309       SIGMA = T
00310 *
00311       RETURN
00312 *
00313 *     End of DLASQ3
00314 *
00315       END