| 328 | |
| 329 | |
| 330 | LTriMatrix::value_type LdLfactor (LTriMatrix& M, std::size_t n) |
| 331 | /* In place modified lower triangular Cholesky factor of a |
| 332 | * Positive definite or semi-definite matrix M |
| 333 | * Reference: A+G p.218 Lower Cholesky algorithm modified for LdL' |
| 334 | * |
| 335 | * Input: M, n=last std::size_t to be included in factorisation |
| 336 | * Output: M as LdL' factor |
| 337 | * strict_lower_triangle(M) = strict_lower_triangle(L) |
| 338 | * diagonal(M) = d |
| 339 | * Return: |
| 340 | * reciprocal condition number, -1 if negative, 0 if semi-definite (including zero) |
| 341 | * ISSUE: This could change to be equivalent to UdUfactor_varient2 |
| 342 | */ |
| 343 | { |
| 344 | std::size_t i,j,k; |
| 345 | LTriMatrix::value_type e, d; |
| 346 | |
| 347 | for (j = 0; j < n; ++j) |
| 348 | { |
| 349 | d = M(j,j); |
| 350 | |
| 351 | // Diagonal element |
| 352 | if (d > 0) |
| 353 | { |
| 354 | // Positive definite |
| 355 | d = 1 / d; |
| 356 | |
| 357 | for (i = j+1; i < n; ++i) |
| 358 | { |
| 359 | e = M(i,j); |
| 360 | M(i,j) = d*e; |
| 361 | for (k = i; k < n; ++k) |
| 362 | { |
| 363 | LTriMatrix::Row Mk(M,k); |
| 364 | Mk[i] -= e*Mk[j]; |
| 365 | } |
| 366 | } |
| 367 | } |
| 368 | else if (d == 0) |
| 369 | { |
| 370 | // Possibly semi-definite, check not negative |
| 371 | for (i = j+1; i < n; ++i) |
| 372 | { |
| 373 | if (M(i,j) != 0) |
| 374 | goto Negative; |
| 375 | } |
| 376 | } |
| 377 | else |
| 378 | { |
| 379 | // Negative |
| 380 | goto Negative; |
| 381 | } |
| 382 | } |
| 383 | |
| 384 | // Estimate the reciprocal condition number |
| 385 | return rcond_internal (diag(M,n)); |
| 386 | |
| 387 | Negative: |
no test coverage detected