Annotation of rpl/lapack/lapack/dlasq4.f, revision 1.21

1.12      bertrand    1: *> \brief \b DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.
1.9       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.17      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.17      bertrand    9: *> Download DLASQ4 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq4.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq4.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq4.f">
1.9       bertrand   15: *> [TXT]</a>
1.17      bertrand   16: *> \endhtmlonly
1.9       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
                     22: *                          DN1, DN2, TAU, TTYPE, G )
1.17      bertrand   23: *
1.9       bertrand   24: *       .. Scalar Arguments ..
                     25: *       INTEGER            I0, N0, N0IN, PP, TTYPE
                     26: *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
                     27: *       ..
                     28: *       .. Array Arguments ..
                     29: *       DOUBLE PRECISION   Z( * )
                     30: *       ..
1.17      bertrand   31: *
1.9       bertrand   32: *
                     33: *> \par Purpose:
                     34: *  =============
                     35: *>
                     36: *> \verbatim
                     37: *>
                     38: *> DLASQ4 computes an approximation TAU to the smallest eigenvalue
                     39: *> using values of d from the previous transform.
                     40: *> \endverbatim
                     41: *
                     42: *  Arguments:
                     43: *  ==========
                     44: *
                     45: *> \param[in] I0
                     46: *> \verbatim
                     47: *>          I0 is INTEGER
                     48: *>        First index.
                     49: *> \endverbatim
                     50: *>
                     51: *> \param[in] N0
                     52: *> \verbatim
                     53: *>          N0 is INTEGER
                     54: *>        Last index.
                     55: *> \endverbatim
                     56: *>
                     57: *> \param[in] Z
                     58: *> \verbatim
1.15      bertrand   59: *>          Z is DOUBLE PRECISION array, dimension ( 4*N0 )
1.9       bertrand   60: *>        Z holds the qd array.
                     61: *> \endverbatim
                     62: *>
                     63: *> \param[in] PP
                     64: *> \verbatim
                     65: *>          PP is INTEGER
                     66: *>        PP=0 for ping, PP=1 for pong.
                     67: *> \endverbatim
                     68: *>
                     69: *> \param[in] N0IN
                     70: *> \verbatim
                     71: *>          N0IN is INTEGER
                     72: *>        The value of N0 at start of EIGTEST.
                     73: *> \endverbatim
                     74: *>
                     75: *> \param[in] DMIN
                     76: *> \verbatim
                     77: *>          DMIN is DOUBLE PRECISION
                     78: *>        Minimum value of d.
                     79: *> \endverbatim
                     80: *>
                     81: *> \param[in] DMIN1
                     82: *> \verbatim
                     83: *>          DMIN1 is DOUBLE PRECISION
                     84: *>        Minimum value of d, excluding D( N0 ).
                     85: *> \endverbatim
                     86: *>
                     87: *> \param[in] DMIN2
                     88: *> \verbatim
                     89: *>          DMIN2 is DOUBLE PRECISION
                     90: *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
                     91: *> \endverbatim
                     92: *>
                     93: *> \param[in] DN
                     94: *> \verbatim
                     95: *>          DN is DOUBLE PRECISION
                     96: *>        d(N)
                     97: *> \endverbatim
                     98: *>
                     99: *> \param[in] DN1
                    100: *> \verbatim
                    101: *>          DN1 is DOUBLE PRECISION
                    102: *>        d(N-1)
                    103: *> \endverbatim
                    104: *>
                    105: *> \param[in] DN2
                    106: *> \verbatim
                    107: *>          DN2 is DOUBLE PRECISION
                    108: *>        d(N-2)
                    109: *> \endverbatim
                    110: *>
                    111: *> \param[out] TAU
                    112: *> \verbatim
                    113: *>          TAU is DOUBLE PRECISION
                    114: *>        This is the shift.
                    115: *> \endverbatim
                    116: *>
                    117: *> \param[out] TTYPE
                    118: *> \verbatim
                    119: *>          TTYPE is INTEGER
                    120: *>        Shift type.
                    121: *> \endverbatim
                    122: *>
                    123: *> \param[in,out] G
                    124: *> \verbatim
1.15      bertrand  125: *>          G is DOUBLE PRECISION
1.9       bertrand  126: *>        G is passed as an argument in order to save its value between
                    127: *>        calls to DLASQ4.
                    128: *> \endverbatim
                    129: *
                    130: *  Authors:
                    131: *  ========
                    132: *
1.17      bertrand  133: *> \author Univ. of Tennessee
                    134: *> \author Univ. of California Berkeley
                    135: *> \author Univ. of Colorado Denver
                    136: *> \author NAG Ltd.
1.9       bertrand  137: *
                    138: *> \ingroup auxOTHERcomputational
                    139: *
                    140: *> \par Further Details:
                    141: *  =====================
                    142: *>
                    143: *> \verbatim
                    144: *>
                    145: *>  CNST1 = 9/16
                    146: *> \endverbatim
                    147: *>
                    148: *  =====================================================================
1.1       bertrand  149:       SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
                    150:      $                   DN1, DN2, TAU, TTYPE, G )
                    151: *
1.21    ! bertrand  152: *  -- LAPACK computational routine --
1.1       bertrand  153: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    154: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    155: *
                    156: *     .. Scalar Arguments ..
                    157:       INTEGER            I0, N0, N0IN, PP, TTYPE
                    158:       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
                    159: *     ..
                    160: *     .. Array Arguments ..
                    161:       DOUBLE PRECISION   Z( * )
                    162: *     ..
                    163: *
                    164: *  =====================================================================
                    165: *
                    166: *     .. Parameters ..
                    167:       DOUBLE PRECISION   CNST1, CNST2, CNST3
                    168:       PARAMETER          ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
                    169:      $                   CNST3 = 1.050D0 )
                    170:       DOUBLE PRECISION   QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
                    171:       PARAMETER          ( QURTR = 0.250D0, THIRD = 0.3330D0,
                    172:      $                   HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
                    173:      $                   TWO = 2.0D0, HUNDRD = 100.0D0 )
                    174: *     ..
                    175: *     .. Local Scalars ..
                    176:       INTEGER            I4, NN, NP
                    177:       DOUBLE PRECISION   A2, B1, B2, GAM, GAP1, GAP2, S
                    178: *     ..
                    179: *     .. Intrinsic Functions ..
                    180:       INTRINSIC          MAX, MIN, SQRT
                    181: *     ..
                    182: *     .. Executable Statements ..
                    183: *
                    184: *     A negative DMIN forces the shift to take that absolute value
                    185: *     TTYPE records the type of shift.
                    186: *
                    187:       IF( DMIN.LE.ZERO ) THEN
                    188:          TAU = -DMIN
                    189:          TTYPE = -1
                    190:          RETURN
                    191:       END IF
1.17      bertrand  192: *
1.1       bertrand  193:       NN = 4*N0 + PP
                    194:       IF( N0IN.EQ.N0 ) THEN
                    195: *
                    196: *        No eigenvalues deflated.
                    197: *
                    198:          IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
                    199: *
                    200:             B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
                    201:             B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
                    202:             A2 = Z( NN-7 ) + Z( NN-5 )
                    203: *
                    204: *           Cases 2 and 3.
                    205: *
                    206:             IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
                    207:                GAP2 = DMIN2 - A2 - DMIN2*QURTR
                    208:                IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
                    209:                   GAP1 = A2 - DN - ( B2 / GAP2 )*B2
                    210:                ELSE
                    211:                   GAP1 = A2 - DN - ( B1+B2 )
                    212:                END IF
                    213:                IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
                    214:                   S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
                    215:                   TTYPE = -2
                    216:                ELSE
                    217:                   S = ZERO
                    218:                   IF( DN.GT.B1 )
                    219:      $               S = DN - B1
                    220:                   IF( A2.GT.( B1+B2 ) )
                    221:      $               S = MIN( S, A2-( B1+B2 ) )
                    222:                   S = MAX( S, THIRD*DMIN )
                    223:                   TTYPE = -3
                    224:                END IF
                    225:             ELSE
                    226: *
                    227: *              Case 4.
                    228: *
                    229:                TTYPE = -4
                    230:                S = QURTR*DMIN
                    231:                IF( DMIN.EQ.DN ) THEN
                    232:                   GAM = DN
                    233:                   A2 = ZERO
                    234:                   IF( Z( NN-5 ) .GT. Z( NN-7 ) )
                    235:      $               RETURN
                    236:                   B2 = Z( NN-5 ) / Z( NN-7 )
                    237:                   NP = NN - 9
                    238:                ELSE
                    239:                   NP = NN - 2*PP
                    240:                   GAM = DN1
                    241:                   IF( Z( NP-4 ) .GT. Z( NP-2 ) )
                    242:      $               RETURN
                    243:                   A2 = Z( NP-4 ) / Z( NP-2 )
                    244:                   IF( Z( NN-9 ) .GT. Z( NN-11 ) )
                    245:      $               RETURN
                    246:                   B2 = Z( NN-9 ) / Z( NN-11 )
                    247:                   NP = NN - 13
                    248:                END IF
                    249: *
                    250: *              Approximate contribution to norm squared from I < NN-1.
                    251: *
                    252:                A2 = A2 + B2
                    253:                DO 10 I4 = NP, 4*I0 - 1 + PP, -4
                    254:                   IF( B2.EQ.ZERO )
                    255:      $               GO TO 20
                    256:                   B1 = B2
                    257:                   IF( Z( I4 ) .GT. Z( I4-2 ) )
                    258:      $               RETURN
                    259:                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
                    260:                   A2 = A2 + B2
1.17      bertrand  261:                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
1.1       bertrand  262:      $               GO TO 20
                    263:    10          CONTINUE
                    264:    20          CONTINUE
                    265:                A2 = CNST3*A2
                    266: *
                    267: *              Rayleigh quotient residual bound.
                    268: *
                    269:                IF( A2.LT.CNST1 )
                    270:      $            S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
                    271:             END IF
                    272:          ELSE IF( DMIN.EQ.DN2 ) THEN
                    273: *
                    274: *           Case 5.
                    275: *
                    276:             TTYPE = -5
                    277:             S = QURTR*DMIN
                    278: *
                    279: *           Compute contribution to norm squared from I > NN-2.
                    280: *
                    281:             NP = NN - 2*PP
                    282:             B1 = Z( NP-2 )
                    283:             B2 = Z( NP-6 )
                    284:             GAM = DN2
                    285:             IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
                    286:      $         RETURN
                    287:             A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
                    288: *
                    289: *           Approximate contribution to norm squared from I < NN-2.
                    290: *
                    291:             IF( N0-I0.GT.2 ) THEN
                    292:                B2 = Z( NN-13 ) / Z( NN-15 )
                    293:                A2 = A2 + B2
                    294:                DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
                    295:                   IF( B2.EQ.ZERO )
                    296:      $               GO TO 40
                    297:                   B1 = B2
                    298:                   IF( Z( I4 ) .GT. Z( I4-2 ) )
                    299:      $               RETURN
                    300:                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
                    301:                   A2 = A2 + B2
1.17      bertrand  302:                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
1.1       bertrand  303:      $               GO TO 40
                    304:    30          CONTINUE
                    305:    40          CONTINUE
                    306:                A2 = CNST3*A2
                    307:             END IF
                    308: *
                    309:             IF( A2.LT.CNST1 )
                    310:      $         S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
                    311:          ELSE
                    312: *
                    313: *           Case 6, no information to guide us.
                    314: *
                    315:             IF( TTYPE.EQ.-6 ) THEN
                    316:                G = G + THIRD*( ONE-G )
                    317:             ELSE IF( TTYPE.EQ.-18 ) THEN
                    318:                G = QURTR*THIRD
                    319:             ELSE
                    320:                G = QURTR
                    321:             END IF
                    322:             S = G*DMIN
                    323:             TTYPE = -6
                    324:          END IF
                    325: *
                    326:       ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
                    327: *
                    328: *        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
                    329: *
1.17      bertrand  330:          IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
1.1       bertrand  331: *
                    332: *           Cases 7 and 8.
                    333: *
                    334:             TTYPE = -7
                    335:             S = THIRD*DMIN1
                    336:             IF( Z( NN-5 ).GT.Z( NN-7 ) )
                    337:      $         RETURN
                    338:             B1 = Z( NN-5 ) / Z( NN-7 )
                    339:             B2 = B1
                    340:             IF( B2.EQ.ZERO )
                    341:      $         GO TO 60
                    342:             DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
                    343:                A2 = B1
                    344:                IF( Z( I4 ).GT.Z( I4-2 ) )
                    345:      $            RETURN
                    346:                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
                    347:                B2 = B2 + B1
1.17      bertrand  348:                IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
1.1       bertrand  349:      $            GO TO 60
                    350:    50       CONTINUE
                    351:    60       CONTINUE
                    352:             B2 = SQRT( CNST3*B2 )
                    353:             A2 = DMIN1 / ( ONE+B2**2 )
                    354:             GAP2 = HALF*DMIN2 - A2
                    355:             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
                    356:                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
1.17      bertrand  357:             ELSE
1.1       bertrand  358:                S = MAX( S, A2*( ONE-CNST2*B2 ) )
                    359:                TTYPE = -8
                    360:             END IF
                    361:          ELSE
                    362: *
                    363: *           Case 9.
                    364: *
                    365:             S = QURTR*DMIN1
                    366:             IF( DMIN1.EQ.DN1 )
                    367:      $         S = HALF*DMIN1
                    368:             TTYPE = -9
                    369:          END IF
                    370: *
                    371:       ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
                    372: *
                    373: *        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
                    374: *
                    375: *        Cases 10 and 11.
                    376: *
1.17      bertrand  377:          IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
1.1       bertrand  378:             TTYPE = -10
                    379:             S = THIRD*DMIN2
                    380:             IF( Z( NN-5 ).GT.Z( NN-7 ) )
                    381:      $         RETURN
                    382:             B1 = Z( NN-5 ) / Z( NN-7 )
                    383:             B2 = B1
                    384:             IF( B2.EQ.ZERO )
                    385:      $         GO TO 80
                    386:             DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
                    387:                IF( Z( I4 ).GT.Z( I4-2 ) )
                    388:      $            RETURN
                    389:                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
                    390:                B2 = B2 + B1
                    391:                IF( HUNDRD*B1.LT.B2 )
                    392:      $            GO TO 80
                    393:    70       CONTINUE
                    394:    80       CONTINUE
                    395:             B2 = SQRT( CNST3*B2 )
                    396:             A2 = DMIN2 / ( ONE+B2**2 )
                    397:             GAP2 = Z( NN-7 ) + Z( NN-9 ) -
                    398:      $             SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
                    399:             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
                    400:                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
1.17      bertrand  401:             ELSE
1.1       bertrand  402:                S = MAX( S, A2*( ONE-CNST2*B2 ) )
                    403:             END IF
                    404:          ELSE
                    405:             S = QURTR*DMIN2
                    406:             TTYPE = -11
                    407:          END IF
                    408:       ELSE IF( N0IN.GT.( N0+2 ) ) THEN
                    409: *
                    410: *        Case 12, more than two eigenvalues deflated. No information.
                    411: *
1.17      bertrand  412:          S = ZERO
1.1       bertrand  413:          TTYPE = -12
                    414:       END IF
                    415: *
                    416:       TAU = S
                    417:       RETURN
                    418: *
                    419: *     End of DLASQ4
                    420: *
                    421:       END

CVSweb interface <joel.bertrand@systella.fr>