SRC/dlarrr.f

   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 *> \date September 2012
  82 *
  83 *> \ingroup auxOTHERauxiliary
  84 *
  85 *> \par Contributors:
  86 *  ==================
  87 *>
  88 *> Beresford Parlett, University of California, Berkeley, USA \n
  89 *> Jim Demmel, University of California, Berkeley, USA \n
  90 *> Inderjit Dhillon, University of Texas, Austin, USA \n
  91 *> Osni Marques, LBNL/NERSC, USA \n
  92 *> Christof Voemel, University of California, Berkeley, USA
  93 *
  94 *  =====================================================================
  95       SUBROUTINE DLARRR( N, D, E, INFO )
  96 *
  97 *  -- LAPACK auxiliary routine (version 3.4.2) --
  98 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  99 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 100 *     September 2012
 101 *
 102 *     .. Scalar Arguments ..
 103       INTEGER            N, INFO
 104 *     ..
 105 *     .. Array Arguments ..
 106       DOUBLE PRECISION   D( * ), E( * )
 107 *     ..
 108 *
 109 *
 110 *  =====================================================================
 111 *
 112 *     .. Parameters ..
 113       DOUBLE PRECISION   ZERO, RELCOND
 114       PARAMETER          ( ZERO = 0.0D0,
 115      $                     RELCOND = 0.999D0 )
 116 *     ..
 117 *     .. Local Scalars ..
 118       INTEGER            I
 119       LOGICAL            YESREL
 120       DOUBLE PRECISION   EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
 121      $          OFFDIG, OFFDIG2
 122
 123 *     ..
 124 *     .. External Functions ..
 125       DOUBLE PRECISION   DLAMCH
 126       EXTERNAL           DLAMCH
 127 *     ..
 128 *     .. Intrinsic Functions ..
 129       INTRINSIC          ABS
 130 *     ..
 131 *     .. Executable Statements ..
 132 *
 133 *     As a default, do NOT go for relative-accuracy preserving computations.
 134       INFO = 1
 135
 136       SAFMIN = DLAMCH( 'Safe minimum' )
 137       EPS = DLAMCH( 'Precision' )
 138       SMLNUM = SAFMIN / EPS
 139       RMIN = SQRT( SMLNUM )
 140
 141 *     Tests for relative accuracy
 142 *
 143 *     Test for scaled diagonal dominance
 144 *     Scale the diagonal entries to one and check whether the sum of the
 145 *     off-diagonals is less than one
 146 *
 147 *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
 148 *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
 149 *     accuracy is promised.  In the notation of the code fragment below,
 150 *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
 151 *     We don't think it is worth going into "sdd mode" unless the relative
 152 *     condition number is reasonable, not 1/macheps.
 153 *     The threshold should be compatible with other thresholds used in the
 154 *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
 155 *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
 156 *     instead of the current OFFDIG + OFFDIG2 < 1
 157 *
 158       YESREL = .TRUE.
 159       OFFDIG = ZERO
 160       TMP = SQRT(ABS(D(1)))
 161       IF (TMP.LT.RMIN) YESREL = .FALSE.
 162       IF(.NOT.YESREL) GOTO 11
 163       DO 10 I = 2, N
 164          TMP2 = SQRT(ABS(D(I)))
 165          IF (TMP2.LT.RMIN) YESREL = .FALSE.
 166          IF(.NOT.YESREL) GOTO 11
 167          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
 168          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
 169          IF(.NOT.YESREL) GOTO 11
 170          TMP = TMP2
 171          OFFDIG = OFFDIG2
 172  10   CONTINUE
 173  11   CONTINUE
 174
 175       IF( YESREL ) THEN
 176          INFO = 0
 177          RETURN
 178       ELSE
 179       ENDIF
 180 *
 181
 182 *
 183 *     *** MORE TO BE IMPLEMENTED ***
 184 *
 185
 186 *
 187 *     Test if the lower bidiagonal matrix L from T = L D L^T
 188 *     (zero shift facto) is well conditioned
 189 *
 190
 191 *
 192 *     Test if the upper bidiagonal matrix U from T = U D U^T
 193 *     (zero shift facto) is well conditioned.
 194 *     In this case, the matrix needs to be flipped and, at the end
 195 *     of the eigenvector computation, the flip needs to be applied
 196 *     to the computed eigenvectors (and the support)
 197 *
 198
 199 *
 200       RETURN
 201 *
 202 *     END OF DLARRR
 203 *
 204       END