File:  [local] / rpl / lapack / lapack / dlarrr.f
Revision 1.19: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:57 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief \b DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DLARRR + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrr.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrr.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrr.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLARRR( N, D, E, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            N, INFO
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       DOUBLE PRECISION   D( * ), E( * )
   28: *       ..
   29: *
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> Perform tests to decide whether the symmetric tridiagonal matrix T
   38: *> warrants expensive computations which guarantee high relative accuracy
   39: *> in the eigenvalues.
   40: *> \endverbatim
   41: *
   42: *  Arguments:
   43: *  ==========
   44: *
   45: *> \param[in] N
   46: *> \verbatim
   47: *>          N is INTEGER
   48: *>          The order of the matrix. N > 0.
   49: *> \endverbatim
   50: *>
   51: *> \param[in] D
   52: *> \verbatim
   53: *>          D is DOUBLE PRECISION array, dimension (N)
   54: *>          The N diagonal elements of the tridiagonal matrix T.
   55: *> \endverbatim
   56: *>
   57: *> \param[in,out] E
   58: *> \verbatim
   59: *>          E is DOUBLE PRECISION array, dimension (N)
   60: *>          On entry, the first (N-1) entries contain the subdiagonal
   61: *>          elements of the tridiagonal matrix T; E(N) is set to ZERO.
   62: *> \endverbatim
   63: *>
   64: *> \param[out] INFO
   65: *> \verbatim
   66: *>          INFO is INTEGER
   67: *>          INFO = 0(default) : the matrix warrants computations preserving
   68: *>                              relative accuracy.
   69: *>          INFO = 1          : the matrix warrants computations guaranteeing
   70: *>                              only absolute accuracy.
   71: *> \endverbatim
   72: *
   73: *  Authors:
   74: *  ========
   75: *
   76: *> \author Univ. of Tennessee
   77: *> \author Univ. of California Berkeley
   78: *> \author Univ. of Colorado Denver
   79: *> \author NAG Ltd.
   80: *
   81: *> \ingroup OTHERauxiliary
   82: *
   83: *> \par Contributors:
   84: *  ==================
   85: *>
   86: *> Beresford Parlett, University of California, Berkeley, USA \n
   87: *> Jim Demmel, University of California, Berkeley, USA \n
   88: *> Inderjit Dhillon, University of Texas, Austin, USA \n
   89: *> Osni Marques, LBNL/NERSC, USA \n
   90: *> Christof Voemel, University of California, Berkeley, USA
   91: *
   92: *  =====================================================================
   93:       SUBROUTINE DLARRR( N, D, E, INFO )
   94: *
   95: *  -- LAPACK auxiliary routine --
   96: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   97: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   98: *
   99: *     .. Scalar Arguments ..
  100:       INTEGER            N, INFO
  101: *     ..
  102: *     .. Array Arguments ..
  103:       DOUBLE PRECISION   D( * ), E( * )
  104: *     ..
  105: *
  106: *
  107: *  =====================================================================
  108: *
  109: *     .. Parameters ..
  110:       DOUBLE PRECISION   ZERO, RELCOND
  111:       PARAMETER          ( ZERO = 0.0D0,
  112:      $                     RELCOND = 0.999D0 )
  113: *     ..
  114: *     .. Local Scalars ..
  115:       INTEGER            I
  116:       LOGICAL            YESREL
  117:       DOUBLE PRECISION   EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
  118:      $          OFFDIG, OFFDIG2
  119: 
  120: *     ..
  121: *     .. External Functions ..
  122:       DOUBLE PRECISION   DLAMCH
  123:       EXTERNAL           DLAMCH
  124: *     ..
  125: *     .. Intrinsic Functions ..
  126:       INTRINSIC          ABS
  127: *     ..
  128: *     .. Executable Statements ..
  129: *
  130: *     Quick return if possible
  131: *
  132:       IF( N.LE.0 ) THEN
  133:          INFO = 0
  134:          RETURN
  135:       END IF
  136: *
  137: *     As a default, do NOT go for relative-accuracy preserving computations.
  138:       INFO = 1
  139: 
  140:       SAFMIN = DLAMCH( 'Safe minimum' )
  141:       EPS = DLAMCH( 'Precision' )
  142:       SMLNUM = SAFMIN / EPS
  143:       RMIN = SQRT( SMLNUM )
  144: 
  145: *     Tests for relative accuracy
  146: *
  147: *     Test for scaled diagonal dominance
  148: *     Scale the diagonal entries to one and check whether the sum of the
  149: *     off-diagonals is less than one
  150: *
  151: *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
  152: *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
  153: *     accuracy is promised.  In the notation of the code fragment below,
  154: *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
  155: *     We don't think it is worth going into "sdd mode" unless the relative
  156: *     condition number is reasonable, not 1/macheps.
  157: *     The threshold should be compatible with other thresholds used in the
  158: *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
  159: *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
  160: *     instead of the current OFFDIG + OFFDIG2 < 1
  161: *
  162:       YESREL = .TRUE.
  163:       OFFDIG = ZERO
  164:       TMP = SQRT(ABS(D(1)))
  165:       IF (TMP.LT.RMIN) YESREL = .FALSE.
  166:       IF(.NOT.YESREL) GOTO 11
  167:       DO 10 I = 2, N
  168:          TMP2 = SQRT(ABS(D(I)))
  169:          IF (TMP2.LT.RMIN) YESREL = .FALSE.
  170:          IF(.NOT.YESREL) GOTO 11
  171:          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
  172:          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
  173:          IF(.NOT.YESREL) GOTO 11
  174:          TMP = TMP2
  175:          OFFDIG = OFFDIG2
  176:  10   CONTINUE
  177:  11   CONTINUE
  178: 
  179:       IF( YESREL ) THEN
  180:          INFO = 0
  181:          RETURN
  182:       ELSE
  183:       ENDIF
  184: *
  185: 
  186: *
  187: *     *** MORE TO BE IMPLEMENTED ***
  188: *
  189: 
  190: *
  191: *     Test if the lower bidiagonal matrix L from T = L D L^T
  192: *     (zero shift facto) is well conditioned
  193: *
  194: 
  195: *
  196: *     Test if the upper bidiagonal matrix U from T = U D U^T
  197: *     (zero shift facto) is well conditioned.
  198: *     In this case, the matrix needs to be flipped and, at the end
  199: *     of the eigenvector computation, the flip needs to be applied
  200: *     to the computed eigenvectors (and the support)
  201: *
  202: 
  203: *
  204:       RETURN
  205: *
  206: *     End of DLARRR
  207: *
  208:       END

CVSweb interface <joel.bertrand@systella.fr>