Return to dlarra.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dlarra.f, revision 1.3

1.1       bertrand    1:       SUBROUTINE DLARRA( N, D, E, E2, SPLTOL, TNRM,
                      2:      $                    NSPLIT, ISPLIT, INFO )
                      3:       IMPLICIT NONE
                      4: *
                      5: *  -- LAPACK auxiliary routine (version 3.2) --
                      6: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      7: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      8: *     November 2006
                      9: *
                     10: *     .. Scalar Arguments ..
                     11:       INTEGER            INFO, N, NSPLIT
                     12:       DOUBLE PRECISION    SPLTOL, TNRM
                     13: *     ..
                     14: *     .. Array Arguments ..
                     15:       INTEGER            ISPLIT( * )
                     16:       DOUBLE PRECISION   D( * ), E( * ), E2( * )
                     17: *     ..
                     18: *
                     19: *  Purpose
                     20: *  =======
                     21: *
                     22: *  Compute the splitting points with threshold SPLTOL.
                     23: *  DLARRA sets any "small" off-diagonal elements to zero.
                     24: *
                     25: *  Arguments
                     26: *  =========
                     27: *
                     28: *  N       (input) INTEGER
                     29: *          The order of the matrix. N > 0.
                     30: *
                     31: *  D       (input) DOUBLE PRECISION array, dimension (N)
                     32: *          On entry, the N diagonal elements of the tridiagonal
                     33: *          matrix T.
                     34: *
                     35: *  E       (input/output) DOUBLE PRECISION array, dimension (N)
                     36: *          On entry, the first (N-1) entries contain the subdiagonal
                     37: *          elements of the tridiagonal matrix T; E(N) need not be set.
                     38: *          On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT,
                     39: *          are set to zero, the other entries of E are untouched.
                     40: *
                     41: *  E2      (input/output) DOUBLE PRECISION array, dimension (N)
                     42: *          On entry, the first (N-1) entries contain the SQUARES of the
                     43: *          subdiagonal elements of the tridiagonal matrix T;
                     44: *          E2(N) need not be set.
                     45: *          On exit, the entries E2( ISPLIT( I ) ),
                     46: *          1 <= I <= NSPLIT, have been set to zero
                     47: *
                     48: *  SPLTOL (input) DOUBLE PRECISION
                     49: *          The threshold for splitting. Two criteria can be used:
                     50: *          SPLTOL<0 : criterion based on absolute off-diagonal value
                     51: *          SPLTOL>0 : criterion that preserves relative accuracy
                     52: *
                     53: *  TNRM (input) DOUBLE PRECISION
                     54: *          The norm of the matrix.
                     55: *
                     56: *  NSPLIT  (output) INTEGER
                     57: *          The number of blocks T splits into. 1 <= NSPLIT <= N.
                     58: *
                     59: *  ISPLIT  (output) INTEGER array, dimension (N)
                     60: *          The splitting points, at which T breaks up into blocks.
                     61: *          The first block consists of rows/columns 1 to ISPLIT(1),
                     62: *          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
                     63: *          etc., and the NSPLIT-th consists of rows/columns
                     64: *          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
                     65: *
                     66: *
                     67: *  INFO    (output) INTEGER
                     68: *          = 0:  successful exit
                     69: *
                     70: *  Further Details
                     71: *  ===============
                     72: *
                     73: *  Based on contributions by
                     74: *     Beresford Parlett, University of California, Berkeley, USA
                     75: *     Jim Demmel, University of California, Berkeley, USA
                     76: *     Inderjit Dhillon, University of Texas, Austin, USA
                     77: *     Osni Marques, LBNL/NERSC, USA
                     78: *     Christof Voemel, University of California, Berkeley, USA
                     79: *
                     80: *  =====================================================================
                     81: *
                     82: *     .. Parameters ..
                     83:       DOUBLE PRECISION   ZERO
                     84:       PARAMETER          ( ZERO = 0.0D0 )
                     85: *     ..
                     86: *     .. Local Scalars ..
                     87:       INTEGER            I
                     88:       DOUBLE PRECISION   EABS, TMP1
                     89: 
                     90: *     ..
                     91: *     .. Intrinsic Functions ..
                     92:       INTRINSIC          ABS
                     93: *     ..
                     94: *     .. Executable Statements ..
                     95: *
                     96:       INFO = 0
                     97: 
                     98: *     Compute splitting points
                     99:       NSPLIT = 1
                    100:       IF(SPLTOL.LT.ZERO) THEN
                    101: *        Criterion based on absolute off-diagonal value
                    102:          TMP1 = ABS(SPLTOL)* TNRM
                    103:          DO 9 I = 1, N-1
                    104:             EABS = ABS( E(I) )
                    105:             IF( EABS .LE. TMP1) THEN
                    106:                E(I) = ZERO
                    107:                E2(I) = ZERO
                    108:                ISPLIT( NSPLIT ) = I
                    109:                NSPLIT = NSPLIT + 1
                    110:             END IF
                    111:  9       CONTINUE
                    112:       ELSE
                    113: *        Criterion that guarantees relative accuracy
                    114:          DO 10 I = 1, N-1
                    115:             EABS = ABS( E(I) )
                    116:             IF( EABS .LE. SPLTOL * SQRT(ABS(D(I)))*SQRT(ABS(D(I+1))) )
                    117:      $      THEN
                    118:                E(I) = ZERO
                    119:                E2(I) = ZERO
                    120:                ISPLIT( NSPLIT ) = I
                    121:                NSPLIT = NSPLIT + 1
                    122:             END IF
                    123:  10      CONTINUE
                    124:       ENDIF
                    125:       ISPLIT( NSPLIT ) = N
                    126: 
                    127:       RETURN
                    128: *
                    129: *     End of DLARRA
                    130: *
                    131:       END