| 258 | |
| 259 | |
| 260 | RowMatrix::value_type UdUfactor_variant2 (RowMatrix& M, std::size_t n) |
| 261 | /* In place modified upper triangular Cholesky factor of a |
| 262 | * Positive definite or semi-definite matrix M |
| 263 | * Reference: A+G p.219 right side of table |
| 264 | * Algorithm has good locality of reference and preferable for large matrices |
| 265 | * Infinity values on the diagonal cannot be factorised |
| 266 | * |
| 267 | * Strict lower triangle of M is ignored in computation |
| 268 | * |
| 269 | * Input: M, n=last std::size_t to be included in factorisation |
| 270 | * Output: M as UdU' factor |
| 271 | * strict_upper_triangle(M) = strict_upper_triangle(U) |
| 272 | * diagonal(M) = d |
| 273 | * strict_lower_triangle(M) is unmodified |
| 274 | * Return: |
| 275 | * reciprocal condition number, -1 if negative, 0 if semi-definite (including zero) |
| 276 | */ |
| 277 | { |
| 278 | std::size_t i,j,k; |
| 279 | RowMatrix::value_type e, d; |
| 280 | if (n > 0) |
| 281 | { |
| 282 | j = n-1; |
| 283 | do { |
| 284 | RowMatrix::Row Mj(M,j); |
| 285 | d = Mj[j]; |
| 286 | |
| 287 | // Diagonal element |
| 288 | if (d > 0) |
| 289 | { // Positive definite |
| 290 | i = j; |
| 291 | do |
| 292 | { |
| 293 | RowMatrix::Row Mi(M,i); |
| 294 | e = Mi[j]; |
| 295 | for (k = j+1; k < n; ++k) |
| 296 | { |
| 297 | e -= Mi[k]*M(k,k)*Mj[k]; |
| 298 | } |
| 299 | if (i == j) { |
| 300 | Mi[j] = d = e; // Diagonal element |
| 301 | } |
| 302 | else { |
| 303 | Mi[j] = e / d; |
| 304 | } |
| 305 | } while (i-- > 0); |
| 306 | } |
| 307 | else if (d == 0) |
| 308 | { // Possibly semi-definite, check not negative, whole row must be identically zero |
| 309 | for (k = j+1; k < n; ++k) |
| 310 | { |
| 311 | if (Mj[k] != 0) |
| 312 | goto Negative; |
| 313 | } |
| 314 | } |
| 315 | else |
| 316 | { // Negative |
| 317 | goto Negative; |
no test coverage detected