| 390 | |
| 391 | |
| 392 | UTriMatrix::value_type UCfactor (UTriMatrix& M, std::size_t n) |
| 393 | /* In place upper triangular Cholesky factor of a |
| 394 | * Positive definite or semi-definite matrix M |
| 395 | * Reference: A+G p.218 |
| 396 | * Strict lower triangle of M is ignored in computation |
| 397 | * |
| 398 | * Input: M, n=last std::size_t to be included in factorisation |
| 399 | * Output: M as UC*UC' factor |
| 400 | * upper_triangle(M) = UC |
| 401 | * Return: |
| 402 | * reciprocal condition number, -1 if negative, 0 if semi-definite (including zero) |
| 403 | */ |
| 404 | { |
| 405 | std::size_t i,j,k; |
| 406 | UTriMatrix::value_type e, d; |
| 407 | |
| 408 | if (n > 0) |
| 409 | { |
| 410 | j = n-1; |
| 411 | do { |
| 412 | d = M(j,j); |
| 413 | |
| 414 | // Diagonal element |
| 415 | if (d > 0) |
| 416 | { |
| 417 | // Positive definite |
| 418 | d = std::sqrt(d); |
| 419 | M(j,j) = d; |
| 420 | d = 1 / d; |
| 421 | |
| 422 | for (i = 0; i < j; ++i) |
| 423 | { |
| 424 | e = d*M(i,j); |
| 425 | M(i,j) = e; |
| 426 | for (k = 0; k <= i; ++k) |
| 427 | { |
| 428 | UTriMatrix::Row Mk(M,k); |
| 429 | Mk[i] -= e*Mk[j]; |
| 430 | } |
| 431 | } |
| 432 | } |
| 433 | else if (d == 0) |
| 434 | { |
| 435 | // Possibly semi-definite, check not negative |
| 436 | for (i = 0; i < j; ++i) |
| 437 | { |
| 438 | if (M(i,j) != 0) |
| 439 | goto Negative; |
| 440 | } |
| 441 | } |
| 442 | else |
| 443 | { |
| 444 | // Negative |
| 445 | goto Negative; |
| 446 | } |
| 447 | } while (j-- > 0); |
| 448 | } |
| 449 |
no test coverage detected