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 209