Return to dlasq3.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence Lapack.

    1: *> \brief \b DLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLASQ3 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
   22: *                          ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
   23: *                          DN2, G, TAU )
   24: * 
   25: *       .. Scalar Arguments ..
   26: *       LOGICAL            IEEE
   27: *       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
   28: *       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
   29: *      $                   QMAX, SIGMA, TAU
   30: *       ..
   31: *       .. Array Arguments ..
   32: *       DOUBLE PRECISION   Z( * )
   33: *       ..
   34: *  
   35: *
   36: *> \par Purpose:
   37: *  =============
   38: *>
   39: *> \verbatim
   40: *>
   41: *> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
   42: *> In case of failure it changes shifts, and tries again until output
   43: *> is positive.
   44: *> \endverbatim
   45: *
   46: *  Arguments:
   47: *  ==========
   48: *
   49: *> \param[in] I0
   50: *> \verbatim
   51: *>          I0 is INTEGER
   52: *>         First index.
   53: *> \endverbatim
   54: *>
   55: *> \param[in,out] N0
   56: *> \verbatim
   57: *>          N0 is INTEGER
   58: *>         Last index.
   59: *> \endverbatim
   60: *>
   61: *> \param[in,out] Z
   62: *> \verbatim
   63: *>          Z is DOUBLE PRECISION array, dimension ( 4*N0 )
   64: *>         Z holds the qd array.
   65: *> \endverbatim
   66: *>
   67: *> \param[in,out] PP
   68: *> \verbatim
   69: *>          PP is INTEGER
   70: *>         PP=0 for ping, PP=1 for pong.
   71: *>         PP=2 indicates that flipping was applied to the Z array   
   72: *>         and that the initial tests for deflation should not be 
   73: *>         performed.
   74: *> \endverbatim
   75: *>
   76: *> \param[out] DMIN
   77: *> \verbatim
   78: *>          DMIN is DOUBLE PRECISION
   79: *>         Minimum value of d.
   80: *> \endverbatim
   81: *>
   82: *> \param[out] SIGMA
   83: *> \verbatim
   84: *>          SIGMA is DOUBLE PRECISION
   85: *>         Sum of shifts used in current segment.
   86: *> \endverbatim
   87: *>
   88: *> \param[in,out] DESIG
   89: *> \verbatim
   90: *>          DESIG is DOUBLE PRECISION
   91: *>         Lower order part of SIGMA
   92: *> \endverbatim
   93: *>
   94: *> \param[in] QMAX
   95: *> \verbatim
   96: *>          QMAX is DOUBLE PRECISION
   97: *>         Maximum value of q.
   98: *> \endverbatim
   99: *>
  100: *> \param[in,out] NFAIL
  101: *> \verbatim
  102: *>          NFAIL is INTEGER
  103: *>         Increment NFAIL by 1 each time the shift was too big.
  104: *> \endverbatim
  105: *>
  106: *> \param[in,out] ITER
  107: *> \verbatim
  108: *>          ITER is INTEGER
  109: *>         Increment ITER by 1 for each iteration.
  110: *> \endverbatim
  111: *>
  112: *> \param[in,out] NDIV
  113: *> \verbatim
  114: *>          NDIV is INTEGER
  115: *>         Increment NDIV by 1 for each division.
  116: *> \endverbatim
  117: *>
  118: *> \param[in] IEEE
  119: *> \verbatim
  120: *>          IEEE is LOGICAL
  121: *>         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
  122: *> \endverbatim
  123: *>
  124: *> \param[in,out] TTYPE
  125: *> \verbatim
  126: *>          TTYPE is INTEGER
  127: *>         Shift type.
  128: *> \endverbatim
  129: *>
  130: *> \param[in,out] DMIN1
  131: *> \verbatim
  132: *>          DMIN1 is DOUBLE PRECISION
  133: *> \endverbatim
  134: *>
  135: *> \param[in,out] DMIN2
  136: *> \verbatim
  137: *>          DMIN2 is DOUBLE PRECISION
  138: *> \endverbatim
  139: *>
  140: *> \param[in,out] DN
  141: *> \verbatim
  142: *>          DN is DOUBLE PRECISION
  143: *> \endverbatim
  144: *>
  145: *> \param[in,out] DN1
  146: *> \verbatim
  147: *>          DN1 is DOUBLE PRECISION
  148: *> \endverbatim
  149: *>
  150: *> \param[in,out] DN2
  151: *> \verbatim
  152: *>          DN2 is DOUBLE PRECISION
  153: *> \endverbatim
  154: *>
  155: *> \param[in,out] G
  156: *> \verbatim
  157: *>          G is DOUBLE PRECISION
  158: *> \endverbatim
  159: *>
  160: *> \param[in,out] TAU
  161: *> \verbatim
  162: *>          TAU is DOUBLE PRECISION
  163: *>
  164: *>         These are passed as arguments in order to save their values
  165: *>         between calls to DLASQ3.
  166: *> \endverbatim
  167: *
  168: *  Authors:
  169: *  ========
  170: *
  171: *> \author Univ. of Tennessee 
  172: *> \author Univ. of California Berkeley 
  173: *> \author Univ. of Colorado Denver 
  174: *> \author NAG Ltd. 
  175: *
  176: *> \date June 2016
  177: *
  178: *> \ingroup auxOTHERcomputational
  179: *
  180: *  =====================================================================
  181:       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
  182:      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
  183:      $                   DN2, G, TAU )
  184: *
  185: *  -- LAPACK computational routine (version 3.6.1) --
  186: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  187: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  188: *     June 2016
  189: *
  190: *     .. Scalar Arguments ..
  191:       LOGICAL            IEEE
  192:       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
  193:       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
  194:      $                   QMAX, SIGMA, TAU
  195: *     ..
  196: *     .. Array Arguments ..
  197:       DOUBLE PRECISION   Z( * )
  198: *     ..
  199: *
  200: *  =====================================================================
  201: *
  202: *     .. Parameters ..
  203:       DOUBLE PRECISION   CBIAS
  204:       PARAMETER          ( CBIAS = 1.50D0 )
  205:       DOUBLE PRECISION   ZERO, QURTR, HALF, ONE, TWO, HUNDRD
  206:       PARAMETER          ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
  207:      $                     ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
  208: *     ..
  209: *     .. Local Scalars ..
  210:       INTEGER            IPN4, J4, N0IN, NN, TTYPE
  211:       DOUBLE PRECISION   EPS, S, T, TEMP, TOL, TOL2
  212: *     ..
  213: *     .. External Subroutines ..
  214:       EXTERNAL           DLASQ4, DLASQ5, DLASQ6
  215: *     ..
  216: *     .. External Function ..
  217:       DOUBLE PRECISION   DLAMCH
  218:       LOGICAL            DISNAN
  219:       EXTERNAL           DISNAN, DLAMCH
  220: *     ..
  221: *     .. Intrinsic Functions ..
  222:       INTRINSIC          ABS, MAX, MIN, SQRT
  223: *     ..
  224: *     .. Executable Statements ..
  225: *
  226:       N0IN = N0
  227:       EPS = DLAMCH( 'Precision' )
  228:       TOL = EPS*HUNDRD
  229:       TOL2 = TOL**2
  230: *
  231: *     Check for deflation.
  232: *
  233:    10 CONTINUE
  234: *
  235:       IF( N0.LT.I0 )
  236:      $   RETURN
  237:       IF( N0.EQ.I0 )
  238:      $   GO TO 20
  239:       NN = 4*N0 + PP
  240:       IF( N0.EQ.( I0+1 ) )
  241:      $   GO TO 40
  242: *
  243: *     Check whether E(N0-1) is negligible, 1 eigenvalue.
  244: *
  245:       IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
  246:      $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
  247:      $   GO TO 30
  248: *
  249:    20 CONTINUE
  250: *
  251:       Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
  252:       N0 = N0 - 1
  253:       GO TO 10
  254: *
  255: *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
  256: *
  257:    30 CONTINUE
  258: *
  259:       IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
  260:      $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
  261:      $   GO TO 50
  262: *
  263:    40 CONTINUE
  264: *
  265:       IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
  266:          S = Z( NN-3 )
  267:          Z( NN-3 ) = Z( NN-7 )
  268:          Z( NN-7 ) = S
  269:       END IF
  270:       T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
  271:       IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN
  272:          S = Z( NN-3 )*( Z( NN-5 ) / T )
  273:          IF( S.LE.T ) THEN
  274:             S = Z( NN-3 )*( Z( NN-5 ) /
  275:      $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
  276:          ELSE
  277:             S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
  278:          END IF
  279:          T = Z( NN-7 ) + ( S+Z( NN-5 ) )
  280:          Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
  281:          Z( NN-7 ) = T
  282:       END IF
  283:       Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
  284:       Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
  285:       N0 = N0 - 2
  286:       GO TO 10
  287: *
  288:    50 CONTINUE
  289:       IF( PP.EQ.2 ) 
  290:      $   PP = 0
  291: *
  292: *     Reverse the qd-array, if warranted.
  293: *
  294:       IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
  295:          IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
  296:             IPN4 = 4*( I0+N0 )
  297:             DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
  298:                TEMP = Z( J4-3 )
  299:                Z( J4-3 ) = Z( IPN4-J4-3 )
  300:                Z( IPN4-J4-3 ) = TEMP
  301:                TEMP = Z( J4-2 )
  302:                Z( J4-2 ) = Z( IPN4-J4-2 )
  303:                Z( IPN4-J4-2 ) = TEMP
  304:                TEMP = Z( J4-1 )
  305:                Z( J4-1 ) = Z( IPN4-J4-5 )
  306:                Z( IPN4-J4-5 ) = TEMP
  307:                TEMP = Z( J4 )
  308:                Z( J4 ) = Z( IPN4-J4-4 )
  309:                Z( IPN4-J4-4 ) = TEMP
  310:    60       CONTINUE
  311:             IF( N0-I0.LE.4 ) THEN
  312:                Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
  313:                Z( 4*N0-PP ) = Z( 4*I0-PP )
  314:             END IF
  315:             DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
  316:             Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
  317:      $                            Z( 4*I0+PP+3 ) )
  318:             Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
  319:      $                          Z( 4*I0-PP+4 ) )
  320:             QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
  321:             DMIN = -ZERO
  322:          END IF
  323:       END IF
  324: *
  325: *     Choose a shift.
  326: *
  327:       CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
  328:      $             DN2, TAU, TTYPE, G )
  329: *
  330: *     Call dqds until DMIN > 0.
  331: *
  332:    70 CONTINUE
  333: *
  334:       CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
  335:      $             DN1, DN2, IEEE, EPS )
  336: *
  337:       NDIV = NDIV + ( N0-I0+2 )
  338:       ITER = ITER + 1
  339: *
  340: *     Check status.
  341: *
  342:       IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
  343: *
  344: *        Success.
  345: *
  346:          GO TO 90
  347: *
  348:       ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 
  349:      $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
  350:      $         ABS( DN ).LT.TOL*SIGMA ) THEN
  351: *
  352: *        Convergence hidden by negative DN.
  353: *
  354:          Z( 4*( N0-1 )-PP+2 ) = ZERO
  355:          DMIN = ZERO
  356:          GO TO 90
  357:       ELSE IF( DMIN.LT.ZERO ) THEN
  358: *
  359: *        TAU too big. Select new TAU and try again.
  360: *
  361:          NFAIL = NFAIL + 1
  362:          IF( TTYPE.LT.-22 ) THEN
  363: *
  364: *           Failed twice. Play it safe.
  365: *
  366:             TAU = ZERO
  367:          ELSE IF( DMIN1.GT.ZERO ) THEN
  368: *
  369: *           Late failure. Gives excellent shift.
  370: *
  371:             TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
  372:             TTYPE = TTYPE - 11
  373:          ELSE
  374: *
  375: *           Early failure. Divide by 4.
  376: *
  377:             TAU = QURTR*TAU
  378:             TTYPE = TTYPE - 12
  379:          END IF
  380:          GO TO 70
  381:       ELSE IF( DISNAN( DMIN ) ) THEN
  382: *
  383: *        NaN.
  384: *
  385:          IF( TAU.EQ.ZERO ) THEN
  386:             GO TO 80
  387:          ELSE
  388:             TAU = ZERO
  389:             GO TO 70
  390:          END IF
  391:       ELSE
  392: *            
  393: *        Possible underflow. Play it safe.
  394: *
  395:          GO TO 80
  396:       END IF
  397: *
  398: *     Risk of underflow.
  399: *
  400:    80 CONTINUE
  401:       CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
  402:       NDIV = NDIV + ( N0-I0+2 )
  403:       ITER = ITER + 1
  404:       TAU = ZERO
  405: *
  406:    90 CONTINUE
  407:       IF( TAU.LT.SIGMA ) THEN
  408:          DESIG = DESIG + TAU
  409:          T = SIGMA + DESIG
  410:          DESIG = DESIG - ( T-SIGMA )
  411:       ELSE
  412:          T = SIGMA + TAU
  413:          DESIG = SIGMA - ( T-TAU ) + DESIG
  414:       END IF
  415:       SIGMA = T
  416: *
  417:       RETURN
  418: *
  419: *     End of DLASQ3
  420: *
  421:       END