You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

dlarrr.f 6.1 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211
  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 June 2017
  82. *
  83. *> \ingroup OTHERauxiliary
  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.7.1) --
  98. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  99. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  100. * June 2017
  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. * .. External Functions ..
  124. DOUBLE PRECISION DLAMCH
  125. EXTERNAL DLAMCH
  126. * ..
  127. * .. Intrinsic Functions ..
  128. INTRINSIC ABS
  129. * ..
  130. * .. Executable Statements ..
  131. *
  132. * Quick return if possible
  133. *
  134. IF( N.LE.0 ) THEN
  135. INFO = 0
  136. RETURN
  137. END IF
  138. *
  139. * As a default, do NOT go for relative-accuracy preserving computations.
  140. INFO = 1
  141. SAFMIN = DLAMCH( 'Safe minimum' )
  142. EPS = DLAMCH( 'Precision' )
  143. SMLNUM = SAFMIN / EPS
  144. RMIN = SQRT( SMLNUM )
  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. IF( YESREL ) THEN
  179. INFO = 0
  180. RETURN
  181. ELSE
  182. ENDIF
  183. *
  184. *
  185. * *** MORE TO BE IMPLEMENTED ***
  186. *
  187. *
  188. * Test if the lower bidiagonal matrix L from T = L D L^T
  189. * (zero shift facto) is well conditioned
  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. RETURN
  200. *
  201. * END OF DLARRR
  202. *
  203. END