Function UdUfactor_variant1

bayes/src/UdU.cpp:193–257

RowMatrix::value_type UdUfactor_variant1 (RowMatrix& M, std::size_t n)
/* In place Modified upper triangular Cholesky factor of a
 *  Positive definite or semi-definite matrix M
 * Reference: A+G p.218 Upper Cholesky algorithm modified for UdU'
 *  Numerical stability may not be as good as M(k,i) is updated from previous results
 *  Algorithm has poor locality of reference and avoided for large matrices
 *  Infinity values on the diagonal can be factorised
 *
 * Strict lower triangle of M is ignored in computation
 *
 * Input: M, n=last std::size_t to be included in factorisation
 * Output: M as UdU' factor
 *    strict_upper_triangle(M) = strict_upper_triangle(U)
 *    diagonal(M) = d
 *    strict_lower_triangle(M) is unmodified
 * Return:
 *    reciprocal condition number, -1 if negative, 0 if semi-definite (including zero)
 */
{
	std::size_t i,j,k;
	RowMatrix::value_type e, d;

	if (n > 0)
	{
		j = n-1;
		do {
			d = M(j,j);

			// Diagonal element
			if (d > 0)
			{	// Positive definite
				d = 1 / d;

				for (i = 0; i < j; ++i)
				{
					e = M(i,j);
					M(i,j) = d*e;
					for (k = 0; k <= i; ++k)
					{
						RowMatrix::Row Mk(M,k);
						Mk[i] -= e*Mk[j];
					}
				}
			}
			else if (d == 0)
			{	// Possibly semi-definite, check not negative
				for (i = 0; i < j; ++i)
				{
					if (M(i,j) != 0)
						goto Negative;
				}
			}
			else
			{	// Negative
				goto Negative;
			}
		} while (j-- > 0);
	}