| 32 | |
| 33 | template <class V> |
| 34 | inline typename V::value_type rcond_internal (const V& D) |
| 35 | /* Estimate the reciprocal condition number of a Diagonal Matrix for inversion. |
| 36 | * D represents a diagonal matrix, the parameter is actually passed as a vector |
| 37 | * |
| 38 | * The Condition Number is defined from a matrix norm. |
| 39 | * Choose max element of D as the norm of the original matrix. |
| 40 | * Assume this norm for inverse matrix is min element D. |
| 41 | * Therefore rcond = min/max |
| 42 | * |
| 43 | * Note: |
| 44 | * Defined to be 0 for semi-definite and 0 for an empty matrix |
| 45 | * Defined to be 0 for max and min infinite |
| 46 | * Defined to be <0 for negative matrix (D element a value < 0) |
| 47 | * Defined to be <0 with any NaN element |
| 48 | * |
| 49 | * A negative matrix may be due to errors in the original matrix resulting in |
| 50 | * a factorisation producing special values in D (e.g. -infinity,NaN etc) |
| 51 | * By definition rcond <= 1 as min<=max |
| 52 | */ |
| 53 | { |
| 54 | // Special case an empty matrix |
| 55 | const std::size_t n = D.size(); |
| 56 | if (n == 0) |
| 57 | return 0; |
| 58 | |
| 59 | Vec::value_type rcond, mind = D[0], maxd = 0; |
| 60 | |
| 61 | for (std::size_t i = 0; i < n; ++i) { |
| 62 | Vec::value_type d = D[i]; |
| 63 | if (d != d) // NaN |
| 64 | return -1; |
| 65 | if (d < mind) mind = d; |
| 66 | if (d > maxd) maxd = d; |
| 67 | } |
| 68 | |
| 69 | if (mind < 0) // matrix is negative |
| 70 | return -1; |
| 71 | // ISSUE mind may still be -0, this is progated into rcond |
| 72 | assert (mind <= maxd); // check sanity |
| 73 | |
| 74 | rcond = mind / maxd; // rcond from min/max norm |
| 75 | if (rcond != rcond) // NaN, singular due to (mind == maxd) == (zero or infinity) |
| 76 | rcond = 0; |
| 77 | assert (rcond <= 1); |
| 78 | return rcond; |
| 79 | } |
| 80 | |
| 81 | template <class V> |
| 82 | inline typename V::value_type rcond_ignore_infinity_internal (const V& D) |
no test coverage detected