Return to dlarra.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    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