Return to dlarrr.f CVS log

Up to [local] / rpl / lapack / lapack

En route pour la 4.0.15 !

    1:       SUBROUTINE DLARRR( N, D, E, INFO )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            N, INFO
   10: *     ..
   11: *     .. Array Arguments ..
   12:       DOUBLE PRECISION   D( * ), E( * )
   13: *     ..
   14: *
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  Perform tests to decide whether the symmetric tridiagonal matrix T
   20: *  warrants expensive computations which guarantee high relative accuracy
   21: *  in the eigenvalues.
   22: *
   23: *  Arguments
   24: *  =========
   25: *
   26: *  N       (input) INTEGER
   27: *          The order of the matrix. N > 0.
   28: *
   29: *  D       (input) DOUBLE PRECISION array, dimension (N)
   30: *          The N diagonal elements of the tridiagonal matrix T.
   31: *
   32: *  E       (input/output) DOUBLE PRECISION array, dimension (N)
   33: *          On entry, the first (N-1) entries contain the subdiagonal
   34: *          elements of the tridiagonal matrix T; E(N) is set to ZERO.
   35: *
   36: *  INFO    (output) INTEGER
   37: *          INFO = 0(default) : the matrix warrants computations preserving
   38: *                              relative accuracy.
   39: *          INFO = 1          : the matrix warrants computations guaranteeing
   40: *                              only absolute accuracy.
   41: *
   42: *  Further Details
   43: *  ===============
   44: *
   45: *  Based on contributions by
   46: *     Beresford Parlett, University of California, Berkeley, USA
   47: *     Jim Demmel, University of California, Berkeley, USA
   48: *     Inderjit Dhillon, University of Texas, Austin, USA
   49: *     Osni Marques, LBNL/NERSC, USA
   50: *     Christof Voemel, University of California, Berkeley, USA
   51: *
   52: *  =====================================================================
   53: *
   54: *     .. Parameters ..
   55:       DOUBLE PRECISION   ZERO, RELCOND
   56:       PARAMETER          ( ZERO = 0.0D0,
   57:      $                     RELCOND = 0.999D0 )
   58: *     ..
   59: *     .. Local Scalars ..
   60:       INTEGER            I
   61:       LOGICAL            YESREL
   62:       DOUBLE PRECISION   EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
   63:      $          OFFDIG, OFFDIG2
   64: 
   65: *     ..
   66: *     .. External Functions ..
   67:       DOUBLE PRECISION   DLAMCH
   68:       EXTERNAL           DLAMCH
   69: *     ..
   70: *     .. Intrinsic Functions ..
   71:       INTRINSIC          ABS
   72: *     ..
   73: *     .. Executable Statements ..
   74: *
   75: *     As a default, do NOT go for relative-accuracy preserving computations.
   76:       INFO = 1
   77: 
   78:       SAFMIN = DLAMCH( 'Safe minimum' )
   79:       EPS = DLAMCH( 'Precision' )
   80:       SMLNUM = SAFMIN / EPS
   81:       RMIN = SQRT( SMLNUM )
   82: 
   83: *     Tests for relative accuracy
   84: *
   85: *     Test for scaled diagonal dominance
   86: *     Scale the diagonal entries to one and check whether the sum of the
   87: *     off-diagonals is less than one
   88: *
   89: *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
   90: *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
   91: *     accuracy is promised.  In the notation of the code fragment below,
   92: *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
   93: *     We don't think it is worth going into "sdd mode" unless the relative
   94: *     condition number is reasonable, not 1/macheps.
   95: *     The threshold should be compatible with other thresholds used in the
   96: *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
   97: *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
   98: *     instead of the current OFFDIG + OFFDIG2 < 1
   99: *
  100:       YESREL = .TRUE.
  101:       OFFDIG = ZERO
  102:       TMP = SQRT(ABS(D(1)))
  103:       IF (TMP.LT.RMIN) YESREL = .FALSE.
  104:       IF(.NOT.YESREL) GOTO 11
  105:       DO 10 I = 2, N
  106:          TMP2 = SQRT(ABS(D(I)))
  107:          IF (TMP2.LT.RMIN) YESREL = .FALSE.
  108:          IF(.NOT.YESREL) GOTO 11
  109:          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
  110:          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
  111:          IF(.NOT.YESREL) GOTO 11
  112:          TMP = TMP2
  113:          OFFDIG = OFFDIG2
  114:  10   CONTINUE
  115:  11   CONTINUE
  116: 
  117:       IF( YESREL ) THEN
  118:          INFO = 0
  119:          RETURN
  120:       ELSE
  121:       ENDIF
  122: *
  123: 
  124: *
  125: *     *** MORE TO BE IMPLEMENTED ***
  126: *
  127: 
  128: *
  129: *     Test if the lower bidiagonal matrix L from T = L D L^T
  130: *     (zero shift facto) is well conditioned
  131: *
  132: 
  133: *
  134: *     Test if the upper bidiagonal matrix U from T = U D U^T
  135: *     (zero shift facto) is well conditioned.
  136: *     In this case, the matrix needs to be flipped and, at the end
  137: *     of the eigenvector computation, the flip needs to be applied
  138: *     to the computed eigenvectors (and the support)
  139: *
  140: 
  141: *
  142:       RETURN
  143: *
  144: *     END OF DLARRR
  145: *
  146:       END