dlarrr.c 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177
  1. /* dlarrr.f -- translated by f2c (version 20061008).
  2. You must link the resulting object file with libf2c:
  3. on Microsoft Windows system, link with libf2c.lib;
  4. on Linux or Unix systems, link with .../path/to/libf2c.a -lm
  5. or, if you install libf2c.a in a standard place, with -lf2c -lm
  6. -- in that order, at the end of the command line, as in
  7. cc *.o -lf2c -lm
  8. Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
  9. http://www.netlib.org/f2c/libf2c.zip
  10. */
  11. #include "f2c.h"
  12. #include "blaswrap.h"
  13. /* Subroutine */ int dlarrr_(integer *n, doublereal *d__, doublereal *e,
  14. integer *info)
  15. {
  16. /* System generated locals */
  17. integer i__1;
  18. doublereal d__1;
  19. /* Builtin functions */
  20. double sqrt(doublereal);
  21. /* Local variables */
  22. integer i__;
  23. doublereal eps, tmp, tmp2, rmin;
  24. extern doublereal dlamch_(char *);
  25. doublereal offdig, safmin;
  26. logical yesrel;
  27. doublereal smlnum, offdig2;
  28. /* -- LAPACK auxiliary routine (version 3.2) -- */
  29. /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  30. /* November 2006 */
  31. /* .. Scalar Arguments .. */
  32. /* .. */
  33. /* .. Array Arguments .. */
  34. /* .. */
  35. /* Purpose */
  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. /* Arguments */
  41. /* ========= */
  42. /* N (input) INTEGER */
  43. /* The order of the matrix. N > 0. */
  44. /* D (input) DOUBLE PRECISION array, dimension (N) */
  45. /* The N diagonal elements of the tridiagonal matrix T. */
  46. /* E (input/output) DOUBLE PRECISION array, dimension (N) */
  47. /* On entry, the first (N-1) entries contain the subdiagonal */
  48. /* elements of the tridiagonal matrix T; E(N) is set to ZERO. */
  49. /* INFO (output) INTEGER */
  50. /* INFO = 0(default) : the matrix warrants computations preserving */
  51. /* relative accuracy. */
  52. /* INFO = 1 : the matrix warrants computations guaranteeing */
  53. /* only absolute accuracy. */
  54. /* Further Details */
  55. /* =============== */
  56. /* Based on contributions by */
  57. /* Beresford Parlett, University of California, Berkeley, USA */
  58. /* Jim Demmel, University of California, Berkeley, USA */
  59. /* Inderjit Dhillon, University of Texas, Austin, USA */
  60. /* Osni Marques, LBNL/NERSC, USA */
  61. /* Christof Voemel, University of California, Berkeley, USA */
  62. /* ===================================================================== */
  63. /* .. Parameters .. */
  64. /* .. */
  65. /* .. Local Scalars .. */
  66. /* .. */
  67. /* .. External Functions .. */
  68. /* .. */
  69. /* .. Intrinsic Functions .. */
  70. /* .. */
  71. /* .. Executable Statements .. */
  72. /* As a default, do NOT go for relative-accuracy preserving computations. */
  73. /* Parameter adjustments */
  74. --e;
  75. --d__;
  76. /* Function Body */
  77. *info = 1;
  78. safmin = dlamch_("Safe minimum");
  79. eps = dlamch_("Precision");
  80. smlnum = safmin / eps;
  81. rmin = sqrt(smlnum);
  82. /* Tests for relative accuracy */
  83. /* Test for scaled diagonal dominance */
  84. /* Scale the diagonal entries to one and check whether the sum of the */
  85. /* off-diagonals is less than one */
  86. /* The sdd relative error bounds have a 1/(1- 2*x) factor in them, */
  87. /* x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative */
  88. /* accuracy is promised. In the notation of the code fragment below, */
  89. /* 1/(1 - (OFFDIG + OFFDIG2)) is the condition number. */
  90. /* We don't think it is worth going into "sdd mode" unless the relative */
  91. /* condition number is reasonable, not 1/macheps. */
  92. /* The threshold should be compatible with other thresholds used in the */
  93. /* code. We set OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds */
  94. /* to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000 */
  95. /* instead of the current OFFDIG + OFFDIG2 < 1 */
  96. yesrel = TRUE_;
  97. offdig = 0.;
  98. tmp = sqrt((abs(d__[1])));
  99. if (tmp < rmin) {
  100. yesrel = FALSE_;
  101. }
  102. if (! yesrel) {
  103. goto L11;
  104. }
  105. i__1 = *n;
  106. for (i__ = 2; i__ <= i__1; ++i__) {
  107. tmp2 = sqrt((d__1 = d__[i__], abs(d__1)));
  108. if (tmp2 < rmin) {
  109. yesrel = FALSE_;
  110. }
  111. if (! yesrel) {
  112. goto L11;
  113. }
  114. offdig2 = (d__1 = e[i__ - 1], abs(d__1)) / (tmp * tmp2);
  115. if (offdig + offdig2 >= .999) {
  116. yesrel = FALSE_;
  117. }
  118. if (! yesrel) {
  119. goto L11;
  120. }
  121. tmp = tmp2;
  122. offdig = offdig2;
  123. /* L10: */
  124. }
  125. L11:
  126. if (yesrel) {
  127. *info = 0;
  128. return 0;
  129. } else {
  130. }
  131. /* *** MORE TO BE IMPLEMENTED *** */
  132. /* Test if the lower bidiagonal matrix L from T = L D L^T */
  133. /* (zero shift facto) is well conditioned */
  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. return 0;
  140. /* END OF DLARRR */
  141. } /* dlarrr_ */