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

1.1       bertrand    1:       SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
                      2:      $                   DN1, DN2, TAU, TTYPE, G )
                      3: *
1.8     ! bertrand    4: *  -- LAPACK routine (version 3.3.1)                                    --
1.1       bertrand    5: *
                      6: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
                      7: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
                      8: *  -- Berkeley                                                        --
                      9: *  -- November 2008                                                   --
                     10: *
                     11: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                     12: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                     13: *
                     14: *     .. Scalar Arguments ..
                     15:       INTEGER            I0, N0, N0IN, PP, TTYPE
                     16:       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
                     17: *     ..
                     18: *     .. Array Arguments ..
                     19:       DOUBLE PRECISION   Z( * )
                     20: *     ..
                     21: *
                     22: *  Purpose
                     23: *  =======
                     24: *
                     25: *  DLASQ4 computes an approximation TAU to the smallest eigenvalue
                     26: *  using values of d from the previous transform.
                     27: *
1.8     ! bertrand   28: *  Arguments
        !            29: *  =========
        !            30: *
1.1       bertrand   31: *  I0    (input) INTEGER
                     32: *        First index.
                     33: *
                     34: *  N0    (input) INTEGER
                     35: *        Last index.
                     36: *
                     37: *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
                     38: *        Z holds the qd array.
                     39: *
                     40: *  PP    (input) INTEGER
                     41: *        PP=0 for ping, PP=1 for pong.
                     42: *
                     43: *  NOIN  (input) INTEGER
                     44: *        The value of N0 at start of EIGTEST.
                     45: *
                     46: *  DMIN  (input) DOUBLE PRECISION
                     47: *        Minimum value of d.
                     48: *
                     49: *  DMIN1 (input) DOUBLE PRECISION
                     50: *        Minimum value of d, excluding D( N0 ).
                     51: *
                     52: *  DMIN2 (input) DOUBLE PRECISION
                     53: *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
                     54: *
                     55: *  DN    (input) DOUBLE PRECISION
                     56: *        d(N)
                     57: *
                     58: *  DN1   (input) DOUBLE PRECISION
                     59: *        d(N-1)
                     60: *
                     61: *  DN2   (input) DOUBLE PRECISION
                     62: *        d(N-2)
                     63: *
                     64: *  TAU   (output) DOUBLE PRECISION
                     65: *        This is the shift.
                     66: *
                     67: *  TTYPE (output) INTEGER
                     68: *        Shift type.
                     69: *
                     70: *  G     (input/output) REAL
                     71: *        G is passed as an argument in order to save its value between
                     72: *        calls to DLASQ4.
                     73: *
                     74: *  Further Details
                     75: *  ===============
                     76: *  CNST1 = 9/16
                     77: *
                     78: *  =====================================================================
                     79: *
                     80: *     .. Parameters ..
                     81:       DOUBLE PRECISION   CNST1, CNST2, CNST3
                     82:       PARAMETER          ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
                     83:      $                   CNST3 = 1.050D0 )
                     84:       DOUBLE PRECISION   QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
                     85:       PARAMETER          ( QURTR = 0.250D0, THIRD = 0.3330D0,
                     86:      $                   HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
                     87:      $                   TWO = 2.0D0, HUNDRD = 100.0D0 )
                     88: *     ..
                     89: *     .. Local Scalars ..
                     90:       INTEGER            I4, NN, NP
                     91:       DOUBLE PRECISION   A2, B1, B2, GAM, GAP1, GAP2, S
                     92: *     ..
                     93: *     .. Intrinsic Functions ..
                     94:       INTRINSIC          MAX, MIN, SQRT
                     95: *     ..
                     96: *     .. Executable Statements ..
                     97: *
                     98: *     A negative DMIN forces the shift to take that absolute value
                     99: *     TTYPE records the type of shift.
                    100: *
                    101:       IF( DMIN.LE.ZERO ) THEN
                    102:          TAU = -DMIN
                    103:          TTYPE = -1
                    104:          RETURN
                    105:       END IF
                    106: *       
                    107:       NN = 4*N0 + PP
                    108:       IF( N0IN.EQ.N0 ) THEN
                    109: *
                    110: *        No eigenvalues deflated.
                    111: *
                    112:          IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
                    113: *
                    114:             B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
                    115:             B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
                    116:             A2 = Z( NN-7 ) + Z( NN-5 )
                    117: *
                    118: *           Cases 2 and 3.
                    119: *
                    120:             IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
                    121:                GAP2 = DMIN2 - A2 - DMIN2*QURTR
                    122:                IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
                    123:                   GAP1 = A2 - DN - ( B2 / GAP2 )*B2
                    124:                ELSE
                    125:                   GAP1 = A2 - DN - ( B1+B2 )
                    126:                END IF
                    127:                IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
                    128:                   S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
                    129:                   TTYPE = -2
                    130:                ELSE
                    131:                   S = ZERO
                    132:                   IF( DN.GT.B1 )
                    133:      $               S = DN - B1
                    134:                   IF( A2.GT.( B1+B2 ) )
                    135:      $               S = MIN( S, A2-( B1+B2 ) )
                    136:                   S = MAX( S, THIRD*DMIN )
                    137:                   TTYPE = -3
                    138:                END IF
                    139:             ELSE
                    140: *
                    141: *              Case 4.
                    142: *
                    143:                TTYPE = -4
                    144:                S = QURTR*DMIN
                    145:                IF( DMIN.EQ.DN ) THEN
                    146:                   GAM = DN
                    147:                   A2 = ZERO
                    148:                   IF( Z( NN-5 ) .GT. Z( NN-7 ) )
                    149:      $               RETURN
                    150:                   B2 = Z( NN-5 ) / Z( NN-7 )
                    151:                   NP = NN - 9
                    152:                ELSE
                    153:                   NP = NN - 2*PP
                    154:                   B2 = Z( NP-2 )
                    155:                   GAM = DN1
                    156:                   IF( Z( NP-4 ) .GT. Z( NP-2 ) )
                    157:      $               RETURN
                    158:                   A2 = Z( NP-4 ) / Z( NP-2 )
                    159:                   IF( Z( NN-9 ) .GT. Z( NN-11 ) )
                    160:      $               RETURN
                    161:                   B2 = Z( NN-9 ) / Z( NN-11 )
                    162:                   NP = NN - 13
                    163:                END IF
                    164: *
                    165: *              Approximate contribution to norm squared from I < NN-1.
                    166: *
                    167:                A2 = A2 + B2
                    168:                DO 10 I4 = NP, 4*I0 - 1 + PP, -4
                    169:                   IF( B2.EQ.ZERO )
                    170:      $               GO TO 20
                    171:                   B1 = B2
                    172:                   IF( Z( I4 ) .GT. Z( I4-2 ) )
                    173:      $               RETURN
                    174:                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
                    175:                   A2 = A2 + B2
                    176:                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 
                    177:      $               GO TO 20
                    178:    10          CONTINUE
                    179:    20          CONTINUE
                    180:                A2 = CNST3*A2
                    181: *
                    182: *              Rayleigh quotient residual bound.
                    183: *
                    184:                IF( A2.LT.CNST1 )
                    185:      $            S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
                    186:             END IF
                    187:          ELSE IF( DMIN.EQ.DN2 ) THEN
                    188: *
                    189: *           Case 5.
                    190: *
                    191:             TTYPE = -5
                    192:             S = QURTR*DMIN
                    193: *
                    194: *           Compute contribution to norm squared from I > NN-2.
                    195: *
                    196:             NP = NN - 2*PP
                    197:             B1 = Z( NP-2 )
                    198:             B2 = Z( NP-6 )
                    199:             GAM = DN2
                    200:             IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
                    201:      $         RETURN
                    202:             A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
                    203: *
                    204: *           Approximate contribution to norm squared from I < NN-2.
                    205: *
                    206:             IF( N0-I0.GT.2 ) THEN
                    207:                B2 = Z( NN-13 ) / Z( NN-15 )
                    208:                A2 = A2 + B2
                    209:                DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
                    210:                   IF( B2.EQ.ZERO )
                    211:      $               GO TO 40
                    212:                   B1 = B2
                    213:                   IF( Z( I4 ) .GT. Z( I4-2 ) )
                    214:      $               RETURN
                    215:                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
                    216:                   A2 = A2 + B2
                    217:                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 
                    218:      $               GO TO 40
                    219:    30          CONTINUE
                    220:    40          CONTINUE
                    221:                A2 = CNST3*A2
                    222:             END IF
                    223: *
                    224:             IF( A2.LT.CNST1 )
                    225:      $         S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
                    226:          ELSE
                    227: *
                    228: *           Case 6, no information to guide us.
                    229: *
                    230:             IF( TTYPE.EQ.-6 ) THEN
                    231:                G = G + THIRD*( ONE-G )
                    232:             ELSE IF( TTYPE.EQ.-18 ) THEN
                    233:                G = QURTR*THIRD
                    234:             ELSE
                    235:                G = QURTR
                    236:             END IF
                    237:             S = G*DMIN
                    238:             TTYPE = -6
                    239:          END IF
                    240: *
                    241:       ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
                    242: *
                    243: *        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
                    244: *
                    245:          IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN 
                    246: *
                    247: *           Cases 7 and 8.
                    248: *
                    249:             TTYPE = -7
                    250:             S = THIRD*DMIN1
                    251:             IF( Z( NN-5 ).GT.Z( NN-7 ) )
                    252:      $         RETURN
                    253:             B1 = Z( NN-5 ) / Z( NN-7 )
                    254:             B2 = B1
                    255:             IF( B2.EQ.ZERO )
                    256:      $         GO TO 60
                    257:             DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
                    258:                A2 = B1
                    259:                IF( Z( I4 ).GT.Z( I4-2 ) )
                    260:      $            RETURN
                    261:                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
                    262:                B2 = B2 + B1
                    263:                IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) 
                    264:      $            GO TO 60
                    265:    50       CONTINUE
                    266:    60       CONTINUE
                    267:             B2 = SQRT( CNST3*B2 )
                    268:             A2 = DMIN1 / ( ONE+B2**2 )
                    269:             GAP2 = HALF*DMIN2 - A2
                    270:             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
                    271:                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
                    272:             ELSE 
                    273:                S = MAX( S, A2*( ONE-CNST2*B2 ) )
                    274:                TTYPE = -8
                    275:             END IF
                    276:          ELSE
                    277: *
                    278: *           Case 9.
                    279: *
                    280:             S = QURTR*DMIN1
                    281:             IF( DMIN1.EQ.DN1 )
                    282:      $         S = HALF*DMIN1
                    283:             TTYPE = -9
                    284:          END IF
                    285: *
                    286:       ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
                    287: *
                    288: *        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
                    289: *
                    290: *        Cases 10 and 11.
                    291: *
                    292:          IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN 
                    293:             TTYPE = -10
                    294:             S = THIRD*DMIN2
                    295:             IF( Z( NN-5 ).GT.Z( NN-7 ) )
                    296:      $         RETURN
                    297:             B1 = Z( NN-5 ) / Z( NN-7 )
                    298:             B2 = B1
                    299:             IF( B2.EQ.ZERO )
                    300:      $         GO TO 80
                    301:             DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
                    302:                IF( Z( I4 ).GT.Z( I4-2 ) )
                    303:      $            RETURN
                    304:                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
                    305:                B2 = B2 + B1
                    306:                IF( HUNDRD*B1.LT.B2 )
                    307:      $            GO TO 80
                    308:    70       CONTINUE
                    309:    80       CONTINUE
                    310:             B2 = SQRT( CNST3*B2 )
                    311:             A2 = DMIN2 / ( ONE+B2**2 )
                    312:             GAP2 = Z( NN-7 ) + Z( NN-9 ) -
                    313:      $             SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
                    314:             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
                    315:                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
                    316:             ELSE 
                    317:                S = MAX( S, A2*( ONE-CNST2*B2 ) )
                    318:             END IF
                    319:          ELSE
                    320:             S = QURTR*DMIN2
                    321:             TTYPE = -11
                    322:          END IF
                    323:       ELSE IF( N0IN.GT.( N0+2 ) ) THEN
                    324: *
                    325: *        Case 12, more than two eigenvalues deflated. No information.
                    326: *
                    327:          S = ZERO 
                    328:          TTYPE = -12
                    329:       END IF
                    330: *
                    331:       TAU = S
                    332:       RETURN
                    333: *
                    334: *     End of DLASQ4
                    335: *
                    336:       END

CVSweb interface <joel.bertrand@systella.fr>