| 80 | |
| 81 | template <class V> |
| 82 | inline typename V::value_type rcond_ignore_infinity_internal (const V& D) |
| 83 | /* Estimate the reciprocal condition number of a Diagonal Matrix for inversion. |
| 84 | * Same as rcond_internal except that elements are infinity are ignored |
| 85 | * when determining the maximum element. |
| 86 | */ |
| 87 | { |
| 88 | // Special case an empty matrix |
| 89 | const std::size_t n = D.size(); |
| 90 | if (n == 0) |
| 91 | return 0; |
| 92 | |
| 93 | Vec::value_type rcond, mind = D[0], maxd = 0; |
| 94 | |
| 95 | for (std::size_t i = 0; i < n; ++i) { |
| 96 | Vec::value_type d = D[i]; |
| 97 | if (d != d) // NaN |
| 98 | return -1; |
| 99 | if (d < mind) mind = d; |
| 100 | if (d > maxd && 1/d != 0) // ignore infinity for maxd |
| 101 | maxd = d; |
| 102 | } |
| 103 | |
| 104 | if (mind < 0) // matrix is negative |
| 105 | return -1; |
| 106 | // ISSUE mind may still be -0, this is progated into rcond |
| 107 | if (maxd == 0) // singular due to maxd == zero (elements all zero or infinity) |
| 108 | return 0; |
| 109 | assert (mind <= maxd); // check sanity |
| 110 | |
| 111 | rcond = mind / maxd; // rcond from min/max norm |
| 112 | // CRITICAL CHECK requires NaN != NaN |
| 113 | if (rcond != rcond) // NaN, singular due to (mind == maxd) == infinity |
| 114 | rcond = 0; |
| 115 | assert (rcond <= 1); |
| 116 | return rcond; |
| 117 | } |
| 118 | |
| 119 | Vec::value_type UdUrcond (const Vec& d) |
| 120 | /* Estimate the reciprocal condition number for inversion of the original PSD |
no test coverage detected