LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasq3.f
Go to the documentation of this file.
00001 *> \brief \b DLASQ3
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASQ3 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
00022 *                          ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
00023 *                          DN2, G, TAU )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       LOGICAL            IEEE
00027 *       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
00028 *       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
00029 *      $                   QMAX, SIGMA, TAU
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       DOUBLE PRECISION   Z( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
00042 *> In case of failure it changes shifts, and tries again until output
00043 *> is positive.
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] I0
00050 *> \verbatim
00051 *>          I0 is INTEGER
00052 *>         First index.
00053 *> \endverbatim
00054 *>
00055 *> \param[in,out] N0
00056 *> \verbatim
00057 *>          N0 is INTEGER
00058 *>         Last index.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] Z
00062 *> \verbatim
00063 *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
00064 *>         Z holds the qd array.
00065 *> \endverbatim
00066 *>
00067 *> \param[in,out] PP
00068 *> \verbatim
00069 *>          PP is INTEGER
00070 *>         PP=0 for ping, PP=1 for pong.
00071 *>         PP=2 indicates that flipping was applied to the Z array   
00072 *>         and that the initial tests for deflation should not be 
00073 *>         performed.
00074 *> \endverbatim
00075 *>
00076 *> \param[out] DMIN
00077 *> \verbatim
00078 *>          DMIN is DOUBLE PRECISION
00079 *>         Minimum value of d.
00080 *> \endverbatim
00081 *>
00082 *> \param[out] SIGMA
00083 *> \verbatim
00084 *>          SIGMA is DOUBLE PRECISION
00085 *>         Sum of shifts used in current segment.
00086 *> \endverbatim
00087 *>
00088 *> \param[in,out] DESIG
00089 *> \verbatim
00090 *>          DESIG is DOUBLE PRECISION
00091 *>         Lower order part of SIGMA
00092 *> \endverbatim
00093 *>
00094 *> \param[in] QMAX
00095 *> \verbatim
00096 *>          QMAX is DOUBLE PRECISION
00097 *>         Maximum value of q.
00098 *> \endverbatim
00099 *>
00100 *> \param[out] NFAIL
00101 *> \verbatim
00102 *>          NFAIL is INTEGER
00103 *>         Number of times shift was too big.
00104 *> \endverbatim
00105 *>
00106 *> \param[out] ITER
00107 *> \verbatim
00108 *>          ITER is INTEGER
00109 *>         Number of iterations.
00110 *> \endverbatim
00111 *>
00112 *> \param[out] NDIV
00113 *> \verbatim
00114 *>          NDIV is INTEGER
00115 *>         Number of divisions.
00116 *> \endverbatim
00117 *>
00118 *> \param[in] IEEE
00119 *> \verbatim
00120 *>          IEEE is LOGICAL
00121 *>         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
00122 *> \endverbatim
00123 *>
00124 *> \param[in,out] TTYPE
00125 *> \verbatim
00126 *>          TTYPE is INTEGER
00127 *>         Shift type.
00128 *> \endverbatim
00129 *>
00130 *> \param[in,out] DMIN1
00131 *> \verbatim
00132 *>          DMIN1 is DOUBLE PRECISION
00133 *> \endverbatim
00134 *>
00135 *> \param[in,out] DMIN2
00136 *> \verbatim
00137 *>          DMIN2 is DOUBLE PRECISION
00138 *> \endverbatim
00139 *>
00140 *> \param[in,out] DN
00141 *> \verbatim
00142 *>          DN is DOUBLE PRECISION
00143 *> \endverbatim
00144 *>
00145 *> \param[in,out] DN1
00146 *> \verbatim
00147 *>          DN1 is DOUBLE PRECISION
00148 *> \endverbatim
00149 *>
00150 *> \param[in,out] DN2
00151 *> \verbatim
00152 *>          DN2 is DOUBLE PRECISION
00153 *> \endverbatim
00154 *>
00155 *> \param[in,out] G
00156 *> \verbatim
00157 *>          G is DOUBLE PRECISION
00158 *> \endverbatim
00159 *>
00160 *> \param[in,out] TAU
00161 *> \verbatim
00162 *>          TAU is DOUBLE PRECISION
00163 *>
00164 *>         These are passed as arguments in order to save their values
00165 *>         between calls to DLASQ3.
00166 *> \endverbatim
00167 *
00168 *  Authors:
00169 *  ========
00170 *
00171 *> \author Univ. of Tennessee 
00172 *> \author Univ. of California Berkeley 
00173 *> \author Univ. of Colorado Denver 
00174 *> \author NAG Ltd. 
00175 *
00176 *> \date April 2012
00177 *
00178 *> \ingroup auxOTHERcomputational
00179 *
00180 *  =====================================================================
00181       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
00182      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
00183      $                   DN2, G, TAU )
00184 *
00185 *  -- LAPACK computational routine (version 3.4.1) --
00186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00188 *     April 2012
00189 *
00190 *     .. Scalar Arguments ..
00191       LOGICAL            IEEE
00192       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
00193       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
00194      $                   QMAX, SIGMA, TAU
00195 *     ..
00196 *     .. Array Arguments ..
00197       DOUBLE PRECISION   Z( * )
00198 *     ..
00199 *
00200 *  =====================================================================
00201 *
00202 *     .. Parameters ..
00203       DOUBLE PRECISION   CBIAS
00204       PARAMETER          ( CBIAS = 1.50D0 )
00205       DOUBLE PRECISION   ZERO, QURTR, HALF, ONE, TWO, HUNDRD
00206       PARAMETER          ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
00207      $                     ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
00208 *     ..
00209 *     .. Local Scalars ..
00210       INTEGER            IPN4, J4, N0IN, NN, TTYPE
00211       DOUBLE PRECISION   EPS, S, T, TEMP, TOL, TOL2
00212 *     ..
00213 *     .. External Subroutines ..
00214       EXTERNAL           DLASQ4, DLASQ5, DLASQ6
00215 *     ..
00216 *     .. External Function ..
00217       DOUBLE PRECISION   DLAMCH
00218       LOGICAL            DISNAN
00219       EXTERNAL           DISNAN, DLAMCH
00220 *     ..
00221 *     .. Intrinsic Functions ..
00222       INTRINSIC          ABS, MAX, MIN, SQRT
00223 *     ..
00224 *     .. Executable Statements ..
00225 *
00226       N0IN = N0
00227       EPS = DLAMCH( 'Precision' )
00228       TOL = EPS*HUNDRD
00229       TOL2 = TOL**2
00230 *
00231 *     Check for deflation.
00232 *
00233    10 CONTINUE
00234 *
00235       IF( N0.LT.I0 )
00236      $   RETURN
00237       IF( N0.EQ.I0 )
00238      $   GO TO 20
00239       NN = 4*N0 + PP
00240       IF( N0.EQ.( I0+1 ) )
00241      $   GO TO 40
00242 *
00243 *     Check whether E(N0-1) is negligible, 1 eigenvalue.
00244 *
00245       IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
00246      $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
00247      $   GO TO 30
00248 *
00249    20 CONTINUE
00250 *
00251       Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
00252       N0 = N0 - 1
00253       GO TO 10
00254 *
00255 *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
00256 *
00257    30 CONTINUE
00258 *
00259       IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
00260      $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
00261      $   GO TO 50
00262 *
00263    40 CONTINUE
00264 *
00265       IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
00266          S = Z( NN-3 )
00267          Z( NN-3 ) = Z( NN-7 )
00268          Z( NN-7 ) = S
00269       END IF
00270       IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
00271          T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
00272          S = Z( NN-3 )*( Z( NN-5 ) / T )
00273          IF( S.LE.T ) THEN
00274             S = Z( NN-3 )*( Z( NN-5 ) /
00275      $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
00276          ELSE
00277             S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
00278          END IF
00279          T = Z( NN-7 ) + ( S+Z( NN-5 ) )
00280          Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
00281          Z( NN-7 ) = T
00282       END IF
00283       Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
00284       Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
00285       N0 = N0 - 2
00286       GO TO 10
00287 *
00288    50 CONTINUE
00289       IF( PP.EQ.2 ) 
00290      $   PP = 0
00291 *
00292 *     Reverse the qd-array, if warranted.
00293 *
00294       IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
00295          IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
00296             IPN4 = 4*( I0+N0 )
00297             DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
00298                TEMP = Z( J4-3 )
00299                Z( J4-3 ) = Z( IPN4-J4-3 )
00300                Z( IPN4-J4-3 ) = TEMP
00301                TEMP = Z( J4-2 )
00302                Z( J4-2 ) = Z( IPN4-J4-2 )
00303                Z( IPN4-J4-2 ) = TEMP
00304                TEMP = Z( J4-1 )
00305                Z( J4-1 ) = Z( IPN4-J4-5 )
00306                Z( IPN4-J4-5 ) = TEMP
00307                TEMP = Z( J4 )
00308                Z( J4 ) = Z( IPN4-J4-4 )
00309                Z( IPN4-J4-4 ) = TEMP
00310    60       CONTINUE
00311             IF( N0-I0.LE.4 ) THEN
00312                Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
00313                Z( 4*N0-PP ) = Z( 4*I0-PP )
00314             END IF
00315             DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
00316             Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
00317      $                            Z( 4*I0+PP+3 ) )
00318             Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
00319      $                          Z( 4*I0-PP+4 ) )
00320             QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
00321             DMIN = -ZERO
00322          END IF
00323       END IF
00324 *
00325 *     Choose a shift.
00326 *
00327       CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
00328      $             DN2, TAU, TTYPE, G )
00329 *
00330 *     Call dqds until DMIN > 0.
00331 *
00332    70 CONTINUE
00333 *
00334       CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
00335      $             DN1, DN2, IEEE, EPS )
00336 *
00337       NDIV = NDIV + ( N0-I0+2 )
00338       ITER = ITER + 1
00339 *
00340 *     Check status.
00341 *
00342       IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
00343 *
00344 *        Success.
00345 *
00346          GO TO 90
00347 *
00348       ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 
00349      $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
00350      $         ABS( DN ).LT.TOL*SIGMA ) THEN
00351 *
00352 *        Convergence hidden by negative DN.
00353 *
00354          Z( 4*( N0-1 )-PP+2 ) = ZERO
00355          DMIN = ZERO
00356          GO TO 90
00357       ELSE IF( DMIN.LT.ZERO ) THEN
00358 *
00359 *        TAU too big. Select new TAU and try again.
00360 *
00361          NFAIL = NFAIL + 1
00362          IF( TTYPE.LT.-22 ) THEN
00363 *
00364 *           Failed twice. Play it safe.
00365 *
00366             TAU = ZERO
00367          ELSE IF( DMIN1.GT.ZERO ) THEN
00368 *
00369 *           Late failure. Gives excellent shift.
00370 *
00371             TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
00372             TTYPE = TTYPE - 11
00373          ELSE
00374 *
00375 *           Early failure. Divide by 4.
00376 *
00377             TAU = QURTR*TAU
00378             TTYPE = TTYPE - 12
00379          END IF
00380          GO TO 70
00381       ELSE IF( DISNAN( DMIN ) ) THEN
00382 *
00383 *        NaN.
00384 *
00385          IF( TAU.EQ.ZERO ) THEN
00386             GO TO 80
00387          ELSE
00388             TAU = ZERO
00389             GO TO 70
00390          END IF
00391       ELSE
00392 *            
00393 *        Possible underflow. Play it safe.
00394 *
00395          GO TO 80
00396       END IF
00397 *
00398 *     Risk of underflow.
00399 *
00400    80 CONTINUE
00401       CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
00402       NDIV = NDIV + ( N0-I0+2 )
00403       ITER = ITER + 1
00404       TAU = ZERO
00405 *
00406    90 CONTINUE
00407       IF( TAU.LT.SIGMA ) THEN
00408          DESIG = DESIG + TAU
00409          T = SIGMA + DESIG
00410          DESIG = DESIG - ( T-SIGMA )
00411       ELSE
00412          T = SIGMA + TAU
00413          DESIG = SIGMA - ( T-TAU ) + DESIG
00414       END IF
00415       SIGMA = T
00416 *
00417       RETURN
00418 *
00419 *     End of DLASQ3
00420 *
00421       END
 All Files Functions